A Nonlinear Hierarchical Model for Longitudinal Data on Manifolds

by   Martin Hanik, et al.

Large longitudinal studies provide lots of valuable information, especially in medical applications. A problem which must be taken care of in order to utilize their full potential is that of correlation between intra-subject measurements taken at different times. For data in Euclidean space this can be done with hierarchical models, that is, models that consider intra-subject and between-subject variability in two different stages. Nevertheless, data from medical studies often takes values in nonlinear manifolds. Here, as a first step, geodesic hierarchical models have been developed that generalize the linear ansatz by assuming that time-induced intra-subject variations occur along a generalized straight line in the manifold. However, this is often not the case (e.g., periodic motion or processes with saturation). We propose a hierarchical model for manifold-valued data that extends this to include trends along higher-order curves, namely Bézier splines in the manifold. To this end, we present a principled way of comparing shape trends in terms of a functional-based Riemannian metric. Remarkably, this metric allows efficient, yet simple computations by virtue of a variational time discretization requiring only the solution of regression problems. We validate our model on longitudinal data from the osteoarthritis initiative, including classification of disease progression.



page 4


A Hierarchical Geodesic Model for Longitudinal Analysis on Manifolds

In many applications, geodesic hierarchical models are adequate for the ...

Modeling Longitudinal Data on Riemannian Manifolds

When considering functional principal component analysis for sparsely ob...

Prediction of the progression of subcortical brain structures in Alzheimer's disease from baseline

We propose a method to predict the subject-specific longitudinal progres...

Longitudinal data analysis using matrix completion

In clinical practice and biomedical research, measurements are often col...

Learning distributions of shape trajectories from longitudinal datasets: a hierarchical model on a manifold of diffeomorphisms

We propose a method to learn a distribution of shape trajectories from l...

Rigid Motion Invariant Statistical Shape Modeling based on Discrete Fundamental Forms

We present a novel approach for nonlinear statistical shape modeling tha...

Modeling Firn Density through Spatially Varying Smoothed Arrhenius Regression

Scientists use firn (compacted snow) density models as a function of dep...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

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 no


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 library

Morphomatics [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.

2.3 Computation

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.

Input: with control points , respectively
    Output: Control points of the discrete -geodesic from to

  for  do
  end for
     for  do
        for  do
        end for
     end for
  until convergence
Algorithm 1 Discrete -geodesic in . Solving spline regression by minimizing (2) is denoted by reg.

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.

Input: with control points , discretization parameter
    Output: Control points of the mean curve

     for  do
     end for
     for  do
        for  do
        end for
     end for
  until convergence
Algorithm 2 Mean trajectory. Computation of Fréchet mean and -geodesic are denoted mean and -geo, respectively.

3 Experimental evaluation

Figure 1: Mean of cubic femoral trends of 22 subjects evaluated at 5 equidistant points. The surface distance to the baseline (value of the computed mean at ) is color coded wherever the distance is larger than 0.25 mm.

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 for

trajectory 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.

HH 19 1 2
DD 2 11 9
HD 4 6 12

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.

6 Acknowledgments

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..