Stochastic Development Regression on Non-Linear Manifolds

TL;DR

This paper presents a novel regression model for data on non-linear manifolds using stochastic development of Euclidean processes, enabling analysis of complex geometric data like anatomical shapes.

Contribution

It introduces a stochastic development-based regression framework for manifold-valued data, incorporating geometric structure and maximum likelihood estimation with Laplace approximation.

Findings

Model effectively captures relationships between manifold data and Euclidean variables.

Simulation studies demonstrate accurate parameter estimation.

Application to Corpus Callosum shapes shows practical utility.

Abstract

We introduce a regression model for data on non-linear manifolds. The model describes the relation between a set of manifold valued observations, such as shapes of anatomical objects, and Euclidean explanatory variables. The approach is based on stochastic development of Euclidean diffusion processes to the manifold. Defining the data distribution as the transition distribution of the mapped stochastic process, parameters of the model, the non-linear analogue of design matrix and intercept, are found via maximum likelihood. The model is intrinsically related to the geometry encoded in the connection of the manifold. We propose an estimation procedure which applies the Laplace approximation of the likelihood function. A simulation study of the performance of the model is performed and the model is applied to a real dataset of Corpus Callosum shapes.

Equations

$$\mathbf{y}=\mathbf{a}+\mathbf{X}\mathbf{b}+\mathbf{\epsilon }.$$

$${\mathbf{y}}_t={\mathbf{a}}_t+{\mathbf{X}}_t\mathbf{b}+{\mathbf{\epsilon }}_t,\text{for}t\in [t_1,t_2].$$

$${\Gamma }_{ij}^k=\frac{1}{2}{g}^{kl}({\partial }_i{g}_{jl}+{\partial }_j{g}_{il}-{\partial }_l{g}_{ij})$$

$$d{U}_t=\sum _{i=1}^d{H}_i(U_t)\circ d{X}_t^i,$$

$$d{z}_t^i=\beta dt+\tilde{W}d{X}_t^i+d{\epsilon }_t.$$

$$y_i={\mathcal{Y}}_i^T+{\tilde{\epsilon }}_i,$$

$$L(\theta ;\mathbf{y})$$

$$L(\theta ;\mathbf{y})\approx f(y|d{\mathbf{x}}_t^o,\theta )p(d{\mathbf{x}}_t^o|\theta )(2\pi {)}^{\frac{m{n}_s}{2}}|\Sigma {|}^{\frac{1}{2}}{n}^{-\frac{m{n}_s}{2}},$$

$$\tilde{W}$$

$${\hat{W}}_{20}=\left(\begin{array}{cc}0.206& 0.136\\ 0.147& 0.322\end{array}\right),{\hat{W}}_{60}=\left(\begin{array}{cc}0.22& 0.11\\ 0.11& 0.21\end{array}\right),{\hat{W}}_{100}=\left(\begin{array}{cc}0.205& 0.104\\ 0.115& 0.214\end{array}\right)$$

Full text

1 Introduction

A main focus in computational anatomy is to study the shape of anatomical objects. Performing statistical analysis of anatomical objects is however challenging due to the non-linear nature of shape spaces. The established statistical theory for Euclidean data does not directly allow us to answer questions like: How does a treatment affect the deformation of an organ? or: Is it possible to categorize sick and healthy patients based on the shape of the subject’s organs?

Shape spaces are typically non-linear and often equipped with manifold structure. Examples of manifold-valued shape data include landmarks, curves, surfaces, and images with warp variation. The lack of vector space structure for manifold-valued data implies that addition and scalar multiplication are not defined. Several concepts in statistics rely on addition and scalar multiplication, these including mean value, variance, and regression models. Hence, in order to make inference on manifold-valued data, generalization of Euclidean statistical theory is necessary.

This paper focuses on generalization of regression models to manifolds. The aim is to model the relation between Euclidean explanatory variables and a manifold-valued response. The regression model has, as an example, applications in computational anatomy [23]. The proposed model can for example be used to analyze how age affects the shape of Corpus Callosum [7].

Several approaches have previously been proposed for defining normal distributions on manifolds [14, 19]. In [19], the distribution is defined based on Brownian motions in $\mathbb{R}^{m}$ and the fact that normal distributions on $\mathbb{R}^{m}$ can be defined as transition distributions of Brownian motions. The normal distribution on the manifold is then defined as the transition distribution of the stochastic development of the Euclidean Brownian motion [10]. The proposed regression model will be defined in a similar manner. The construction can be considered intrinsic as it only depends on the connection of the manifold, e.g. the Levi-Civita connection of a Riemannian manifold. It does not rely on linearization of the manifold, and it naturally includes the effect of curvature in the mapping of the stochastic processes.

In Euclidean linear regression, the relation between explanatory variables, $\boldsymbol{X}$, and a response variable, $\boldsymbol{y}$, is modeled by an affine function of $\boldsymbol{X}$,

[TABLE]

Due to the lack of vector space structure, alternatives for modeling relations between the given variables, $\boldsymbol{X}$ and $\boldsymbol{y}$, are needed in the non-linear situation. Several ideas have previously been introduced and a selection of these will be described in Section 2.

In this paper, the regression model is considered as a transported linear regression defined in $\mathbb{R}^{m}$. This approach is inspired by the transport of normal distributions defined in [19]. Notice that the linear regression model $(\ref{ereg})$ can be generalized to situations in which several observations are observed over time,

[TABLE]

Our approach suggests to define the regression model by transportation of stochastic processes, $Z_{t}=\boldsymbol{a}_{t}+\boldsymbol{X}_{t}\boldsymbol{b}+\boldsymbol{\varepsilon}_{t}$, in $\mathbb{R}^{m}$ on to the manifold in order to obtain the relation to the response variable, $\boldsymbol{y}$ (see Figure 1).

The paper will be structured as follows. In Section 2, we give a discussion on previous methods developed for regression on manifolds. Section 3 presents a short description of development of stochastic paths from a Euclidean space to the manifold. Section 4 introduces the proposed model, followed by a description of the estimation procedure in Section 5. In Section 6 and 7, illustrative examples are considered for the application and performance of the model. The paper is ended by a discussion of the defined model in Section 8.

2 Background

Multiple approaches have been proposed for generalizing regression models to non-linear manifolds. The methods consider the regression problem in different situations. In this paper we will consider the case of Euclidean exaplanatory variables and a manifold-valued response. There have been several works describing regression models for manifold-valued data in other situations [4, 1, 13, 21].

Regression models for describing the relation between a manifold-valued response and Euclidean explanatory variables have also previously been introduced. Examples include [12] in which an extrinsic regression model is introduced, and [15], which defines an intrinsic regression model where the parameter vector is estimated by minimizing the total sum of squares based on the Riemannian manifold distance. Another example is the geodesic regression model introduced in [7], which is a generalization of the linear regression model in Euclidean spaces. The relation is here modeled by a geodesic described by an initial velocity dependent on an explanatory variable and a starting point on the manifold.

In this paper, we will take a different view on how to relate the response and explanatory variables. Instead of considering the relation as being modeled by geodesics on the manifold as in [7], we will describe the relation by stochastic paths transported from the space of explanatory variables to the manifold. By defining the regression model using stochastic paths, we are able to model non-geodesic relations, incorporate several explanatory variables, and consider random effects in the model. Non-geodesic relations have been considered by others before. An example is [17] in which the geodesic regression model from [7] is generalized in order to model more complex shape changes. The regression function is in this case fitted by piecewise cubic splines that describes the variation of one explanatory variable. In [9], a regression model is introduced, in which the non-geodesic relation is obtained by time-warping. Others have proposed to model the non-geodesic relation by either a generalized polynomial regression model or by non-linear kernel-based regression [8, 25, 3, 2, 6]. On the contrary, [16] introduces the Hierarchical Geodesic Model which are able to consider several explanatory variables including random variables, but assumes nested observations and does only consider geodesic relations. A regression model, which incorporates both a non-geodesic relation and several explanatory variables, is proposed in [5]. This work defines an intrinsic regression model on Riemannian symmetric spaces, in which the regression function is obtained by minimizing the conditional mean of residuals defined by the log-map.

In addition to describing the proposed model, we perform estimation of model parameters by maximum likelihood using the transition density on the manifold. The model does not linearize the manifold as in many of the local regression models, but instead take into account the curvature of the manifold at each point as encoded in the connection through the mapping of the stochastic process.

3 Stochastic Development

In this section we give a brief description of stochastic development of curves in $\mathbb{R}^{m}$ to the manifold. The reader is referred to [10, 20, 18] for a deeper description of this concept.

Let $\mathcal{M}$ be a $d$-dimensional manifold provided with a connection $\nabla$ and metric $g$. The connection is necessary for transportation of tangent vectors along curves on the manifold. A frequently used connection is the Levi-Civita connection coming from a Riemannian structure on $\mathcal{M}$. Let $\partial_{i}$ for $i=1,\ldots,d$ denote a coordinate frame on $\mathcal{M}$ and let $dx^{i}$ be the corresponding dual frame. A connection $\nabla$ is given in terms of its Christoffel symbols defined by $\nabla_{\partial_{i}}\partial_{j}=\Gamma_{ij}^{k}\partial_{k}$. For the Levi-Civita connection, the Christoffel symbols are given by

[TABLE]

in which $g_{ij}$ is the components of $g$ in the coordinate basis, i.e. $g=g_{ij}dx^{i}dx^{j}$, and $g^{ij}$ is the inverse components.

Consider the frame bundle $\mathcal{FM}$ being the set of tuples $(y,\nu)$ in which $y\in\mathcal{M}$ and $\nu$ is a frame for the tangent space $T_{y}\mathcal{M}$. Let $\pi\colon\mathcal{FM}\to\mathcal{M}$ be the projection map given by $\pi(y,\nu)=y$ for $(y,\nu)\in\mathcal{FM}$. A smooth curve $U_{t}$ on $\mathcal{FM}$ is a smooth selection of frames, i.e. for every $t\in I$, $U_{t}=(y_{t},\nu_{t})$ in which $\nu_{t}\colon\mathbb{R}^{d}\to T_{\pi(U_{t})}\mathcal{M}$ is a frame.

Given a connection $\nabla$, the tangent space of the frame bundle, $T\mathcal{FM}$, splits into a horizontal and a vertical part, $T\mathcal{FM}=H\mathcal{FM}\oplus V\mathcal{FM}$. The horizontal subspace explains infinitesimal changes of the base point on the manifold. On the other hand, tangent vectors in $V\mathcal{FM}$ describe changes of the frame $\nu$ keeping the base point fixed. Given a tangent vector $v\in T_{y}\mathcal{M}$ and a frame $\nu$, a vector in $H_{(y,\nu)}\mathcal{FM}$ can be defined by horizontal lift. The horizontal lift of a tangent vector $v$ is the unique horizontal vector $w\in H_{(y,\nu)}\mathcal{FM}$, satisfying $\pi_{\star}w=v$, where $\pi_{\star}\colon H_{(y,\nu)}\mathcal{FM}\to T_{y}\mathcal{M}$ is induced by the projection $\pi$. The horizontal lift of $v$ will be denoted $h_{l}(v)$.

Consider a probability space $(\Omega,\mathcal{F},P)$ and a stochastic process $X_{t}\colon\Omega\to\mathcal{W}(\mathbb{R}^{m})$, where $\mathcal{W}(\mathbb{R}^{m})$ denotes the path space of $\mathbb{R}^{m}$. The stochastic development of $X_{t}$ to $\mathcal{FM}$ can be defined as a solution, $U_{t}$, of the Stratonovich stochastic differential equation,

[TABLE]

where $\circ$ symbolizes a Stratonovich stochastic differential equation. The vector fields $H_{1},\ldots,H_{d}$ denotes a basis for the horizontal subspace of $T\mathcal{FM}$. Given a point $u=(y,\nu)\in\mathcal{FM}$, $H_{i}$ are defined as $H_{i}(u)=h_{l}(\nu(e_{i})),\ i=1,\ldots,d$, where $e_{1},\ldots,e_{d}$ is the canonical basis for $\mathbb{R}^{d}$. A path $Y_{t}$ on the manifold $\mathcal{M}$ can then be obtained by the projection of $U_{t}$ onto $\mathcal{M}$ by the projection map $\pi$, i.e. $Y_{t}=\pi(U_{t})$.

Consider two processes $X^{1}_{t},X^{2}_{t}$ in $\mathbb{R}^{m}$, $t\in[0,T]$ for $T>0$, for which $X^{1}_{0}=X^{2}_{0}=\boldsymbol{x}_{0}$ and $X^{1}_{T}=X_{T}^{2}$. If $Y^{1}_{t}$, $Y^{2}_{t}$ denotes the stochastic development of $X^{1}_{t}$ and $X^{2}_{t}$ respectively on $\mathcal{M}$, then it does not in general hold that $Y^{1}_{T}=Y^{2}_{T}$ on $\mathcal{M}$ due to the curvature of the manifold.

4 Model

Let $\mathcal{M}$ be a $d$-dimensional manifold embedded in the ambient space $\mathbb{R}^{k}$ for some $k\geq d$ and consider a response variable $y$ in $\mathcal{M}$. Let $\nu_{y_{0}}\colon\mathbb{R}^{d}\to T_{y_{0}}\mathcal{M}$ be a frame for the tangent space at a reference point $y_{0}\in\mathcal{M}$. Assume that $\boldsymbol{y}_{1},\ldots,\boldsymbol{y}_{n}\in\mathbb{R}^{k}$ are $n$ realizations of $y\in\mathcal{M}$ and let $\boldsymbol{x}_{i}=(x^{1}_{i},\ldots,x^{m}_{i})\in\mathbb{R}^{m}$ denote the vector of explanatory variables for the $i$’th observation. Notice that the realizations of $y$ are assumed to lie in the ambient space $\mathbb{R}^{k}$ and not required to be in $\mathcal{M}$. This construction allows for observations measured with noise which are not necessarily observed as elements of $\mathcal{M}$.

The strategy of the proposed model is to define stochastic processes according to the generalized linear regression in $(\ref{treg})$ and transport these to the manifold by stochastic development. All stochastic processes are defined for $t\in[0,T]$ for a $T>0$. Consider for each observation $i$ the stochastic process $z_{t}^{i}\colon\Omega\to\mathcal{W}(\mathbb{R}^{m})$, solution to the stochastic differential equation,

[TABLE]

The first term, $\beta dt$, is a fixed drift for $\beta\in\mathbb{R}^{m}$. $\tilde{W}dX_{t}^{i}$ is the dependence of the explanatory variables with $X_{t}^{i}\colon\Omega\to\mathcal{W}(\mathbb{R}^{m})$ being a stochastic process satisfying $X_{0}^{i}(\omega)=0$ and $X_{T}^{i}(\omega)=\boldsymbol{x}_{i}$ for $\omega\in\Omega$. The matrix $\tilde{W}$ is a $m\times m$-dimensional matrix with columns relating to the basis vectors of the frame $\nu_{y_{0}}$ on $\mathcal{M}$. Consider the matrix $W$ with columns consisting of basis vectors of $\nu_{y_{0}}$. If $\mathcal{M}$ has a Riemannian metric, then $W=U\tilde{W}$, in which $U$ denotes a $d\times m$ orthonormal matrix with respect to the metric. Notice that this model can incorporate both fixed and random explanatory variables. If the $j$’th explanatory variable, $x_{i}^{j}$, is a random effect, $X_{t}^{ij}$ is modeled as a Brownian bridge, while it for fixed effects are modeled as a constant drift. The random error, $\varepsilon_{t}$, is modeled as a multidimensional Brownian motion on $\mathbb{R}^{m}$.

The $i$’th observation $y_{i}$ is modeled as a noisy endpoint of the stochastic development of $z_{t}^{i}$. If $m<d$ only a reduced frame $\tilde{\nu}_{y_{0}}$ is used for the stochastic development of $z_{t}^{i}$. The reduced frame is considered as we are only interested in the effect of frame vectors associated to the explanatory variables. The basis vectors of $\tilde{\nu}_{y_{0}}$ corresponds to the columns of $W$. Given the reference point $y_{0}\in\mathcal{M}$, define stochastic processes $Y_{t}^{i}$ as the stochastic development of $z_{t}^{i}$. Let $\mathcal{Y}^{T}_{i}\colon\Omega\to\mathcal{M}$ be a random variable following the distribution of endpoints of the stochastic development $Y_{t}^{i}$. Then

[TABLE]

where $\tilde{\varepsilon}_{i}\sim\mathcal{N}(\boldsymbol{0},\tau^{2}\mathbb{I}_{d})$ represents the random measurement error that pulls the realization, $y_{i}$, from the manifold. In Figure 2, the two steps of the model are illustrated. First, the stochastic development of $z_{t}^{i}$ are defined on the frame bundle and finally, this stochastic development is projected to the manifold.

Notice that in the case $\mathcal{M}=\mathbb{R}^{k}$ with the standard connection on $\mathbb{R}^{k}$, the proposed model reduces to the regular regression model for data in $\mathbb{R}^{k}$. Assume $y\in\mathbb{R}^{k}$ and that $X_{t}^{i}$ is a vector from [math] to $\boldsymbol{x}_{i}$. Then $\beta$ and $y_{0}$ relates to the intercept, $W$ is the matrix of regression coefficients and $\varepsilon_{t}$ and $\tilde{\varepsilon}$ the iid. random noise.

5 Estimation

The reference point $y_{0}$, the matrix $W$, the drift $\beta$, and the variance parameter $\tau^{2}$ are the parameters of the model. These parameters can be estimated in several ways. This section describes a Laplace approximation of the marginal likelihood function which are used for finding optimal parameter estimates. We could alternatively use a Monte Carlo EM based procedure using simulations of the missing data, $Y_{t}^{i}$ for $t\in[0,T]$, to optimize the complete data likelihood. This will be considered in future works.

Laplace approximation can be used to determine a linear approximation of a non-linear likelihood function [11]. Let $\theta$ denote the vector of parameters, and $d\boldsymbol{x}_{t}$ a discretization of the process $X_{t}$ at $n_{s}+1$ time-points. Hence $d\boldsymbol{x}_{t}$ is a vector of length $n\cdot m\cdot n_{s}$, in which $n_{s}$ denotes the number of time steps, $n$ the number of observations, and $m$ the number of explanatory variables. Let $f(y\lvert\theta)$ be the conditional density of the response $y\in\mathcal{M}$ given $\theta$ and $p(d\boldsymbol{x}_{t}\lvert\theta)$ the density of the discretization of $X_{t}$ given $\theta$. To find the optimal parameter vector, $\theta$, the following likelihood has to be optimized,

[TABLE]

where $h(d\boldsymbol{x}_{t})=-\frac{1}{n}\log f(y\lvert d\boldsymbol{x}_{t},\theta)-\frac{1}{n}\log p(d\boldsymbol{x}_{t}\lvert\theta)$. The Laplace approximation of $L$ is then given by

[TABLE]

in which $d\boldsymbol{x}^{o}_{t}=\text{argmax}_{d\boldsymbol{x}_{t}}\{-h(d\boldsymbol{x}_{t})\}$ and $\Sigma=\left(D^{2}h(d\boldsymbol{x}_{t})\right)^{-1}$, the inverse of the Hessian of $h(d\boldsymbol{x}_{t})$. The approximated likelihood is then optimized wrt. $\theta$ to obtain the estimated parameters. In the following simulation study, the Laplace approximation is used for parameter estimation. The code for the estimation algorithm as well as the simulation study below was implemented in Theano [22]. The code is available at https://bitbucket.org/stefansommer/theanodiffgeom.

6 Simulation Study

This section investigates properties of the model on simulated synthetic data. Two setups will be introduced, both considering landmark representations of shapes. The data are assumed to lie in a manifold defined in the LDDMM (Large Deformation Diffeomorphic Metric Mapping) framework [24].

In the LDDMM framework, deformations of shapes are modeled as smooth flows which are solutions to ordinary differential equations defined by vector fields. A point $q\in\mathcal{M}$ is a finite number of landmarks, $q=(x_{1}^{1},x_{1}^{2},\ldots,x_{n_{l}}^{1},x_{n_{l}}^{2})$. The metric on $\mathcal{M}$ is given by $g(v,w)=\sum_{i,j}^{n_{l}}vK^{-1}(x_{i},x_{j})w$, where $K^{-1}$ denotes the inverse of a kernel $K$. In this simulation study $K$ is the Gaussian kernel with standard deviation, $\sigma=0.5$. Based on this metric the Levi-Civita connection can be obtained by calculating the Christoffel symbols defined in $(\ref{eq:Chris})$.

To begin with, we consider estimation of $\tilde{W}$ and $y_{0}$ and investigate the performance of the estimation procedure. The shapes that will be considered consists of $8$ landmarks generated from the unit circle with landmarks located at $0,\frac{\pi}{4},\frac{\pi}{2},\ldots,\frac{3\pi}{2},\frac{7\pi}{4}$ radians. The center plot of Figure 3 shows the unit circle with the chosen frame for each landmark. The number of explanatory variables are set to $m=2$ and the variables are drawn from a normal distribution with mean [math] and standard deviation $2$. The other parameters are set to

[TABLE]

In Figure 3 is shown an example of simulated observations as well as the sample paths $X_{t}^{i}$. A total of 50 datasets were sampled, in which each consisted of 20 observations. For each simulated dataset, the $\tilde{W}$ matrix was estimated. Each of the estimated distrubtions for the entries of $\tilde{W}$ are shown in Figure 4. By the results, we conclude that the estimated parameters are fairly stable between the different simulations and that the true values are well centered in each distribution. For this simulation, the estimation procedure is thus able to estimate the true $\tilde{W}$ parameters that were specified in the model.

Three similar datasets, as explained above, were sampled with different number of observations, 20, 60 and 100 respectively. The matrix $\tilde{W}$ as well as the reference point $y_{0}$ were estimated for each of the three datasets. In this case, the estimated $\tilde{W}$ matrix was found to be

[TABLE]

while the estimated reference points are shown in Figure 5. By increasing the number of observations, we conclude that the estimated parameters $\tilde{W}$ and $y_{0}$ converge towards the true parameters.

In the second study, we consider the problem of estimating the frame matrix $U$. In this case, each observation consists of 3 landmarks that were generated from a setup shown in Figure 6. We only consider one explanatory variable, meaning that only one frame vector has to be estimated for each landmark. The true frame vectors for each landmark was set to a vertical unit vector. In the estimation procedure, the frame vectors were initialized with the Euclidean linear regression estimate. In Figure 6 is shown the true (red), the initial (green) and the estimated frame (blue) for each landmark. The estimation procedure converges to a good estimate of the true frame. Estimation of the initial frame was considered for different number of observations, but the estimated frame did not seem to converge for increasing number of observations. The difference in the parameter estimates might therefore be a result of either the linear approximation of the likelihood or that the optimal solution of the initial frame is not unique.

7 Data Example

We now apply the model to a real dataset consisting of landmark representations of Corpus Callosum (CC) shapes. The model is used to describe the effect of age on CC shapes. The manifold considered is the same as that introduced in Section 6, but in this case $\sigma=0.1$. Again the Levi-Civita connection is used.

A subset of the CC dataset is plotted in Figure 7. For model fitting, a dataset of $20$ CC shapes was considered with age values ranging from $22$ to $78$. The model was fitted to CC shapes represented by a subset of $20$ landmarks. We did not incorporate a drift term in the model, and only the frame and $\tilde{W}$ has been estimated. The refrence point was set to the mean shape (Figure 7) and $\tau=0.1$.

The estimated frame for the $20$ landmarks are shown in Figure 7 on top of the mean shape. The weight matrix was estimated as $\tilde{W}=-0.0002$. Given the low estimate of $\tilde{W}$ and hence a small frame matrix $W$, the result of this experiment suggests a low age effect on CC for these data.

8 Discussion

A method was proposed for modeling the relation between a manifold-valued response and Euclidean explanatory variables. The relation was modeled by transport of stochastic paths from $\mathbb{R}^{m}$ to the manifold. The stochastic paths defined on $\mathbb{R}^{m}$ was given as solutions to a stochastic differential equation with a contribution from a fixed drift, a stochastic process related to the explanatory variables, and a random noise assumed to follow a multidimensional Brownian motion. The response variable was then modeled as a noisy observation of a stochastic variable following the distribution of the endpoints of the transported process. The proposed model is intrinsic and based on a connection on the manifold without making linearization of the non-linear space. Moreover, a likelihood based estimation procedure were described using Laplace approximation of the marginal likelihood. We experimentally illustrated the model and the parameter estimation using a simulation study and a real data example.

Other procedures could be used for estimation of parameters. As an example, the Monte Carlo EM procedure could be used to optimize the complete data likelihood based on simulations of the missing data. Another example is to approximate the distribution of the response by moment matching.

An interesting problem to investigate is how to make variable selection in the model. As the contribution from the explanatory variables is defined in comparison with the frame basis vectors, one idea is to exclude those explanatory variables which corresponds to frame vectors parallel to the curve. These frame vectors will not contribute to the stochastic development and hence will not be important for explaining the relation to the response variable.

An important assumption of the manifold considered, is that the manifold is equipped with a connection. In this paper, the Levi-Civita connection was used, but several other connections could have been chosen. It would be interesting to explore how the choice of connection affects the model.

As it is possible to transport stochastic paths from a manifold to a Euclidean space, the model could be generalized to handle situations in which a Euclidean response variable is compared to manifold-valued explanatory variables. Based on such a model, one might be able to make categorization of individuals based on manifold-valued shapes.