A Geodesic Mixed Effects Model in Kendall’s Shape Space

Abstract

In many applications, geodesic hierarchical models are adequate for the study of temporal observations. We employ such a model derived for manifold-valued data to Kendall’s shape space. In particular, instead of the Sasaki metric, we adapt a functional-based metric, which increases the computational efficiency and does not require the implementation of the curvature tensor. We propose the corresponding variational time discretization of geodesics and apply the approach for the estimation of group trends and statistical testing of 3D shapes derived from an open access longitudinal imaging study on osteoarthritis.

1 Introduction

Analysis of time-dependent shape data has become increasingly important for a wide range of applications. For individual biological changes, subject-specific smooth regression models are adequate. The obtained trajectories provide the possibility to estimate data values at unobserved times and to compare trends even with unbalanced data. Due to the non-linear structure of Kendall’s shape space, common statistical tools derived for Euclidean spaces are not applicable. Within the previous years, many novel approaches have been presented for the statistical analysis of time-dependent manifold-valued data [2, 8, 11, 15]. The derived schemes benefit from a compact encoding of constraints and exhibit a superior consistency as compared to their Euclidean counterparts.

In this work, we employ the generative hierarchical approach, based on geodesic analysis, introduced in [3] and [11] with some modifications. In the first stage, inner-individual changes are modeled as geodesic trends, which in the second phase are considered as disturbances of a population-averaged geodesic trend. For an introduction and overview over longitudinal analysis based on mixed effects models we refer to [5].

We use the approach, to approximate the observed temporal shape data by geodesics in the shape space and to estimate the overall trends within the groups on this basis. Geodesic models are attractive because they have a compact representation and enable computational efficiency. Moreover inconsistencies, e.g. due to acquisition noise and reconstruction errors, are minimized via geodesic regression.

Analysis of geodesic trends requires a notion of distance that is consistent with the Riemannian metric of shape space. State-of-the-art approaches parametrize geodesics as points in the tangent bundle of the shape space [11]. While the Sasaki metric is a natural metric on the tangent bundle, its geodesic computations require time-discrete approximation schemes involving the Riemannian curvature tensor. This not only incurs high computational costs but also impacts numerical stability. We consider a novel approach that overcomes these shortcomings. To this end, we identify elements of the tangent bundle with vector fields along the geodesic trend. This provides a notion of a canonical metric that is motivated from a functional view of parameterized curves in the shape space [16]. Considering the space of the geodesics as a submanifold in the space of shape trajectories, this allows in particular the use of a naturally induced distance. The corresponding shortest path, log map and average geodesic, can be computed by variational time-discretization. Remarkably, the underlying energy function allows for fast and simple evaluation increasing computational efficiency. In particular, it neither requires curvature computation nor decomposition in horizontal and vertical components.

The main challenge in the first stage, viz. geodesic regression, is to determine parallel transport and Jacobi fields in Kendall’s shape space. While these important geometric quantities are not given by closed form expressions, efficient approximation schemes have been presented [9, 10, 12].

Using the derived metric for geodesic trends, we obtain a notion of mean, covariance, and Mahalanobis distance. This allowed us to develop a statistical hypothesis test for comparing the group-wise mean trends. Non-parametric permutation tests are applied to test for significance of estimated differences in group trends. We perform this in terms of a manifold-valued Hotelling $t^2$ statistic described in [3] by applying it to the tangent bundle. As example application we demonstrate the methodology on the long term study of incident knee osteoarthritis (OA).

This paper is organized as follows. In Sect. 2, after a short overview of Kendall’s shape space, we present main mathematical tools corresponding to the geometry of the shape space and its tangent bundle, particularly parallel transport, geodesic regression and mean geodesic. Essential for these considerations is to encode the geometry of the shape space via computations in the pre-shape space. In Sect. 3 we present the application of our approach to femur data from an epidemiological longitudinal study dealing with osteoarthritis and discuss the numerical results.

2 Kendall’s Shape Space

In the following we present a brief overview of Kendall’s shape space and its tangent bundle as well as main quantities which will be employed for geodesic analysis and statistics.

2.1 Preliminaries

For a comprehensive introduction to Kendall’s shape space and details on the subjects of this section, we refer to [7] and [12]. For the relevant tools from Riemannian geometry, we refer to [4] and [13].

Let M(m, k) denote the space of real $m\times k$ matrices endowed with its canonical scalar product given by $\langle x,y\rangle ={{\,\mathrm{trace}\,}}(xy^t)$, and the induced Frobenius norm. We call the set of k-ad of landmarks in $\mathbb {R}^m$ after removing translations and scaling the pre-shape space and identify it with $\mathcal {S}_m^k~{:=}~\{x\in M(m,k):\, \sum _{i=1}^kx_i=0,\,\Vert x\Vert =1\}$ endowed with the spherical Procrustes metric $d(x,y)=\arccos (\langle x,y\rangle )$. A shape is a pre-shape with rotations removed. More precisely, the left action of $\mathrm {SO}_m$ on $\mathcal {S}_m^k$ given by $(R,x)\mapsto Rx$ defines an equivalence relation given by $x\sim y$ if and only if $y=Rx$ for some $R\in \mathrm {SO}_m$. Kendall’s shape space is defined as $\varSigma _m^k=\mathcal {S}_m^k/\mathord {\sim }$. Now, denoting the canonical projection of $\sim$ by $\pi$, the induced distance between any two shapes $\pi (x)$ and $\pi (y)$ is given by

$$\begin{aligned} d_{\varSigma }(x,y)= \min _{R\in \mathrm {SO}_m} d(x,Ry)=\arccos \sum _{i=1}^m\lambda _i, \end{aligned}$$

where $\lambda _1\ge \cdots \ge \lambda _{m-1}\ge |\lambda _m|$ denote the pseudo-singular values of $yx^t$. Note that for simplicity of notation, we have identified shapes and their representing pre-shapes in the definition of $d_{\varSigma }$. Moreover, for $k\ge 3$, the shape space $\varSigma _1^k$ (resp. $\varSigma _2^k$) is isometric to the sphere (resp. projective space). We call $x,y\in \mathcal {S}_m^k$ well positioned and write $x{\mathop {\sim }\limits ^{\omega }}y$ if and only if $yx^t$ is symmetric and $d(x,y)=d_{\varSigma }(x,y)$. For each $x,y\in \mathcal {S}_m^k$, there exists an optimal rotation $R\in \mathrm {SO}_m$ such that $x{\mathop {\sim }\limits ^{\omega }}Ry$. Due to [7] the horizontal and vertical spaces at $x\in \mathcal {S}_m^k$ read

$$\begin{aligned} \mathrm {Hor}_x&=\{u\in M(m,k-1):\,ux^t=xu^t\text { and } \langle x,u\rangle =0\},\\ \mathrm {Ver}_x&=\{Ax:\,A+A^t=0\}. \end{aligned}$$

A smooth curve is called horizontal if and only if its tangent field is horizontal. Geodesics in the shape space are equivalence classes of horizontal geodesics. Now, let $\exp$ and $\log$ denote the exponential and log map of the pre-shape sphere. For $x{\mathop {\sim }\limits ^{\omega }}y$ the geodesic from x to y given by

$$\begin{aligned} \varPhi (x,y,t)~{:=}~\exp _x(t\log _x y)=\frac{\sin ((1-t)\varphi )}{\sin \varphi }x+\frac{\sin (t\varphi )}{\sin \varphi }y \end{aligned}$$

(1)

with $\varphi =\arccos (\langle x,y\rangle ),\,0\le t\le 1$, is horizontal. Hence $\varPhi$ realizes the minimizing geodesic from $\pi (x)$ to $\pi (y)$. We recall, that pre-shapes with rank $\ge m-1$, denoted by $\mathcal {S}$, constitute an open and dense subset of $\mathcal {S}_m^k$ and the restriction of the quotient map $\pi$ to $\mathcal {S}$ is a Riemannian submersion with respect to the metric induced by the ambient Euclidean space. Key quantities of the shape space geometry such as parallel transport, Jacobi fields and Fréchet mean can be computed by horizontal lifting to $\mathcal {S}$ (and extension to $\mathcal {S}_m^k$). We refer the reader to [12] for corresponding results.

2.2 Geodesic Regression

The first stage of the employed model is geodesic regression, which we summarize below. We recall that diameter of $\varSigma _m^k$ is $\pi /2$. Now, let $\varOmega _x$ be a neighborhood of $x\in \mathcal {S}_m^k$ with radius smaller then $\pi /4$ and such that for any $y\in \varOmega _x$, denoting the pseudo-singular values of $yx^t$ by $\lambda _i$, $\lambda _{m-1}+\lambda _m> 0$ holds. Then, due to [7, 6.6], for $y\in \varOmega _x$ the optimal rotation R with $x{\mathop {\sim }\limits ^{\omega }}Ry$ is unique. Hence the function $\varOmega _x\ni y\mapsto \omega (x,y)~{:=}~Ry$ is well-defined.

Now, consider scalars $0=t_1< t_2< \cdots < t_N=1$ and distinct pre-shapes $q_1,\cdots ,q_N$. Geodesic regression aims at finding a geodesic curve in shape space that best fits the data $\pi (q_i)$ at $t_i$ in a least-squares sense. Computationally, we employ the parametrization given by (1) to determine the corresponding horizontal geodesic $\varPhi (x^*,y^*,.)$, where $(x^*,y^*)~{:=}~\text {arg min} F$ and

$$ F(x,y)~{:=}~\sum _{i=1}^Nd^2(q_i,\omega (q_i,\varPhi (x,y,t_i,))),\, x{\mathop {\sim }\limits ^{\omega }}y, $$

with $(q_1,\omega (q_1,q_N))$ as initial guess.

We recall that the significance of the regression model can be measured by coefficient of determination denoted $R^2$. To compute it, let $F_{min}:=F(x^*,y^*)$ and denote the minimum of $G(x)~{:=}~\sum _{i=1}^N d^2(q_i,\omega (x,q_i))$ by $G_{min}$. Then $R^2=1 - \frac{F_{min}}{G_{min}}$.

2.3 Tangent Bundle and Mean Geodesic

Geodesic mixed effects models and particularly mean geodesic (group trend) require a notion of distance for the tangent bundle consistent with the Riemannian metric of the shape space. In the following, we present a brief introduction to a natural choice for such a distance provided by the Sasaki metric employed in  [11]. Then, we propose an alternative $L^2$-type approach and induced variational time-discrete geodesics.

In the sequel, $I:=[0,1]$. Let (M, g) be an $\mu$-dimensional Riemannian manifold and $\delta$ a Riemannian metric on the tangent bundle TM with the canonical projection $\tau$. Identifying a geodesic $\gamma$ with $(\gamma (0),\dot{\gamma }(0))$, mean of geodesics is determined by $\delta$. A prominent natural choice for $\delta$ is the Sasaki metric. It is uniquely determined by the following properties (cf. [14]): (a) $\tau$ becomes a Riemannian submersion. (b) The restriction of $\delta$ to any tangent space coincides with g. (c) Parallel vector fields along arbitrary curves in M are orthogonal to their fibers, i.e., for any curve $\gamma$ and parallel vector field v along it, $\dot{\gamma }\perp T_\gamma M$.

Let $\eta :=(p,u):I\rightarrow TM$ be a curve. $\tau$ being a Riemannian submersion, $T_{\eta }TM$ enjoys an orthogonal decomposition in vertical (viz. kernel of $d_{\eta }\tau$) and horizontal subspaces, both $\mu$-dimensional. Identifying each of them with $T_{p}M$, the Sasaki metric at $\eta$ is induced by the quadratic form $\Vert v\Vert ^2+\Vert w\Vert ^2$, where $v=p^\prime$ and $w=u^\prime$. Denoting the covariant derivative and curvature tensor of g by $\nabla$ and R, Sasaki geodesics are given by

$$\begin{aligned} \nabla _v v&= -R(u,w,v),\\ \nabla _v w&= 0. \end{aligned}$$

Algorithms for the computation of the exponential and log map as well as mean geodesic with respect to Sasaki metric, and also an application to corpus callosum longitudinal data as trajectories in Kendall’s shape space are given in [11]. In this case $m=2$ (planar shapes), the shape space can be identified with the complex projective space and the Riemannian curvature tensor is explicitly given in terms of the canonical complex structure and the curvature tensor of the pre-shape space. For $m\ge 3$, computation of R is more delicate.

Next, we present the proposed approach to employ a metric on the tangent bundle. Fix $s\in I$ and let $\gamma _s:J\rightarrow M$ be the geodesic emerging from p(s) with initial velocity $u(s)=\dot{\gamma _s}(0)$. Note that $\gamma _s$ is horizontal. Now, let $\xi _s$ be a vector field along $\gamma _s$. Then by $\int _Jg(\xi _s (t),\xi _s(t))\,dt$ a quadratic form at (p(s), u(s)) is given, which in turn defines a metric, again denoted by $\delta$ in the space of geodesics. Let $H=H(s,t)$ be a path in TM with $\alpha :=H(0,.)$ and $\beta :=H(1,.)$ arc length parametrized geodesics. The energy of H induced by the above quadratic form reads $E(H)=\frac{1}{2}\int _0^1\int _0^1g(H^\prime (s,t),H^\prime (s,t))\,dt\,ds$ ($\xi _s=H^\prime (s,.)$). Let $H^*$ denote the minimizer of E restricted to paths through geodesics, i.e. H(s, .) geodesic for all $s\in I$. Next, we construct time-discrete paths $H_n$ to approximate $H^*$. E achieves its minimum over all paths connecting $\alpha$ to $\beta$ in M if H(., t) is a geodesic for all t (for a proof, we refer the reader to [16, Theorem 3.2]). Now, suppose that $\alpha$ and $\beta$ are close enough, and let $\varPhi$ denote the arc length parametrized geodesic from x to y and $(x_n^*,y_n^*)$ the minimizer of

$$\begin{aligned} E(x,y):=\sum _{i=0}^{n-1}\int _0^1 d^2(\varPhi (x_i,y_i,t),\varPhi (x_{i+1},y_{i+1},t))\,dt \end{aligned}$$

(2)

over $M^n\times M^n$, where $x_0=\alpha (0),y_0=\alpha (1)$, $x_{n}=\beta (0),y_{n}=\beta (1)$. A natural choice for the initial values $x^0$ and $y^0$ is given by the equidistant partition $x^0_i=\varPhi (x_0,x_{n},t_i)$, $y^0_i=\varPhi (y_0,y_{n},t_i)$ with $t_i=\frac{i}{n}$. Then, the desired discrete shortest path reads $H_n=\left( \varPhi (x^*_0,y^*_0,.),\cdots ,\varPhi (x^*_n,y^*_n,.)\right)$.

Notion and computation of the corresponding tangent bundle’s log map is immediate and its linearization ($n=1$) reads $\mathrm {Log}_\alpha \beta$. Similarly, the induced mean of n dense enough geodesics $\gamma _i:I\rightarrow M$ is the geodesic with initial- and endpoints x and y minimizing

$$\begin{aligned} E(x,y)=\sum _i^n\int _0^1d^2(\varPhi (x,y,t),\varPhi (x_i,y_i,t))\,dt \end{aligned}$$

(3)

over $M^2$, where $x_i=\gamma (0)$ and $y_i=\gamma (1)$. A natural choice for the initial value is given by the point-wise mean, viz. $(x ^0,y^0):=$ mean of $(x_i,y_i)$, $i=1,\cdots ,n$. In the sequel, we call $\delta$ the functional-based $L^2$-metric.

Figure 1 shows a geodesic in $TS^2$ as shortest path connecting two $S^2$-geodesics for the Sasaki- and the proposed functional-based metric. For none of them footpoint curves constitute geodesics. However, the functional-based one is closer to the more intuitive shortest path given by the simple point-wise construction $H(s,t)=\varPhi (\alpha (t),\beta (t),s)$.

Fig. 1.

Minimal geodesic in the tangent bundle identified as shortest path connecting two geodesics (red) with respect to Sasaki (left) and functional-based $L^2$-metric (right). (Color figure online)

Fig. 2.

Red geodesics are generated by randomly perturbing the endpoints of the black one. Their mean geodesic with respect to Sasaki (left) and functional-based $L^2$-metric (right) are blue. (Color figure online)

Note that computations of the log map and mean with respect to $\delta$ neither involve the curvature tensor nor decomposition in horizontal and vertical components.

Figure 2 shows the result of an experiment with synthetic spherical data. Geodesics were generated by randomly sampling endpoints following a normal distribution. Computed Sasaki- and functional-mean geodesics consistently provide adequate approximations of the true mean geodesic.

3 Application

In this section, we present the input data, our approach for the estimation of group-wise trends based on regression model, the Hotelling $T^2$ test for group differences, and numerical results (Fig. 3).

3.1 Data Description

We apply the derived scheme to the analysis of group differences in longitudinal femur shapes of subjects with incident and developing osteoarthritis (OA) versus normal controls. The dataset is derived from the Osteoarthritis Initiative (OAI), which is a longitudinal study of knee osteoarthritis maintaining (among others) clinical evaluation data and radiological images from 4,796 men and women of age 45–79. The data are available for public access at http://www.oai.ucsf.edu/. From the OAI database, we determined three groups of shapes trajectories: HH (healthy, i.e. no OA), HD (healthy to diseased, i.e. onset and progression to severe OA), and DD (diseased, i.e. OA at baseline) according to the Kellgren–Lawrence score [6] of grade 0 for all visits, an increase of at least 3 grades over the course of the study, and grade 3 or 4 for all visits, respectively. We extracted surfaces of the distal femora from the respective 3D weDESS MR images (0.37$\times$0.37 mm matrix, 0.7 mm slice thickness) using a state-of-the-art automatic segmentation approach [1]. For each group, we collected 22 trajectories (all available data for group DD minus a record that exhibited inconsistencies, and the same number for groups HD and HH, randomly selected), each of which comprises shapes of all acquired MR images, i.e. at baseline, the 1-, 2-, 3-, 4- and 6-year visits. In a supervised post-process, the quality of segmentations as well as the correspondence of the resulting meshes (8,988 vertices) were ensured.

Fig. 3.

Longitudinal femoral data are divided in groups based on Kellgren–Lawrence score (right). Temporal observations are mapped to trajectories in the shape space (left). An important task is to estimate overall trends within groups.

Fig. 4.

Principal components for Group-wise trends estimated as means of fitted geodesics.

3.2 Numerical Results

We applied the geodesic regression approach to the femoral trajectories described above and represented in Kendall’s shape space. For details and computational aspects, we refer to [12]. Note that geodesic representation provides a less cluttered visualization of the trajectory population making it easier to identify trends within as well as across groups. For the statistical testing of group differences, we employ the manifold-valued Hotelling $T^2$ test described in [11] and present the formulas used therein for the convenience of the reader. Let $x={x_1,\cdots ,x_{n_1}}$ and $y={y_1,\cdots ,y_{n_2}}$ two samples with corresponding Fréchet means $\bar{x}$ and $\bar{y}$, $v_x=\log _{\bar{x}}\bar{y}$, $v_y=\log _{\bar{y}}\bar{x}$. Then the individual group covariances are given by

$$\begin{aligned} W_x&=\frac{1}{n_1}\sum _{i=1}^{n_1}(\log _{\bar{x}}x_i)(\log _{\bar{x}}x_i)^t \\ W_y&=\frac{1}{n_2} \sum _{i=1}^{n_2}(\log _{\bar{y}}y_i)(\log _{\bar{y}}y_i)^t \end{aligned}$$

and the sample $T^2$ statistic reads

$$ t^2=\frac{1}{2}( v_x^tW^{-1}_x v_x + v_y^tW^{-1}_y v_y ). $$

For the estimations of the log map and mean, we employed (2) and (3). We found $t^2$-values 0.0012, 0.000703 and 0.000591 for HH vs. DD, DD vs. HD and HH vs. HD with corresponding p-values 0, 0.011 and 0.033. For the computation of the statistical significance, i.e. p-values, we randomly permuted group memberships of the subject-specific geodesic trends, each identified with its (initial-, endpoint), 1, 000 times. The results reveal clear differences between the group-wise average geodesics demonstrating the descriptiveness of the proposed approach. In particular, the results confirm the obvious differences in group-average trends depicted in the low-dimensional visualization in Fig. 4.

4 Conclusion

We presented a modification of the geodesic hierarchical model introduced in [3] and [11] by employing a discrete geodesic for the tangent bundle of the shape space instead of Sasaki geodesics. Our approach does not involve the Riemannian curvature tensor and provides an efficient approximation. Moreover, we estimated average geodesics and group trends for the example application of femoral longitudinal data incorporating Kendall’s shape space. Furthermore, we employed a manifold-valued Hotelling $t^2$ test confirming that the model well distinguishes the group differences. There are several potential direction for future work. First, it would be interesting to perform the computations utilizing the Sasaki metric and to compare with our numerical result. Second, we would like to extend our approach to independently test for systematic differences in intercept and slope of the trends. Finally, an extension of the approach to higher-dimensional parameters would allow to take further effects into account providing more insight on more complex phenomena.