In medicine, longitudinal studies are of great importance, as they provide information on developmental phenomena that may allow improved prognoses and thus more targeted therapies. To address the problem of strong correlation between measurements taken from the same individual at different times, multivariate hierarchical models were developed for data from Euclidean spaces. In recent years, the enormous potential of longitudinal studies has also come into focus in the analysis of shape and appearance data [GerigFishbaughSadeghi2016]
. In order to obtain as much information as possible from such data, it is necessary to leave the realm of Euclidean vector spaces and turn to methods from Riemannian geometry, as curved manifolds are their natural domains. That is, standard methods from multivariate statistics for analyzing longitudinal data must be transferred to these more general spaces. For hierarchical models this has been done only partially so far. Current parametric hierarchical models for manifold-valued data are based either on geodesics (i.e., generalized straight lines)[MuralidharanFletcher2012, Singh_ea2013, Singh_ea2016, BoneColliotDurrleman2018] or general trajectories [ChakrabortyBanerjeeVemuri2017, SrivastavaKlassen2016]; further, an approach based on nonparametric curves was presented in [CampbellFletcher2018]. As an alternative to the geodesic models, the authors of [NavayazdaniHegevonTycowicz2019]
proposed to use a different Riemannian structure that allows faster algorithms for high-dimensional data. However, to our knowledge, there is noin-between
hierarchical model for parametric trends of higher-order (i.e., based on generalized polynomials) that (a) allows to model more complicated trends and (b) is efficient enough to handle large data bases due to a small number of degrees of freedom. Since many phenomena, e.g., cyclic motion of cardiac anatomy, are only inadequately characterized by geodesic models, such a model could be useful in numerous applications.
Therefore, in this paper, we propose a higher order hierarchical model for the analysis of longitudinal manifold-valued data. By modelling, in a first step, subject-wise trends as splines consisting of generalized Bézier curves [PopielNoakes2007] (i.e., generalized polynomials), we are able to capture a vast range of phenomena—not least periodic ones. To obtain the trajectories from the data, we rely on regression with Bézier splines developed in [Hanik_ea2020]. In the second step, the obtained trajectories are considered as perturbations of a common mean (curve). To this end, we adapt the functional-based metric from [NavayazdaniHegevonTycowicz2019] to compare the obtained subject trajectories within the space of Bézier splines. Thereby, we are able to analyze the inter-subject variability without interference from correlated measurements. For efficient calculations in the obtained space of Bézier splines, we rely on the geodesic calculus introduced in [Heeren_ea2018, RumpfWirth2014] and also used in [NavayazdaniHegevonTycowicz2019]
. An implementation of the presented approach is publicly available as part of the open-source libraryMorphomatics [Morphomatics].
We validate our model using data from the Osteoarthritis Initiative. To the best of our knowledge, this is the first time that whole trajectories (representing progression of osteoarthritis) and not only momentary states are classified.
2 Hierarchical model
2.1 Bézier curves in Riemannian geometry
For preparation we first recall some facts from Riemannian geometry and about Bézier curves on manifolds. In the following “smooth” means “infinitely often differentiable”.
A Riemannian manifold is a differentiable manifold together with a Riemannian metric that assigns to each tangent space a smoothly varying scalar product; the metric also induces a distance function . For every Riemannian manifold there is a unique Levi-Civita connection . Given two vector fields on , it allows us to differentiate along ; the result is again a vector field, which we denote by . The connection allows us to define geodesics (i.e., generalized straight lines). A curve is called geodesic if its acceleration vanishes identically, i.e., , where . It is a useful fact that each element of has a so-called normal convex neighbourhood . Any two points can be joined by a unique length-minimizing geodesic that does not leave . Throughout this paper, we always assume to work in a convex neighbourhood in order to use this property.
Next, we recall manifold-valued Bézier curves and splines. For clarity, we restrict the domain of definition of Bézier curves to in this work; in general, the curves can always be reparametrized.
A set of control points defines a continuously differentiable Bézier curve of order according to the generalized de Casteljau algorithm
by ; see Ref. [PopielNoakes2007]. Several such curves, say , can be joined to a continuously differentiable spline [GousenbourgerMassartAbsil2019]: For let be the control points of the curves with
for all but and . Then, the corresponding Bézier spline is defined by
If and the first and last segment of are at least cubic, can be closed. Then, is and closed if and only if (1) extends cyclically as discussed in [Hanik_ea2020].
In the following, we set if is not closed ( if is closed) and denote the set of distinct control points of by . In the non-closed case this means
while is left out for closed .
Regression with intrinsic Bézier splines [Hanik_ea2020] models the relationship between an independent scalar variable and an -valued dependent variable as a Bézier spline (with a fixed number of control points). That is, for data pairs , the minimizer (represented by its control points) of the sum-of-squared energy
models the relationship between and .
2.2 The Model
In this section we introduce the nonlinear hierarchical model. For this, we define the set of Bézier splines in a normal convex neighbourhood with a given number of segments and degrees:
In the following we assume and to be fixed and, hence, omit the indices for readability.
Consider that for subjects measurements of an independent scalar variable and a manifold-valued dependent variable are given, that is,
(Note that the number of measurements can be different for each subject.) Such data can, for example, arise in a longitudinal study that observes shape developments in several individuals.
In a first step, we model the individual trends by regression with Bézier splines of fixed type; that is, for each we perform spline regression with respect to the data and, thus, obtain Bézier splines , that represent the intra-subject trends.
In the second step, we model the individual trends as perturbations of a common mean trajectory. We do this by considering as a submanifold of the manifold of all smooth curves in . Through the identification of the curves and their control points, tangent vectors at the control points naturally translate into generalized Jacobi fields along the curves [BergmanGousenbourger2018]. The metric from Ref. [SrivastavaKlassen2016, Sec. 3.3] can thus be restricted to to yield a natural Riemannian metric. In particular, it induces a natural distance between two Bézier splines that can be efficiently evaluated using variational time-discretization [RumpfWirth2014], as described in the following.
Let . A path between and through may be represented as a parametrized surface in because it induces a map with and . A geodesic between and is then defined as the minimizer of the path energy Discretizing in and identifying splines with their control points, we obtain a discrete -geodesic between and as the minimizer of the discrete path energy The integral can be evaluated using a suitable quadrature rule. In order to approximate the minimizer of the discrete energy, we extend the iterative procedure from Ref. [NavayazdaniHegevonTycowicz2019] to our setting: We compute the discrete -geodesics between two curves by iteratively performing spline regression. First, we initialize the control points of the inner curves equidistantly along the geodesics that connect the corresponding control points of and . Then, the inner curves are updated so that they lie “in the middle” of their neighbors; to this end, we replace them with the result of spline regression with respect to data points that are given by (equidistant) evaluations of the neighboring curves. The procedure is summarized in Alg. 1.
Next, we discuss the computation of the discrete -mean of curves , with which we approximate the common mean (curve). It is the spline minimizing , s.t. , , where denotes the control points of the discrete geodesic between and . It can be computed with an alternating optimization scheme. As an initialization of the control points of , we choose the means of the corresponding control points of the data curves. Then, in alternating fashion, discrete geodesics towards the mean are computed and, subsequently, the mean is updated by spline regression, since it has to lie “in the middle” of the innermost elements of the discrete geodesics. The procedure is summarized in Alg. 2.
3 Experimental evaluation
Hierarchical models provide a principled way of analyzing longitudinal data. To demonstrate this, we perform group-wise analysis of femoral shape trajectories that have been collected from the Osteoarthritis Initiative (OAI). Furthermore, in order to demonstrate that the presented method is not limited to the estimation of average, group-level trends, we derive a statistical descriptor for shape trajectories in terms of the principal component scores (i.e., the coefficients encoding the trajectories within the basis of principal modes) and use it fortrajectory classification.
The OAI is a longitudinal study of knee OA comprising (among others) clinical evaluation data and radiological images from 4,796 men and women of age 45-79 publicly available at https://nda.nih.gov/oai/. 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 [KellgrenLawrence1957] 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. Using an automatic segmentation approach [Ambellan_ea2019], we extracted surfaces of the distal femora from the respective 3D weDESS MR images (0.370.37 mm matrix, 0.7 mm slice thickness). For each group, we assembled 22 trajectories (all available data for group DD except one subject 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-, 6-, and 8-year visits. As notion of shape space we employ the differential coordinates model [vonTycowicz_ea2018] that allows for closed-form evaluation of Riemannian operations and therefore facilitates fast and numerically robust processing.
As a first application, we estimate a hierarchical model for the HD group. We employ cubic Bézier curves to model the individual trends. This choice is motivated by the findings in [Hanik_ea2020], where cubic models were found to adequately capture the inherent nonlinear shape developments due to OA in a cross-sectional regression-based analysis. Time discrete computations are performed based on 2-geodesics—employing finer discretizations have been found to provide no further improvements for the dataset under study. The estimated group-level trend is visualized in Fig. 1. The determined shape changes consistently expose OA related malformations of the femur, most prominently changes along the ridge of the cartilage plate that are characteristic regions for osteophytic growth. Note, that only minute bone remodeling can be observed for the first half of the captured interval, whereas bone malformations develop more rapidly after four years time. This behavior suggests that there are nontrivial higher order phenomena involved, which geodesic models cannot adequately describe.
For the classification we train a simple support vector machine (linear kernel) on 65-dimensional descriptors (coefficients w.r.t. the PGA modes from an approximated Gram matrix[Heeren_ea2018]) in a leave-one-out cross-validation setup.
The percentage of correctly classified trajectories is 64%. The corresponding confusion matrix is given as inset. Performing the same experiment with a Euclidean model[Cootes_ea1995] results in 59% correct classifications demonstrating the advantage of our Riemannian model.
4 Conclusion and Future Work
We presented a hierarchical statistical model that is based on intrinsic, higher-order Bézier splines and thus allows analyzing a wide range of phenomena with non-monotonous shape changes. To the best of our knowledge this is the first Riemannian model that is neither bound to constraints of geodesicity nor based on non-parametric designs.
A promising direction for future work is to extend the hierarchical model presented to account for subject-specific shifts in the stage of evolution. Adapting the proposed distance to partial trajectories could improve the estimation of group-level trends from longitudinal observations with highly varying inter-individual age range or disease stage coverage.
Furthermore, hierarchical models provide a principled way of analyzing longitudinal data. In this regard, the presented method is not limited to estimating the average trends at, group-level, but can also be employed to assess the variance and principal modes of the distribution of the studied trajectories. This opens up a multitude of applications such as hypothesis testing and Bayesian reconstruction, which we will address in the future.
5 Compliance with Ethical Standards
This research study was conducted retrospectively using human subject data made available in open access. Ethical approval was *not* required as confirmed by the license attached with the open access data.
We are grateful for the funding by DFG111Deutsche Forschungsgemeinschaft (DFG) through Germany’s Excellence Strategy – The Berlin Mathematics Research Center MATH+ (EXC-2046/1, project ID: 390685689) and BMBF222Bundesministerium für Bildung und Forschung (BMBF) through BIFOLD - The Berlin Institute for the Foundations of Learning and Data (ref. 01IS18025A and ref 01IS18037A) as well as for the provision of the data set by the OAI333OAI is a public-private partnership comprised of five contracts (N01-AR-2-2258; N01-AR-2-2259; N01-AR-2-2260; N01-AR-2-2261; N01-AR-2-2262) funded by the National Institutes of Health, a branch of the Department of Health and Human Services, and conducted by the OAI Study Investigators. Private funding partners include Merck Research Laboratories; Novartis Pharmaceuticals Corporation, GlaxoSmithKline; and Pfizer, Inc. Private sector funding for the OAI is managed by the Foundation for the National Institutes of Health. This manuscript was prepared using an OAI public use data set and does not necessarily reflect the opinions or views of the OAI investigators, the NIH, or the private funding partners..