**Multivariate setting** We extended the previous experiment to a two-dimensional setting. This means that the function $g$ is more complex as both its input and output are two-dimensional now. The results are shown in Figure 4.15. On the individual level the trajectories learned are quite good. However the learned function is not matching the ground truth. This is due to the breakpoints being concentrated in the space, which does not allow radial basis functions to cover enough space.

To remedy this problem, we changed the centers of the radial basis functions from the data points themselves to regularly spaced points on a grid. This is a relaxation of the representer theorem, but allows for more flexibility in the learning of the function. This also avoids ill-conditioning that sometimes happen when centers of the radial-basis functions are too close to each other. With a regular grid, the learning is much better, and we can even approximate complicated functions, as shown in Figure 4.16.

### 4.3.3 Discussion

These experiments show that in a simple linear case we can perform a reliable estimation of complex treatment function for the change of slope. The extension to the case where the function also takes the treatment dose and over covariates as input is straightforward. However, if we increase the number of input variables in the learned treatment function, or if the observations dimension increases, the RKHS model requires an exponential amount of points to achieve proper learning. This is very limiting as in practice the amount of data used to effectively learn the treatment function is the total number of breakpoints over the population. As our targeted application concerns clinical trials, this number will be limited to the number of patients in the treated arm, which is far less than the ten thousand points used in our simulated experiments.

This burden can be alleviated by choosing other forms for the treatment function, which is the next step towards a practical application of the Piece-wise Geodesic model. The most simple choice would be a linear function. Any parametric family of functions with a low number of parameters can be used as long as we can compute the gradient of the function with respect to its parameters. This is necessary to perform the gradient descent, as described in the estimation algorithm for the Piece-wise Geodesic model in the next chapter. On this occasion we will also discuss the other limitation of the model: the difficulty to learn the treatment function in a non-Euclidean manifold, namely the logistic curves manifold.
<figure><img src="image_52.png" /><figcaption>Figure 4.15: Piece-wise Geodesic model results on simulated data in dimension 2. The results are shown after 1000 iterations of MCMC-SAEM, with an optimization of the weights $\omega_{i,j}$ every 10 iterations. The first row corresponds to three individuals. Both biomarkers are shown on the same plot, the true trajectory is shown in red with the observations in blue and orange, while the learned model is shown in green. On the bottom plot we represent the function $g$ learned by the Piece-wise Geodesic model against the ground truth with the repartition of data points, which are the centers of the radial basis functions composing the learned $g$ function.</figcaption></figure>
<figure><img src="image_53.png" /><figcaption>Figure 4.16: Piece-wise Geodesic model results on simulated data in dimension 2. We plotted the function *g* learned by the Piece-wise Geodesic model against the ground truth with the repartition of data points (in blue) and the centers of the RBF (in orange).</figcaption></figure>
## 4.4 Conclusion

Treatment effect modelling introduces a new layer of complexity atop the natural disease history modelling task. Being able to approximate the state of the patient, had he not followed a treatment, is permitted by the disease progression models. Based on this latent progression, which can be seen as a digital twin of the patient, we showed that we were able to reconstruct the treatment effect under some conditions on the disease progression model. We believe that this technique has the potential to one day replace placebo arms in some clinical trials, thereby resolving the ethical dilemma of not giving a potentially life-saving drug to half the participants in the trial. We do not claim to have proven the reliability of the method yet, there is still plenty of work to be done. But we hope that we have laid the foundations for future research in this direction.

Current research in neurodegenerative diseases is hoping to find cures or disease-modifying treatments to slow the progression. Treatments with long-term modification of a patient's trajectory call for other types of model than the regressions used before. Estimating if a treatment reduced the rate of progression is easily assessed in a linear setting, less so in the non-linear case. We proposed a new model based on the DCM, called Piece-wise Geodesic, where we introduced a treatment function accounting for the changes of trajectory after every new prescription. We showed that this function could be learned with synthetic data. We have not yet been able to apply this model to real-life data because there is no such treatment for the most common neurodegenerative diseases, not mentioning the fact that the study needs to provide sufficient data for the model to be estimated. However β-amyloïd targeting therapies in Alzheimer's disease have made significant progress in the last two years, with the first drug having recently been approved by the US Food and Drug Administration (FDA)<fnref n="37" />. We can thus expect data from cohorts following patients over the next few years. The Piece-wise Geodesic model is thus an anticipation of the needs of the field, and we hope that it will be applicable in a near future. In the meantime, we provide some ideas for additional developments: simplifying the parametrization of the treatment function, allowing the model to have discontinuities at the breakpoints and introduce an individual latent variable measuring how well a patient responds to treatment.

<footnote n="37">³⁷O. o. t. Commissioner, *FDA Converts Novel Alzheimer's Disease Treatment to Traditional Approval*, en, Jul. 2023.</footnote>
<figure><img src="image_54.png" /></figure>
# Chapter 5

## Estimation

This chapter is structured differently from the other chapters of this thesis. Indeed this is not the presentation of a complete work on a particular subject with some applications to clinical data. Instead it is a patchwork of small parts from all the models presented in the previous chapters. It is however unified by the theme: the estimation algorithm for the DCM and all its proposed variants. The goal was to regroup all the questions about optimization and parameters estimation in the same chapter, so that everything is simpler to understand.

The main algorithm is a variant of the Expectation-Maximization algorithm. We will first present this algorithm, which was chosen for the estimation of the DCM model in the original work<fnref n="1" />. Then we will present the modifications that we brought to it to adapt to each particular model. In order we will describe variants for the Geodesic Bending model, for the mixture models, and finally for the piece-wise geodesic model.

### Contents of the chapter

<table><tr><td>5.1</td><td>The MCMC-SAEM algorithm</td><td>140</td></tr><tr><td>5.2</td><td>Variants of the MCMC SAEM</td><td>143</td></tr><tr><td>5.2.1</td><td>MMCMC SAEM for mixture models</td><td>143</td></tr><tr><td>5.2.2</td><td>Alternating maximization algorithm</td><td>146</td></tr><tr><td>5.2.3</td><td>Gradient maximization step</td><td>148</td></tr><tr><td>5.2.4</td><td>Discussion</td><td>150</td></tr></table>

<footnote n="1">¹J.-B. Schiratti, S. Allassonnière, A. Routier, et al., “A Mixed-Effects Model with Time Reparametrization for Longitudinal Univariate Manifold-Valued Data,” en, in *Information Processing in Medical Imaging*, S. Ourselin, D. C. Alexander, C.-F. Westin, et al., Eds., ser. Lecture Notes in Computer Science, Cham: Springer International Publishing, 2015, pp. 564–575.</footnote>
## 5.1 The MCMC-SAEM algorithm

The main goal of the estimation algorithm is to find the model parameters $\theta$ that maximize the likelihood $q$:

$$ q(\mathbf{y}, \theta | \theta_{\text{hyper}}) = \int q(y, \mathbf{z}, \theta | \theta_{\text{hyper}}) d\mathbf{z} \qquad (5.1.1) $$

$$ = \int q(y | \mathbf{z}, \theta, \theta_{\text{hyper}}) q(\mathbf{z}, \theta | \theta_{\text{hyper}}) d\mathbf{z} \qquad (5.1.2) $$

where $\mathbf{z} = (\mathbf{z}_{pop}, (\mathbf{z}_i)_i)$ are the latent parameters, composed of the population and individual parameters. However this integral is not always tractable since the latent variables are unknown. However the Expectation-Maximization (EM) algorithm, as introduced by<fnref n="2" />, was designed for such cases. This iterative algorithm alternates between an estimation of the latent variables, thus allowing to compute the expected likelihood, and a maximization of the model parameters based on this expectation. The EM algorithm thereby produces a maximum likelihood estimate or, in the case of a Bayesian model, a maximum a posteriori. The algorithm is described in Algorithm 1.

<caption>**Algorithm 1: Expectation-Maximization**</caption>

1. $\theta^{(0)} \leftarrow \theta_0$ $k \leftarrow 0$ **while** $k \le n$ **do**
2. &nbsp;&nbsp;&nbsp;&nbsp;$k \leftarrow k + 1$
3. &nbsp;&nbsp;&nbsp;&nbsp;Expectation step: Compute $Q(\theta, \theta^{(k)}) = \mathbb{E}_{p(z|y;\theta^{(k)})} (\log p(y, z; \theta))$
4. &nbsp;&nbsp;&nbsp;&nbsp;Maximization step: Update $\theta^{(k+1)} = \underset{\theta}{\operatorname{argmax}} Q(\theta, \theta^{(k)})$

5. **end**

**Result:** return $\theta^{(k)}$

At iteration $k$, the E-step computes the expectation of the likelihood over the latent variables with the model parameters $\theta^{(k)}$ being fixed. Then the expectation being a function of the model parameters $\theta$, the M-step finds the best parameters.

In linear cases the EM is sufficient, as in<fnref n="3" />. However in most non-linear cases, the expectation is intractable due to the integration over the latent variables. This is where the Stochastic Approximation EM shines<fnref n="4" />. The expectation in the E-step is not computed explicitly, but approximated by drawing samples of latent variables $z^{(k)} \sim p(z|y; \theta^{(k)})$. The algorithm only requires one sample for $\mathbf{z}$ per step, which avoids computing a full Monte-Carlo approximation of the integral in the expectation formula. The convergence towards the true integral is ensured by a Robbins-Monro scheme<fnref n="5" />, where after drawing the sample $z^{(k)}$ and computing $\log p(y, z^{(k)}; \theta)$, one computes the approximated expectation $Q_k(\theta) = Q_{k-1}(\theta) + \epsilon_k (\log p(y, z^{(k)}; \theta) - Q_{k-1})$, with $\sum_k \epsilon_k = \infty$ and $\sum_k \epsilon_k^2 < \infty$.

<footnote n="2">²A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum Likelihood from Incomplete Data Via the EM Algorithm,” en, *Journal of the Royal Statistical Society: Series B (Methodological)*, vol. 39, no. 1, pp. 1–22, 1977.</footnote>

<footnote n="3">³G. Verbeke and E. Lesaffre, “A Linear Mixed-Effects Model With Heterogeneity in the Random-Effects Population,” *Journal of the American Statistical Association*, vol. 91, no. 433, pp. 217–221, 1996.</footnote>

<footnote n="4">⁴B. Delyon, M. Lavielle, and E. Moulines, “Convergence of a Stochastic Approximation Version of the EM Algorithm,” *The Annals of Statistics*, vol. 27, no. 1, pp. 94–128, 1999.</footnote>

<footnote n="5">⁵H. Robbins and S. Monro, “A Stochastic Approximation Method,” *The Annals of Mathematical Statistics*, vol. 22, no. 3, pp. 400–407, Sep. 1951.</footnote>
However an additional difficulty lies in the fact that $p(z | y; \theta^{(k)})$ might be unknown. The Bayes rule gives us that:

$$p(z | y; \theta^{(k)}) = \frac{p(y | z; \theta^{(k)})p(z; \theta^{(k)})}{p(y; \theta^{(k)})} = \frac{p(y | z; \theta^{(k)})p(z; \theta^{(k)})}{\int p(y | z; \theta^{(k)})p(z; \theta^{(k)})dz} \quad (5.1.3)$$

where $p(y | z; \theta^{(k)})$ is given by the model itself, $p(z; \theta^{(k)})$ is a known prior on the latent variables and the denominator is a constant. We thus know the distribution of the latent variables conditioned on the observations up to the normalization constant. The Metropolis-Hastings algorithm<fnref n="6" /> can then be used to sample from this distribution. The Metropolis-Hastings algorithm is based on a Markov chain, hence the name Monte-Carlo Markov Chain SAEM. This algorithm is proven to converge to a local maximum of the observed likelihood $p(y|\theta)$ (see<fnref n="7" />,<fnref n="8" />,<fnref n="9" />). However the convergence has only been proven for the curved exponential family, which is the family of distributions for which the log-likelihood can be written as: $\forall \theta \in \Theta, \log q(\mathbf{y}, \mathbf{z}; \theta) = -\Phi(\theta) + \langle S(\mathbf{y}, \mathbf{z}), g(\theta) \rangle$ where $\Phi$ and $g$ are smooth functions, $S$ are a function of $\mathbf{z}$, $\mathbf{y}$ called the sufficient statistics. To be noted, all the models presented in this work fall into the curved exponential family, except for the ordinal model.

The sufficient statistics are to be understood as a summary of the required information from the latent variables $\mathbf{z}$ and the observations $\mathbf{y}$. Once the summary statistics are computed, one can derive the log-likelihood linearly from it by definition. Therefore the Robbins-Monro convergence scheme for the computation of the expectation is equivalent to applying the Robbins-Monro scheme directly on the sufficient statistics, as described in the algorithm below. Then the maximization step consists in finding the argmax of the expected log-likelihood based on the sufficient statistics. The maximization in the base DCM model formulation is direct, so the new parameters $\theta^{(k)}$ are computed explicitly in closed form.

The pseudo-code is given in the complete Algorithm 2.

$\pi$ is a proposition law that is known and from which we can easily sample. Typically we use a Gaussian centered on $z^{(k-1)}$ with an adaptive variance<fnref n="10" />, i.e. we automatically scale the variance so that the acceptance rate stays around 30%<fnref n="11" />. In practice we use a Metropolis-Hastings within Gibbs sampler, which corresponds to a sequential sampling of the coordinates of $\mathbf{z}$.

The algorithm requires a burn-in phase, which corresponds to a certain number of iterations $n_{burn-in}$ during which $\forall k \le n_{burn-in}, \epsilon_k = 1$. This corresponds to a memory-less state, allowing first the Markov chain to converge. Once the Markov chain is stable, the sampling step corresponds to drawing the latent variables from their posterior distribution. The Robbins-Monro scheme in the burn-out phase only allows for the convergence of the model parameters to a maximum a posteriori.

We refer the reader to the appendix 15 for complete formulas of the log-likelihood. We did not

<footnote n="6">⁶W. K. Hastings, “Monte Carlo sampling methods using Markov chains and their applications,” *Biometrika*, vol. 57, no. 1, pp. 97–109, Apr. 1970.</footnote>

<footnote n="7">⁷E. Kuhn and M. Lavielle, “Coupling a stochastic approximation version of EM with an MCMC procedure,” fr, *ESAIM: Probability and Statistics*, vol. 8, pp. 115–131, 2004.</footnote>

<footnote n="8">⁸S. Allassonnière and E. Kuhn, “Stochastic algorithm for Bayesian mixture effect template estimation,” en, *ESAIM: Probability and Statistics*, vol. 14, pp. 382–408, Jan. 2010.</footnote>

<footnote n="9">⁹S. Allassonnière and E. Kuhn, “Convergent stochastic Expectation Maximization algorithm with efficient sampling in high dimension. Application to deformable template model estimation,” en, *Computational Statistics & Data Analysis*, vol. 91, pp. 4–19, Nov. 2015.</footnote>

<footnote n="10">¹⁰H. Haario, E. Saksman, and J. Tamminen, “An Adaptive Metropolis Algorithm,” *Bernoulli*, vol. 7, no. 2, pp. 223–242, 2001.</footnote>

<footnote n="11">¹¹S. Chib and E. Greenberg, “Understanding the Metropolis-Hastings Algorithm,” *The American Statistician*, vol. 49, no. 4, pp. 327–335, Nov. 1995.</footnote>
<table><caption>**Algorithm 2:** Monte Carlo Markov Chain Stochastic Approximation Expectation-Maximization</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>θ<sup>(0)</sup> ← θ<sub>0</sub> S<sub>0</sub> ← 0</td></tr><tr><td>2</td><td>z<sup>(0)</sup> ← z<sub>0</sub></td></tr><tr><td>3</td><td>(ε<sub>k</sub>)<sub>k≥0</sub> (for the Robbins-Monro scheme)</td></tr><tr><td>4</td><td>k ← 0</td></tr><tr><td>5</td><td>while k ≤ n do</td></tr><tr><td>6</td><td>k ← k + 1</td></tr><tr><td>7</td><td>Simulation step: Given (y, θ<sup>(k)</sup>), sample z<sup>(k)</sup> :</td></tr><tr><td>8</td><td>Draw z* ~ π(z<sup>(k-1)</sup>)</td></tr><tr><td>9</td><td>Compute acceptance ratio: α = max(<sup>q(y,z<sup>*</sup>|θ)</sup>⁄<sub>q(y,z<sup>(k-1)</sup>|θ)</sub>, 1)</td></tr><tr><td>10</td><td>Accept with probability α: z<sup>(k)</sup> ← z*</td></tr><tr><td>11</td><td>or reject: z<sup>(k)</sup> ← z<sup>(k-1)</sup></td></tr><tr><td>12</td><td>Stochastic Approximation step : S<sub>k</sub> ← S<sub>k-1</sub> + ε<sub>k</sub> (S(y, z<sup>(k)</sup>) - S<sub>k-1</sub>)</td></tr><tr><td>13</td><td>Maximization step : θ<sub>k</sub> = argmax<sub>θ</sub> (− log Φ(θ) + &lang;S̃<sub>k</sub>, g(θ)&rang;)</td></tr><tr><td>14</td><td>end</td></tr><tr><td>15</td><td>Result: return θ<sup>(k)</sup></td></tr></tbody></table>

derive the sufficient statistics here as it has already been done previously in<fnref n="12" /> and<fnref n="13" />. We will give one example to understand how it works:

The log-likelihood term for the prior of the time-shifts is $\log q_{prior}((\tau_i)_i | \theta) = -N \log(\sigma_{\tau} \sqrt{2\pi}) - \frac{1}{2\sigma_{\tau}^2} \sum_{i=1}^{N} (\tau_i - \tau_i)^2$. The summary statistics associated to the model parameters $\bar{\tau}$ and $\sigma_{\tau}$ are $S_1 = \sum_{i=1}^{N} \tau_i$ and $S_2 = \sum_{i=1}^{N} \tau_i^2$. From these sufficient statistics we can derive the explicit maximization rules: $\bar{\tau} \leftarrow \frac{1}{N} S_1$ and $\sigma_{\tau} \leftarrow \frac{1}{N} (S_2 - 2S_1 + \bar{\tau})$.

In its implementation, there is no convergence criterion for the algorithm. We currently fix an arbitrary number of iterations which is large enough for the sampling Markov chain to stabilize. However some improvements could be made by using a stopping criterion such as the Gelman and Rubin convergence criterion for MCMC<fnref n="14" />.

We will now delve into the difference brought by the several models.

<footnote n="12">¹²I. Koval, “Learning Multimodal Digital Models of Disease Progression from Longitudinal Data: Methods & Algorithms for the Description, Prediction and Simulation of Alzheimer’s Disease Progression,” en, Ph.D. dissertation, Institut Polytechnique de Paris, Jan. 2020.</footnote>

<footnote n="13">¹³J.-B. Schiratti, “Methods and algorithms to learn spatio-temporal changes from longitudinal manifold-valued observations,” en, Ph.D. dissertation, Université Paris Saclay (COMUE), Jan. 2017.</footnote>

<footnote n="14">¹⁴A. Gelman and D. B. Rubin, “Inference from Iterative Simulation Using Multiple Sequences,” *Statistical Science*, vol. 7, no. 4, pp. 457–472, Nov. 1992.</footnote>
## 5.2 Variants of the MCMC SAEM

### 5.2.1 MMCMC SAEM for mixture models

We saw that the MCMC SAEM was the base estimation algorithm for the DCM model. Moreover the EM algorithm is also a staple for mixture model estimation. Combining the estimation of a mixture model and the non-linear mixed-effect models has already been studied<fnref n="15" />. However the theoretical properties are often not enough as the MCMC SAEM tends, in practice, to suffer from local maxima attraction. In particular, changing cluster assignments is known to be challenging in such methods; it is a problem referred to as trapping states. Adding a tempering scheme has been shown to alleviate this issue by flattening the target distribution and thus easing the exploration of the parameter space<fnref n="16" />. A mixture model on top of a disease progression model was also proposed in<fnref n="17" />, the clustering during estimation was handled with hard labels and a probability for individuals to switch from one cluster to another.

We remind the reader of our formulation of the mixture of DCM models:

$$Q(\mathbf{y} | \mathbf{z}; \theta^1, ..., \theta^L) = \sum_{c=1}^{L} \pi^c q(\mathbf{y} | \mathbf{z}; \theta^c)$$

and the mixture of individual parameters within the DCM model:

$$q(\begin{bmatrix} \tau_i \\ \xi_i \end{bmatrix} | \theta) = \sum_{c=1}^{L} \mathbb{1}_{c_i=c} \phi\left(\begin{bmatrix} \tau_i \\ \xi_i \end{bmatrix}; \begin{bmatrix} \bar{\tau}^c \\ \bar{\xi}^c \end{bmatrix}, \Sigma_{\tau, \xi}^c\right)$$

which have both the same log-likelihood for the observations once we introduce the latent classes of the individuals $c_i$:

$$q(\mathbf{y}_i | \mathbf{z}_i; \theta^1, ..., \theta^L) = \sum_{c=1}^{L} \mathbb{1}(c_i = c) q(\mathbf{y}_i | \mathbf{z}_i; \theta^c) \qquad (5.2.1)$$

In the estimation, we compute the latent variables $\pi_i^c = \mathbb{P}(c_i = c | \mathbf{z}_i, \mathbf{z}_{pop}, \theta)$. Contrarily to the other latent variables, they can be directly computed in the expectation step, as:

$$\pi_{K+1,i}^c = \frac{\pi_K^c q(y_i|\theta_K^c, z_{K,i}^c, z_{K,pop}^c)) p(z_{K,i}^c|\theta_K^c, z_{K,pop}^c))}{\sum_j \pi_K^j q(y_i|\theta_K^j, z_{K,i}^j, z_{K,pop}^j)) p(z_{K,i}^j|\theta_K^j, z_{K,pop}^j))} \qquad (5.2.2)$$

Contrarily to<fnref n="18" /> where the latent variable is directly $c_i$, forcing individuals to attach to a cluster, our use of soft cluster labelling avoids cluster freeze and trapping states. We provide below the algorithm for the mixture of DCM models.

<footnote n="15">¹⁵M. Lavielle and C. Mbogning, “An improved SAEM algorithm for maximum likelihood estimation in mixtures of non linear mixed effects models,” *Statistics and Computing*, vol. 24, no. 5, pp. 693–707, Sep. 2014.</footnote>

<footnote n="16">¹⁶S. Allassonnière and J. Chevallier, “A new class of stochastic EM algorithms. Escaping local maxima and handling intractable sampling,” en, *Computational Statistics & Data Analysis*, vol. 159, p. 107–159, Jul. 2021.</footnote>

<footnote n="17">¹⁷V. Debavelaere, S. Durrleman, S. Allassonnière, et al., “Learning the Clustering of Longitudinal Shape Data Sets into a Mixture of Independent or Branching Trajectories,” en, *International Journal of Computer Vision*, vol. 128, no. 12, pp. 2794–2809, Dec. 2020.</footnote>

<footnote n="18">¹⁸V. Debavelaere, S. Durrleman, S. Allassonnière, et al., “Learning the Clustering of Longitudinal Shape Data Sets into a Mixture of Independent or Branching Trajectories,” en, *International Journal of Computer Vision*, vol. 128, no. 12, pp. 2794–2809, Dec. 2020.</footnote>
<caption>**Algorithm 3: M-MCMC-SAEM estimation**</caption>

1 Initialization: $\pi_0$, $(\theta_0^c)$, $(\mathbf{z}_{0,i}^c)$, $(\pi_{0,i}^c)$
2 for $K = 0...N$ do
3 &nbsp;&nbsp;&nbsp;&nbsp;for $c = 1...L$ do
4 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;• **E step**
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute probability of individual $i$ to belong to cluster $c$: $\pi_{K+1,i}^c$
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**Population parameters estimation**
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Sample $z_{*,pop}^c$
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute acceptance ratio $\alpha = 1 \wedge \frac{q(\mathbf{y}|\theta_{K}^c, (z_{K,i}^c)_i, z_{*,pop}^c)q((z_{K,i}^c)_i|\theta_{K}^c, z_{*,pop}^c)}{q(\mathbf{y}|\theta_{K}^c, (z_{K,i}^c)_i, z_{K,pop}^c)q((z_{K,i}^c)_i|\theta_{K}^c, z_{K,pop}^c)}$
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Set $z_{K+1,pop}^c = z_{*,pop}^c$ with probability $\alpha$ else $z_{K+1,pop}^c = z_{K,pop}^c$
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**Individual parameters estimation**
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Sample $(z_{*,i}^c)_i$
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute acceptance ratios $\alpha_i = 1 \wedge \frac{q(y_i|\theta_{K}^c, z_{*,i}^c, z_{K+1,pop}^c)q(z_{*,i}^c|\theta_{K}^c, z_{K+1,pop}^c)}{q(y_i|\theta_{K}^c, z_{K,i}^c, z_{K+1,pop}^c)q(z_{K,i}^c|\theta_{K}^c, z_{K+1,pop}^c)}$
13 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;∀$i$, set $z_{K+1,i}^c = z_{*,i}^c$ with probability $\alpha_i$ else $z_{K+1,i}^c = z_{K,i}^c$
14 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute sufficient statistics $(S^c(y_i, \mathbf{z}_{K+1}))_i$
15 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;• **M step**
16 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Update $\theta_{K+1}^c = \underset{\theta \in \Theta}{\operatorname{argmax}} \log q(\mathbf{y}, z_{K+1,pop}^c, (z_{K+1,i}^c)_i, \theta)$
17 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute $\pi_K^c = \frac{1}{n} \sum_i \pi_{K,i}^c$
18 end
19 end

The challenging part is the update of model parameters within each cluster. In the case of a single model, the updates of the parameters $\theta$ of the model are computed in a closed form and are only a function of the total sufficient statistic, which writes $S(\mathbf{y}, \mathbf{z}_{K+1}) = \sum_i S(y_i, \mathbf{z}_{K+1})$ since the observations $(y_i)_{i}$ are supposed independent. When mixtures are involved, the log-likelihood changes to $\log q(\mathbf{y}, \mathbf{z}^c, \theta^c) = \sum_i \log (\sum_c \pi_i^c q(y_i, z_i^c, z_{pop}^c, \theta^c))$. However since the EM computes an expectation of the log-likelihood, we can condition on the true latent class $\pi_i$ of each individual $i$ such that $\pi_i^c = \mathbf{P}(\pi_i = c)$:

$$
\begin{align*}
\mathbf{E}_{\mathbf{z}}\left(\log\left(\sum_c \pi_i^c q(y_i, z_i^c, z_{pop}^c, \theta^c)\right)\right) &= \mathbf{E}_{\pi}\left(\mathbf{E}_{\mathbf{z}}\left(\log\left(\sum_c \mathbf{1}_{\pi_i=c} q(y_i, z_i^c, z_{pop}^c, \theta^c)\right)\middle|\pi_i\right)\right) \\
&= \mathbf{E}_{\pi}\left(\mathbf{E}_{\mathbf{z}}\left(\sum_c \mathbf{1}_{\pi_i=c} \log\left(q(y_i, z_i^c, z_{pop}^c, \theta^c)\right)\middle|\pi_i\right)\right) \\
&= \mathbf{E}_{\mathbf{z}}\left(\sum_c \pi_i^c \log\left(q(y_i, z_i^c, z_{pop}^c, \theta^c)\right)\right)
\end{align*}
$$

Thus we obtain the total sufficient statistics for each cluster $S^c(\mathbf{y}, \mathbf{z}_{K+1}) = \sum_i \pi_i^c S(y_i, \mathbf{z}_{K+1})$ which can be used to update $\theta^c$. This formula is intuitive: each sufficient statistic for each cluster weights the contribution of the individual data by the probability of this individual being in the cluster. This result also proves that the mixture model still belongs to the curved exponential family. Therefore, the convergence of this algorithm is guaranteed by the MCMC-SAEM proof<fnref n="19" />.

<footnote n="19">¹⁹E. Kuhn and M. Lavielle, “Coupling a stochastic approximation version of EM with an MCMC procedure,” fr, *ESAIM: Probability and Statistics*, vol. 8, pp. 115–131, 2004.</footnote>
In the case of the individual parameters mixture model, the main difference is that the sampling of the population parameters is shared between clusters. The maximization of $\bar{\tau}^c$, $\bar{\xi}^c$, $\Sigma_{\tau,\xi}^c$ is done separately for each clusters using the weighted sufficient statistics $S^c(y, z_{K+1})$, but is shared for the rest of the model parameters. The constraint on the $\xi$ to be centered 3.3.3 (for identifiability with $v$) is solved in the M-part by shifting $\tilde{v}$ by the mean of the $\xi_i$ and substracting it from the $\xi_i$. Intuitively we just scale the average speed of the trajectory so that the acceleration factor is still centered around 1.

## Tempered scheme

Even with theoretical guarantees, the convergence might be very slow in practice. The estimation of individual parameters in a non-linear single model may be challenging. Building a mixture on top of it further adds to the difficulty. Practical experiences, described in the results section, show a high reliability on initialization. One of the bottleneck is the amount of cluster regularization contained in the term $p(z_i^c|\theta^c, z_{pop}^c)$ in the likelihood. Until the cluster population parameters stabilize, we do not want to restrict the exploration of the individual parameters space.

We thus propose to use a tempered scheme for the Gibbs sampler, mimicking simulated annealing. Tempered MCMC-SAEM has been shown to converge$^{20}$. In a tempered scheme, inverse temperature comes as a multiplier of the log-likelihood of the model. We realized our model had overly constraining regularization, outweighing the data term in the likelihood. Thus we applied temperature only to the regularization term of the log-likelihood in the sample. In practice, we propose to simply replace $q$ by $\tilde{q} = q^{1/T}$.

For the temperature scheme, we opted for a sinusoidal pattern with a decreasing hull (sine cardinal) following$^{21}$:

$$T(\kappa) = 1 + b \frac{\sin(\kappa)}{\kappa}, \quad \kappa = \Delta + 2\pi \frac{K}{p} \qquad (5.2.3)$$

where $b$ can be seen as the amplitude of the oscillations, $\Delta$ as a phase delay, and $p$ as a period, $K$ is the iteration number. The use of this tempered scheme allows for an alternation between exploratory phases (high temperature) and exploitation phases (low temperature). The values of the hyperparameters were empirically set to $b = 1, \Delta = 0, p = N_{iter}/100$ after several attempts.

## Initialization method

With the same practical considerations, we had to find an initialization method which would improve the odds of convergence towards the global optimum. Random or manual initialization of the model parameters $(\theta_0^c)_c$ was not satisfying, as we will see in the experiments.

We opted for an initialization as close to the clusters as possible. We first fit a disease course mapping model (without mixture) on the whole data. This yields a set of parameters $(\theta_{init}, z_{pop,init}, (z_{i,init}))$. Then we use a Gaussian mixture model (GMM) on the estimated individual parameters $z_{i,init}$. The GMM allows us to identify L clusters within the individual data. The parameters of each mode c in the GMM combined with the population parameters $\theta_{init}$ produce new population parameters for a disease course mapping model which represent the mode: for instance

<footnote>$^{20}$S. Allassonnière and J. Chevallier, “A new class of stochastic EM algorithms. Escaping local maxima and handling intractable sampling,” en, *Computational Statistics & Data Analysis*, vol. 159, p. 107 159, Jul. 2021.</footnote>

<footnote>$^{21}$S. Allassonnière and J. Chevallier, “A new class of stochastic EM algorithms. Escaping local maxima and handling intractable sampling,” en, *Computational Statistics & Data Analysis*, vol. 159, p. 107 159, Jul. 2021.</footnote>
the new $t_0^c$ is the shift of $t_{0,init}$ by the mean value of $\tau_i$ in the mode. This provides initialization parameters for as many disease course mapping models as there are modes in the GMM, and we use these initialization parameters for the clusters of the mixture model.

Please note that the use of a *post hoc* clustering on the estimated parameters is not equivalent to the joint estimation in the mixture model, so this initialization method does not directly give the right cluster parameters. Indeed, we can only produce clusters such that $\forall c, v^c \in \text{Span}(v_{init})$, which is a strong limit to the posterior analysis of a one-class model.

### Simulated data: Multivariate model

In the previous experiments 3.4.1, the results of the mixture DCM model showed consistent performance on a simplified univariate model. However when we add multiple sources, the number of parameters increases drastically and convergence might be more challenging in practice. We generated a dataset as in the 3.4.1 experiment but with two sources and three-dimensional observations. We first fitted a single disease course mapping model on the generated dataset, yielding model parameters $(\theta_{init}, z_{pop,init}, (z_{i,init})_i)$. We then fitted several mixture of DCM models to illustrate our practical improvements. We compared the performance of different initialization methods. In each case, we fitted a model with and without the tempered scheme. Final results are shown in Table 5.1. Each model estimation takes about 5 minutes for 4,000 iterations.

<table><caption>Table 5.1: Performances of models. Single: single disease course mapping model; Random: random initialization; Init: initialization for both clusters at $(\theta_{init}, z_{pop,init})$; GMM: the initialization method described previously using a GMM on $(z_{i,init})_i$; True: perfect initialization at the ground truth parameters.</caption><thead><tr><th>Model</th><th>ROC AUC</th><th>Fit (log q)</th><th>Regularization (log p)</th></tr></thead><tbody><tr><td>Single</td><td>-</td><td>42288</td><td>-3662</td></tr><tr><td>Random non-tempered</td><td>0.50</td><td>25922</td><td>-1439</td></tr><tr><td>Random tempered</td><td>0.50</td><td>36813</td><td>-3622</td></tr><tr><td>Init non-tempered</td><td>0.53</td><td>36523</td><td>-3637</td></tr><tr><td>Init tempered</td><td>0.56</td><td>40034</td><td>-3882</td></tr><tr><td><strong>GMM non-tempered</strong></td><td><strong>0.99</strong></td><td><strong>45787</strong></td><td><strong>-2914</strong></td></tr><tr><td><strong>GMM tempered</strong></td><td><strong>0.99</strong></td><td><strong>46871</strong></td><td><strong>-3068</strong></td></tr><tr><td>True non-tempered</td><td>1.00</td><td>53613</td><td>-2735</td></tr><tr><td>True tempered</td><td>1.00</td><td>53571</td><td>-2280</td></tr></tbody></table>

The table shows that the tempered version allows for a better fit at the cost of regularization overall. Random initialization does not work. "Init" shows a better fit but is not able to cluster the individuals. "GMM" initialization is the best method in absence of better heuristics.

Practical applications in multivariate settings also emphasized the need for a pertinent initialization. Otherwise the mixture of DCM models is prone to diverging.

### 5.2.2 Alternating maximization algorithm

In the Geodesic Bending model, we have a new model parameter: the metric $g$, or more specifically its parametrization. The metric $g$ was optimized via the learning of a diffeomorphism $f$ in a reproducing kernel hilbert space with kernel $k$. The maximization step can be stated as an optimization
problem. If the population and random effects $\theta = (\theta_0, (z_i))$ are estimated, the minimization of the error of reconstruction (in the gaussian noise case) can be performed by solving the following problem:

$$f^* \in \underset{f \in \bar{B}_H(0,c)}{\operatorname{argmin}} L(f), L(f) = \sum_{i,j} |y_{i,j} - \gamma_i(t_{i,j}) - f(\gamma_i(t_{i,j}))|^2 \quad (5.2.4)$$

Thanks to the representer theorem, we have:

$$\forall x \in f^*(x) := \sum_{i,j} k(x_{i,j}, x)w_{i,j}$$

where $x_{i,j} = \gamma_i(t_{i,j})$. With this parametric form for $f^*$, the previous problem (5.2.4) is equivalent to:

$$\underset{W}{\operatorname{minimize}} \quad \|U - K_X W\|^2 \quad (5.2.5a)$$

$$\text{subject to} \quad W^T K_X W \le c_0, \quad (5.2.5b)$$

where the previous variables are vectorized:

$$W := (w_{i,j})_{i,j}, \quad U := (y_{i,j} - \gamma_i(t_{i,j}))_{i,j}, \quad K_X := (k(x_{i,j}, x_{k,l}))_{((i,j),(k,l))}$$

and $c_0$ is the constant ensuring that $f$ is a diffeomorphism.

The minimization of a quadratic function under constraints (5.2.5a) can be solved numerically within a reasonable time as soon as the matrix $K_X$ is not in high dimensions. The total number of visits in a medical cohort can reach up to dozens of thousands and the points $x_{i,j}$ are close to each other, which often results in the matrix being ill-conditioned. To overcome these obstacles, we assume that $y_{i,j} \in [0, 1]^d$ and we create a grid of control points in $[0, 1]^d$ with step-size $\sigma$ if $\sigma$ is the kernel variance (there exists other methods$^{22}$, but this one is more close to the LDDMM framework$^{23}$ and gives satisfying results$^{24}$). We only keep the grid points which are $\sigma$-close to the observations, thus removing the useless ones. We note $N_c$ the number of control points, which can increase exponentially with the dimension according to the previous construction: $N_c(\sigma) \approx \frac{1}{\sigma^d}$. If $d$ is too large, we can choose another method to obtain a reasonable number of control points, a simple way to do so is to subsample the points $(x_{i,j})$. In all the following, the grid method is used because of its practical upsides (no randomness, no risk of ill-conditioning of the matrix $K_X$).

Even though we can compute the argmax for $f$, it would be too costly to do so at every iteration of the MCMC SAEM. Therefore we only maximize $f$ every $n_{MCMC}$ iterations. This initially was called the alternating maximization algorithm, as the procedure alternates between MCMC SAEM iterations that only optimize the base DCM model parameters, and the metric optimization. This metric optimization is then repeated for each diffeomorphic composition. In the original work$^{25}$,

<footnote>$^{22}$F. Liu, X. Huang, Y. Chen, et al., “Random Features for Kernel Approximation: A Survey on Algorithms, Theory, and Beyond,” *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 44, no. 10, pp. 7128–7148, Oct. 2022.</footnote>

<footnote>$^{23}$S. Joshi and M. Miller, “Landmark matching via large deformation diffeomorphisms,” *IEEE Transactions on Image Processing*, vol. 9, no. 8, pp. 1357–1370, Aug. 2000.</footnote>

<footnote>$^{24}$T. Dao and C. D. Sa, “Gaussian Quadrature for Kernel Features,” en,</footnote>

<footnote>$^{25}$S. Gruffaz, P.-E. Poulet, E. Maheux, et al., “Learning Riemannian metric for disease progression modeling,” in *Advances in Neural Information Processing Systems*, M. Ranzato, A. Beygelzimer, Y. Dauphin, et al., Eds., vol. 34, Curran Associates, Inc., 2021, pp. 23 780–23 792.</footnote>
the number of composition is not selected on any precise criterion, but one could simply compute the Bayesian Information Criterion to decide if the added parameters are worth it. The algorithm is described below.

<table><caption>**Algorithm 4:** Geodesics Bending (Alternating maximization algorithm)</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>Initialization: g<sub>init</sub>, θ<sub>init</sub>, N<sub>comp</sub>, k, N<sub>c</sub>, n<sub>MCMC</sub></td></tr><tr><td>2</td><td>g &larr; g<sub>init</sub></td></tr><tr><td>3</td><td>φ &larr; id</td></tr><tr><td>4</td><td>θ &larr; θ<sub>init</sub></td></tr><tr><td>5</td><td>for l = 1 to N<sub>comp</sub> do</td></tr><tr><td>6</td><td>&nbsp;&nbsp;&nbsp;&nbsp;φ<sub>l</sub><sup>0</sup> &larr; id</td></tr><tr><td>7</td><td>&nbsp;&nbsp;&nbsp;&nbsp;for k = 1 to N<sub>c</sub> do</td></tr><tr><td>8</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;g<sub>l</sub><sup>k</sup> &larr; g<sub>l</sub><sup>φ<sub>l</sub><sup>k</sup> &circ; φ</sub></td></tr><tr><td>9</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Run the MCMC-SAEM for n<sub>MCMC</sub> iterations with metric g to estimate θ<sub>l</sub> θ &larr; θ<sub>l</sub></td></tr><tr><td>10</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Solve the optimization problem (5.2.5a) with parameters θ to estimate f<sub>l</sub><sup>k</sup></td></tr><tr><td>11</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;φ<sub>l</sub><sup>0</sup> &larr; id + f<sub>l</sub><sup>k</sup></td></tr><tr><td>12</td><td>&nbsp;&nbsp;&nbsp;&nbsp;end</td></tr><tr><td>13</td><td>&nbsp;&nbsp;&nbsp;&nbsp;φ &larr; φ<sub>l</sub><sup>N<sub>c</sub></sup> &circ; φ</td></tr><tr><td>14</td><td>&nbsp;&nbsp;&nbsp;&nbsp;g &larr; g<sup>φ</sup></td></tr><tr><td>15</td><td>end</td></tr><tr><td colspan="2">Result: return θ, g</td></tr></tbody></table>

We start with a simple metric $g_{init}$ (usually the Euclidean metric or the metric resulting from the push-forward of the logit), and we estimate the mixed-effect parameters $\theta$ with the MCMC-SAEM algorithm. Then, with $\theta$ fixed, we estimate the metric by numerically solving the optimization problem. We repeat these two steps $N_c$ time to obtain one diffeomorphism $\phi_1$ that we compose with the identity. Then we loop on this diffeomorphism learning, and each new one is composed with the previous one. We repeat this loop $N_{comp}$ times. The hyperparameters $k$, $N_{comp}$ and $N_c$ are tuned to make a compromise between loss minimization, computation time and number of parameters.

For $N_{comp}$ compositions, we have $g = g^{\Phi_{N_{comp}}}$ with $\Phi_{N_{comp}} = \phi_{N_{comp}} \circ \dots \circ \phi_1$ where $\phi_i = + \sum_{j=1}^{m} k(x_j^i, \cdot) w_j^i, (x_j^i)_j$ the $m$ control points and $(w_j^i)_j$ their associated weights. This structure resembles the deep neural networks considering $m$ as the width and $N_{comp}$ as the depth. The total time complexity is $O(n_{MCMC}dN_c m N_{comp}^3)$ which reduces in practice the choice of $N_{comp}$.

In the end the estimation algorithm for Geodesic Bending is a MCMC-SAEM with us taking some liberties in the frequency of the maximization rule.

### 5.2.3 Gradient maximization step

In the piece-wise geodesic model we introduced a new function, the treatment function $f$, which we will generally parametrize with $\theta_f$. These parameters become model parameters if we maximize them during the EM, or population parameters if we prefer to sample them. As a general rule of thumb, we prefer avoid sampling parameters as a stochastic search is less optimal than a concrete formula. A model with a higher number of sampled parameters will require more overall MCMC SAEM steps to converge. This is particularly true when the choice of parametrization of $f$ involves
a high number of parameters, such as the RKHS formulation used in 4.3.2. Sampling with a simple Metropolis-Hastings algorithm in such high dimensional cases could be complicated due to the energetic landscape being mostly flat.

Then the question becomes: can we solve the maximization problem ?

$$ \theta_f^* = \underset{\theta_f}{\operatorname{argmax}} Q(\theta_f, \mathbf{y}, \theta, \mathbf{z}) \qquad (5.2.6) $$

Let $q_{prior}$ be a Gaussian prior added to the weights $\theta_f$. If we are in the continuous setting, with a Gaussian noise model, maximizing the likelihood becomes equivalent to the regularized least-square problem:

$$ \theta_f^* = \underset{\theta_f}{\operatorname{argmin}} \sum_{i,j,k} \frac{1}{2\sigma_k^2} (y_{ijk} - \eta_k(\psi_i(t_{ij}); \mathbf{w}_i))^2 + \lambda ||\theta_f||^2 \qquad (5.2.7) $$

We actually used a mixed L1 and L2 regularization in practice, which is equivalent to a mix between a Gaussian and a Laplace prior.

The problem (5.2.7) is not trivial. One simple solution would be to be able to compute the gradient of the objective function and to use it with a gradient descent for example. However the gradient cannot be formulated explicitly. Here is why. After every j-th breakpoint of individual $i$ $b_{i,j}$, the reference curve follows a trajectory that we can generally write with the Riemannian Exponential as : $t \mapsto \operatorname{Exp}_{p_{i,j}}(tv_{i,j})$. The next point $p_{i,j+1}$ is therefore equal to $\operatorname{Exp}_{p_{i,j}}((b_{i,j+1} - b_{i,j})v_{i,j})$, and the next velocity $v_{i,j+1}$ is $f(p_{i,j+1}, \tilde{v}_{i,j+1}, X_{i,j+1})$ where $\tilde{v}_{i,j+1}$ is the derivative at point $\mathbf{p}_{i,j+1}$, and thus a function of $p_{i,j}$ and $v_{i,j}$. So if we assume $f$ to be a function of the position $p_{i,j+1}$, and of the velocity $v_{i,j+1}$, then the velocity after the next break $v_{i,j+1}$ depends on $\theta_f$, $p_{i,j}$ and $v_{i,j}$. And the breakpoint after that will add another layer of complexity by applying $f$ one more time on a position that already depends on $f$. And so on and so forth.

Every model prediction at time $t_{ij}$ will thus depend on the number of breakpoints before $t_{ij}$. For every $n$ breakpoints before $t_{ij}$ the function $f$ would have been applied $n$ times in the graph of operations leading to $\eta_k(\psi_i(t_{ij}); \mathbf{w}_i)$. This is reminiscent of recurrent neural networks where the same function is applied several times in a row. In the same way we will use backpropagation to compute the gradient of the objective with respect to the treatment function parameters $\theta_f$, and then perform gradient descent to find a local optimum.

Performing a full gradient descent would be way too long at each iteration of the MCMC-SAEM. We thus relax the maximization step by only performing one gradient descent step during the maximization. This only improves the expected likelihood $Q$ instead of maximizing it, but this is still sufficient to guarantee the convergence of the algorithm. This version of the EM is called Generalized EM, and is proven to converge in$^{26}$. As for Geodesic Bending, the optimization operation in the M-step is still costly so we only perform the gradient step every $n_{MCMC}$ MCMC SAEM iterations.

Computing the gradient via backpropagation however suffers from the same issues as in the field of Deep Learning with vanishing gradients. This is the case when using non-linear "link" functions, for instance in the logistic curves case. The logistic has almost a null gradient on both extremities, which leads to these vanishing gradient issues. This problem is the reason why all the experiments in the piece-wise geodesic model section 4.3.2 were performed in the linear case.

<footnote>$^{26}$R. Harpaz and R. Haralick, "TR-2006001: The EM Algorithm As a Lower Bound Optimization Technique," Computer Science Technical Reports, Jan. 2006.</footnote>
Another possibility that we did not explore but which is very promising is to allow the model to learn the breakpoints $b_{ij}$ by setting them as latent variables, and sampling them during the MCMC-SAEM.

### 5.2.4 Discussion

The MCMC SAEM is a powerful algorithm to estimate model parameters while also providing sampling for the latent variables. This proves to be very useful as we are in a Bayesian setting, the Markov chain samplers converge towards the posterior distribution of the latent variables, thus we can easily access it.

Even if the Metropolis-Hastings algorithm guarantees the sampler to converge towards the target distribution, in practice we can run into computational issues. We could easily end up in a problematic case where the proposal is never accepted due to the high dimension of $z$, even though we have an adaptive variance proposal. Indeed the chance to sample in the right direction to increase the log-likelihood becomes exponentially smaller in high dimensions, making the Metropolis-Hastings almost useless in practice. This is why we use the Gibbs variation, where each coordinate of the latent variables is sampled separately. However the process can be very slow, as the progress is still dictated by a random walk search.

The Metropolis-Hastings is simple to implement but has long been subclassed by other samplers in the field of MCMC sampling. Namely the Hamiltonian Monte-Carlo (HMC) sampler<fnref n="27" /> uses a dynamical approach to speed-up sampling using gradient information. The state-of-the-art sampler currently is a variant of the HMC called the No-U-turn sampler<fnref n="28" />. This last one is more complicated to implement.

Another sampler using the gradient to improve sampling is the Metropolis Adjusted Langevin Algorithm (MALA) and its Riemannian improvements. In the MALA, the proposal is informed by the local gradient of the parameters to improve the convergence towards the mode of the distribution. This proves to be very useful in high dimensional settings where the random walk has a hard time finding its way to the mode. In<fnref n="29" /> the authors propose to improve both the MALA and the HMC samplers by using the natural gradient instead of the standard gradient. The natural gradient corresponds to the gradient in a Riemannian manifold where the metric is the Fisher information matrix of the parameters we wish to sample. Both samplers show an increased convergence rate compared to their non-Riemannian counterpart. Since we already computed the gradients and the Fisher information of the DCM models 15, it would be a simple future improvement to implement these samplers instead of the Metropolis-Hastings within Gibbs.

<footnote n="27">²⁷M. Betancourt, *A Conceptual Introduction to Hamiltonian Monte Carlo*, Jul. 2018.</footnote>

<footnote n="28">²⁸M. D. Homan and A. Gelman, “The No-U-turn sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo,” *The Journal of Machine Learning Research*, vol. 15, no. 1, pp. 1593–1623, Jan. 2014.</footnote>

<footnote n="29">²⁹M. Girolami and B. Calderhead, “Riemann manifold Langevin and Hamiltonian Monte Carlo methods,” en, *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, vol. 73, no. 2, pp. 123–214, 2011.</footnote>
Conclusion and perspectives

The work presented in this thesis has pushed the boundaries of a previous disease progression model, namely the Disease Course Mapping model. Building upon a geometric model designed to disentangle spatial and temporal effects from a longitudinal dataset, we provided some necessary extensions to apply the model to a wider range of neurodegenerative diseases' cohorts. The model is able to provide insights at the global population level, thereby improving our understanding of the shared patterns of progression across all the individual trajectories. Additionally the hierarchical structure allows for a personalization of the trajectories at an individual level, which is crucial for prognosis or clinical trial enhancement.

We introduced Geodesic Bending as a method to palliate the rigidity in the original DCM model. This technique, inspired from the large deformation diffeomorphic metric mapping field, allows the model to adapt the shape of the trajectory to the data by learning an implicit Riemannian metric. In the applications to Alzheimer's disease, we have shown that this method outperforms the state-of-the-art in predicting the progression of cognitive decline. The average trajectory learned by the model was more in adequacy with the known progression of the disease compared to the logistic curves model. More interestingly the deformations broke the monotonic paradigm of the DCM model. The trajectory remains injective in the multivariate setting, but each biormarker's progression profile can now be learned without the monotonous assumption. This greatly expands the application possibilities of disease progression modelling while retaining the benefits of the interpretability unfolding from the spatio-temporal description.

We then transcended the continuous framework that was inherent to the design of the DCM model with the trajectories on the Riemannian manifold. We used the logistic curves formulation of the DCM to describe the evolution of a latent progression while a Bernoulli and a cumulative logit model described the discrete observations. Even though the model improvement is very simple, the implications in practice are manifold. Symptoms occurrence, disease staging and rating scales are now available in the DCM model. This opens the door to symptom prediction or staging diagnosis. Furthermore, the possibility to model items of a psychometric scale directly instead of the aggregated score allows to describe the disease at the finer-grained level, yielding even more individual variability.

Across our experiments with neurodegenerative data we have felt the urge to improve the model to deal with the heterogeneity of patients' profiles. The random effects of the mixed-effect model structure only account for a natural variability, *random* as implied by the name. However this variability can stem from non-random sources, which can be described by external covariates or by disease subtypes. Focusing on the latter we introduced a mixture option in the DCM model, either at the individual parameters level or at the global progression level. Both rely on the use of latent classes for the individuals. The mixture model aims at jointly estimating these latent classes
while learning the average disease progression of each class. The two variants each have their own preferred use case, allowing to uncover new patterns of variability across the data.

In the last part of our work we focused on what we believe to be the next area of research for disease progression models. Treatment modelling's importance is going to surge as new therapeutic approaches are progressively delivering results. Applications of disease progression models are plural, including in clinical trials. We proposed one such application, which paves the way for the design of clinical trials with only one treated arm by leveraging a synthetic control arm. In our proposed framework, a disease progression model is used to predict the natural disease progression, i.e. under no treatment, allowing a treatment model to regress the difference with the observed progression of the treated patient. We showcased our approach using Parkinson's disease data to analyse the effect of the dopaminergic treatment over motor symptoms, quality of life and the levels of dopamine transporter in the putamen. Finally we introduced the Piece-wise Geodesic model, a Riemannian based model to account for disease-modifying treatment effect over the trajectories. This model has been tested on synthetic data and still awaits the proper dataset to be applied on real-life data.

All of the design of those models and extensions is subject to the estimation condition. We had to design them so that identifiability was ensured, while making sure that the new parameters could be learned properly. In the estimation chapter, we explained our modelling choices from the point of view of the MCMC-SAEM algorithm. We provided insights into why we chose certain parameters to be latent variables, which ones we sampled and which ones were directly maximized, or simply updated following a gradient step in the case of the Piece-wise Geodesic model. We identified the practical issues preventing convergence and discussed about empirical solutions.

Our improvements to the DCM framework all contribute to design better, larger, more comprehensive models for neurodegenerative diseases. By doing so, we indirectly increase our understanding of disease progression, we embrace the heterogeneity and help in the development of a cure. This has been highlighted in our many applications, where we tried to stress the clinical implications of the results. During this thesis we have striven to communicate with clinicians, on the occasion of conferences 15 or by exchanging with neurologist from our institute, so that they benefit from our work and that we could improve our model based on their feedback.

## Limitations and perspectives

Precise limitations have been highlighted at the end of each dedicated chapter. We will try to summarize them:

* Geodesic Bending: the flexibility of the model comes at a higher computational cost. The method is more data-hungry due to a higher risk of overfitting on few samples.

* Binary model: models using purely symptom features do not perform well in terms of prediction, the model should include for informative features, i.e. continuous markers, to guide the trajectory.

* Ordinal model: the model also requires more data, specifically all the levels of the ordinal items should be observed. As the ordinal models tend to use more features, their complexity increases, especially computationally.

* Mixture model: estimation of multivariate mixture models remains complex, and convergence is hard to reach in some cases. The use of better samplers should improve in this matter.
* Additional treatment effect framework: the treatment effect estimation relies on the predicting accuracy of the disease progression model. However training such a model on pre-treatment observations only is usually prone to fail long-term predictions. The solution would be to use another cohort to learn the natural disease history

* Piece-wise Geodesic model: the current parametrization makes the model very data dependent. The model also needs to be confronted to real-life data.

The main challenge often resides in the estimation of those methods. As we increase the complexity more data is required to properly estimate the models. Some sparsity in the parametrization of the learned function, for instance in Geodesic Bending or in the Piece-wise Geodesic model, could help regarding this issue. We believe that one major improvement could be the implementation of better samplers.

Some further investigation into the compatibility of the models should be done. For instance crossing Geodesic Bending with the mixture model could lead to learning different shapes of progression for the subtypes. We mentioned that symptoms were difficult to model because many were influenced by treatments, so it would make sense to use the treatment model to assess what would have been the progression without treatment and how much the treatment affected the symptom. Finally we hope that our first steps in the field of treatment modelling will be followed as this will prove to be useful in the years to come.
<figure><img src="image_55.png" /></figure>
# Publications and scientific communication

## Scientific articles

### Journal articles

* Pierre-Emmanuel Poulet, Stanley Durrleman. Multivariate disease progression modelling with longitudinal ordinal data. Statistics in Medicine, 2023; 1- 20. doi: 10.1002/sim.9770

* Moulaire, P., Poulet, P.E., Petit, E., Klockgether, T., Durr, A., Ashisawa, T., Tezenas du Montcel, S. (2023), Temporal Dynamics of the Scale for the Assessment and Rating of Ataxia in Spinocerebellar Ataxias. Mov Disord, 38: 35-44. https://doi.org/10.1002/mds.29255

### Conference proceedings papers

* Poulet, PE., Durrleman, S. (2021). Mixture Modeling for Identifying Subtypes in Disease Course Mapping. In: Feragen, A., Sommer, S., Schnabel, J., Nielsen, M. (eds) Information Processing in Medical Imaging. IPMI 2021. Lecture Notes in Computer Science(), vol 12729. Springer, Cham. https://doi.org/10.1007/978-3-030-78191-0_44. Poster presentation at the virtual conference event

* Samuel Gruffaz, Pierre-Emmanuel Poulet, Etienne Maheux, Bruno Jedynak, Stanley Durrleman. Learning Riemannian metric for disease progression modeling. NeurIPS 2021 - Thirty-fifth Conference on Neural Information Processing Systems, Dec 2021, link

* (in preparation) Pierre-Emmanuel Poulet, Bruno Jedynak, Stanley Durrleman. A Bayesian symptomatic treatment model using latent disease progression, 2023

## Conference communications

* Pierre-Emmanuel Poulet, Jean-Christophe Corvol, Stanley Durrleman. A latent mixed-effects model for longitudinal categorical data in Parkinson disease. CompAge 2020 - Computational approaches for ageing and age-related diseases 2020, Sep 2020, Paris, France. Poster presentation at the virtual conference event
* Pierre-Emmanuel Poulet, Stanley Durrleman, Jean Christophe Corvol. Predicting symptom onset in Parkinson's disease with latent mixed-effect model. AD/PD 2021 - 15th International Conference on Alzheimer's Parkinson's Diseases, Mar 2021, Barcelone, Spain. Oral presentation at the virtual conference event

* Pierre-Emmanuel Poulet, Bruno Jedynak, Stanley Durrleman. Learning treatment effect in neurodegenerative diseases with a Bayesian mixed-effect model. 43rd Annual Conference of the International Society for Clinical Biostatistics (ISCB), Aug 2022, Newcastle upon Tyne, United Kingdom. Poster presentation, award for best poster.

* Pierre-Emmanuel Poulet, Bruno Jedynak, Zoltan Mari, Stanley Durrleman. Non-parametric modeling for dopaminergic treatment effect in Parkinson's disease. AD/PD 2023 - 17th International Conference on Alzheimer's Parkinson's Diseases, Mar 2023, Goteborg, Sweden. Poster presentation

* (in preparation) Pierre-Emmanuel Poulet, Bruno Jedynak, Stanley Durrleman. Piece-wise Geodesic: a longitudinal disease-modifying treatment effect model. 44th Annual Conference of the International Society for Clinical Biostatistics (ISCB), Aug 2023, Milan, Italy. Poster presentation.

Teaching

* Teaching mission, Sorbonne university, Computational department, Jussieu. List of courses:

- 2020: "Elements of programming" License 1 (equivalent: bachelor 1) introductory course to Python, practical sessions

- 2021: "Elements of programming" License 1, practical sessions

- 2022: "Elements of programming" License 1, group classes and practical sessions

- 2022: "Data science" License 1 introductory course to data science, practical sessions

* "Disease course mapping with longitudinal data" practical session workshop at AI4health Virtual Winter school, with Arnaud Valladier, Igor Koval, Etienne Maheux, Juliette Ortholand and Stanley Durrleman, Jan 2021

* "Disease course mapping with longitudinal data" practical session workshop at AI4health Winter school, with Némo Fournier, Etienne Maheux and Juliette Ortholand, Jan 2022

* "Disease course mapping with longitudinal data" practical session workshop at AI4health Summer school, with Sofia Kaisaridi and Juliette Ortholand, Jul 2023
Software development

Our code contributions are all published in the Python open-source library Leaspy:
https://leaspy.readthedocs.io/en/stable/

The library is still expanding with many new variants, including all the work presented in this
manuscript. The Gitlab link to the code is:
https://gitlab.com/icm-institute/aramislab/leaspy.

The Leaspy library was released in 2017 under the GNU GPLv3 license, and benefits from the work of various contributors. All the code is peer-reviewed by the developers before being released in the stable version. During my thesis I contributed as a developer of the library, with commits corresponding to the discrete data models in the **master** branch, the collaborative implementation of the **GeodesicBending** branch with Samuel Gruffaz, the whole **mixture_model** branch and the whole **treatment_modelling** branch.

The stable version of Leaspy (i.e. the pip release and the master branch of the Git) already
include the binary and ordinal DCM models, which I developed and implemented. A minimal
example of code usage is given below:

```python
1 # load packages
2 import pandas as pd
3 from leavy import Leaspy, Data, AlgorithmSettings
4
5 # load data
6 df = pd.read_csv("path_to_my_data.csv", index_col = ["ID", "TIME"])
7
8 # Data format:
9 # ID      TIME      Y1  Y2
10 # 01      60.       0.  1.
11 # 01      65.       ... ...
12
13 leaspy_data = Data.from_dataframe(df)
14
15 # Instantiating a logistic model
16 leaspy = Leaspy("logistic", noise_model = "gaussian") # for continuous data
17 # or
18 leaspy = Leaspy("logistic", noise_model = "bernoulli") # for binary data
19 # or
20 leaspy = Leaspy("logistic", noise_model = "ordinal") # for ordinal data
21
22 # Calibration
23 algoSettings = AlgorithmSettings('mcmc_saem', n_iter = 1024)
24 leaspy.fit(leaspy_data, algoSettings)
25
26 # Personalization
```
27 settingsPersonalization = AlgorithmSettings("scipy_minimize", use_jacobian=True)
28 individual_parameters = leaspy.personalize(leaspy_data, settingsPersonalization)
29
30 # Prediction
31 values = leaspy.estimate(timepoints, individual_parameters)

An introduction to Leaspy with hands-on tutorial notebooks can be found at:
https://gitlab.com/icm-institute/aramislab/disease-course-mapping-tutorials/-/tree/master/challenges

Our personal code for Geodesic Bending, mixture models and Piece-wise geodesic model is still
waiting for a stable release. Leaspy is undergoing a deep refactoring for a second version as of the
time of writing this thesis, and will probably be released in the year following this work.
# Appendix

## Riemannian geometry

This section is an introduction to the concepts of Riemannian geometry used in this thesis. We don't claim it to be a comprehensive overview of Riemannian geometry, for this we refer the reader to the excellent book by Manfredo Do Carmo: *Riemannian Geometry*<sup>30</sup>. Riemannian geometry is a powerful tool to study structured manifolds. The applications are ubiquitous: for instance medical images only span a subset of all the possible images, and such subset might be of a much lower dimension with a specific geometry. We will therefore often see the manifold $\mathcal{M}$ as an open subset of $\mathbb{R}^D$. We will first provide all the definitions required to reach the notion of Riemannian manifold.

**Definition 1 (Manifold).** A topological space $\mathcal{M}$ is a $d$-dimensional manifold iff, $\forall p \in \mathcal{M}, \exists U(p)$ a neighbourhood and a homeomorphism $x: U(p) \to x(U(p)) \subset \mathbb{R}^d$.

$(U, x)$ is called a local coordinate chart. A collection of charts which covers $\mathcal{M}$ is called an atlas. A manifold is called smooth iff for any two charts $(U, x), (V, y)$ in an atlas, $y \circ x^{-1}$ is $C^\infty$. The dimension $d$ of an embedded manifold is smaller than the space in which it is embedded $d \le D$.

Intuitively, a smooth manifold of dimension $d$ is a topological space where each point can be assigned coordinates, and these coordinates are differentiable.

**Definition 2 (Tangent space).** The tangent space at $p$ in a smooth manifold $\mathcal{M}$ is $\{v/\exists \text{acurvec}: (-\epsilon, \epsilon) \to \mathcal{M}/v = \dot{c}(0) \text{and } p = c(0)\}$. It is noted $TM_p$.

The tangent space is therefore the vector space of all possible derivatives of curves passing through $p$. This space is distinct in each point of the manifold, but if we write $(x_1, ..., x_d)$ the local coordinates in the neighbourhood of $p$, we can define the basis $(\frac{\partial}{\partial_1}(p), ..., \frac{\partial}{\partial_d}(p))$ of the tangent space $TM_p$.

**Definition 3 (Vector field).** A vector field $X$ is a smooth function mapping points to vectors in the tangent space $X: p \in \mathcal{M} \mapsto X(p) \in TM_p$.

In a local chart coordinates, if the coordinates of a point $p$ are $(x_1, ..., x_d)$, we can write the vector field as:

$$X(p) = \sum_{i=1}^{d} X(x_1, ..., x_d) \frac{\partial}{\partial_i}(p)$$

We can now define:

<footnote n="30"><sup>30</sup>M. Do Carmo, *Riemannian Geometry*, en.</footnote>
**Definition 4 (Riemannian manifold).** If $\forall p \in \mathcal{M}, g_p(\cdot, \cdot)$ is a scalar product on $TM_p$ and : $p \to g_p$ is smooth, then $(\mathcal{M}, g)$ is a Riemannian manifold. $g$ is called the Riemannian metric.

One way to construct a metric is to use push-forward metrics, as is done in the Geodesic Bending model:

**Definition 5 (Pushforward metrics).** Provided $(\mathcal{M}, g)$ a Riemannian space, $N$ a manifold and $\phi: \mathcal{M} \to N$ a $C^1$ diffeomorphism, we can equip $N$ with the Riemannian metric $g^\phi$ defined as:

$$ \forall p \in N, \forall w, v \in TN_p, g_p^\phi(w, v) = g_{\phi^{-1}(p)}(d\phi^{-1}(p).w, d\phi^{-1}(p).v) $$

The use of push-forward metrics allows to carry a metric from a manifold to another through a diffeomorphism. In the most simple case the starting manifold is the Euclidean space equipped with the Euclidean metric, from which we simply apply a diffeomorphism.

An important notion in Riemannian geometry is the geodesic, which is the "shortest" path (in the sense of the Riemannian metric) between two points of the manifold. The interest of defining push-forward metrics is that the geodesics then are obtained by simply applying the diffeomorphism to the geodesics in the starting manifold.

Let's define the geodesics:

**Definition 6 (Distance and geodesic).** Given a Riemannian manifold, we define the Riemannian distance $d$ such that

$$ \forall p_1, p_2 \in \mathcal{M}, d(p_1, p_2) = \inf_{\gamma: [0,1] \to \mathcal{M}, \gamma(0)=p_1, \gamma(1)=p_2} \int_0^1 \sqrt{g_{\gamma(t)}(\dot{\gamma}(t), \dot{\gamma}(t))} dt \quad (.0.8) $$

If $\gamma^*$ minimize the previous problem, $\gamma^*$ is called a geodesic.

In the Euclidean space, geodesics are straight lines such that $\ddot{\gamma} = 0$. However we cannot yet compute second order derivatives in the Riemannian manifold since the tangent spaces are different in every point. Indeed, if we have a trajectory $\gamma:] -1, 1[ \to \mathcal{M}$, we can compute the first order derivative for each time $\dot{\gamma}(t) = \lim_{\epsilon \to 0} \frac{\gamma(t+\epsilon)-\gamma(t)}{\epsilon}$ and this element belongs to the tangent space $TM_{\gamma(t)}$ by definition. However we cannot compute $\ddot{\gamma}(t) = \lim_{\epsilon \to 0} \frac{\dot{\gamma}(t+\epsilon)-\dot{\gamma}(t)}{\epsilon}$ since $\dot{\gamma}(t+\epsilon)$ and $\dot{\gamma}(t)$ belong to two possibly different vector spaces. In Euclidean geometry this is not an issue as the tangent spaces are all the same.

This is where the notion of connection comes into play:

To do it, we reduce the problem to the vectors of the basis of each tangent space $(\partial_i(p))$. A connection defines how the basis is varying locally. Precisely:

**Definition 7 (Connection).** An affine connection $\nabla$ maps two vector fields $v : p \in \mathcal{M} \to v(p) \in TM_p$, $w : p \in \mathcal{M} \to w(p) \in TM_p$ to a single vector field $\nabla_v w$ with the following properties:

1. [Left-linearity] $\forall X_1, X_2$ and $Y$ vector fields:

$$ \forall f, g \in C^\infty(\mathcal{M}, \mathbb{R}), \forall p \in \mathcal{M}, \nabla_{fX_1+gX_2} Y(p) = f(p)\nabla_X Y(p) + g(p)\nabla_{X_2} Y(p) $$

2. [Right-Leibniz-linearity] $\forall X, Y$ f-vector fields:

$$ \forall f \in C^\infty(\mathcal{M}, \mathbb{R}), \forall p \in \mathcal{M}, \nabla_X(fY)(p) = f(p)\nabla_{X_1}Y(p) + X(f)(p)Y(p) $$

where $X(f)(p) = \lim_{\epsilon \to 0} \frac{f(\gamma(\epsilon))-f(p)}{\epsilon}$ and $\gamma:] -1, 1[ \to \mathcal{M}$ is a differentiable curve such that $\gamma(0) = p, \dot{\gamma}(0) = X(p)$.
3. [R-right-linearity] $\forall X, Y_1, Y_2$ vector fields:

$$ \forall a, b \in \mathbb{R}, \forall p \in \mathcal{M}, \nabla_X(aY_1 + bY_2)(p) = a\nabla_X Y_1(p) + b\nabla_X Y_2(p) $$

The connection can be defined by its action on the basis $(\partial_i(p))$. If we write $\nabla_{\partial_i(p)} \partial_j(p) = \sum_k \Gamma_{i,j}^k(p) \partial_k(p)$, the coefficients $\Gamma_{i,j}^k(p) \in \mathbb{R}$ are characterizing the connection thanks to the linearity properties. They are called the Christoffel's symbols.

We can now compute the second order derivative by considering that $\nabla_{\dot{\gamma}}\dot{\gamma}$ is the analogue of $\ddot{\gamma}$. There exists a unique specific connection such that the solutions of $\nabla_{\dot{\gamma}}\dot{\gamma} = 0$ are exactly the geodesics. This connection is called the Levi-Civita connection, and it depends on the metric $g$ of the manifold. We will admit these results here. More details on this connection in the book<fnref n="31" />.

We decompose the coordinates of $\dot{\gamma}$ as $\dot{\gamma}(t) = \sum_i \dot{\gamma}_i(t)e_i(\gamma(t))$. Then:

$$ \nabla_{\dot{\gamma}} \dot{\gamma}(\gamma(t)) = \sum_k \left( \ddot{\gamma}_k(t) + \sum_{i,j} \Gamma_{i,j}^k(\gamma(t)) \dot{\gamma}_i(t) \dot{\gamma}_j(t) \right) e_k(\gamma(t)) $$

where $\ddot{\gamma}_k(t) = \lim_{\epsilon \to 0} \frac{\dot{\gamma}_k(t+\epsilon) - \dot{\gamma}_k(t)}{\epsilon}$.

In the end the geodesics equation can be reformulated as:

$$ \ddot{\gamma}_k(t) + \sum_{i,j} \Gamma_{i,j}^k(\gamma(t)) \dot{\gamma}_i(t) \dot{\gamma}_j(t) = 0 \qquad (.0.9) $$

This second order Ordinary Differential Equations (ODE) has solutions under some regularity conditions on the Christoffel's symbols. Geodesics can thus be reduced to the Cauchy conditions in $0 (\gamma(0) = p, \dot{\gamma}(0) = v)$.

We can now define an essential tool of Riemannian geometry:

**Definition 8** (Riemannian Exponential). The Riemannian Exponential is defined for $p \in \mathcal{M}$, $t_0 \in \mathbb{R}$ and $v \in T\mathcal{M}_p$ as the mapping:

$$ \mathrm{Exp}_{p,t_0}(v)(\cdot) : \begin{cases} I \to \mathcal{M} \\ t \mapsto \mathrm{Exp}_{p,t_0}(v)(t) \end{cases} $$

where $\mathrm{Exp}_{p,t_0}(v)(t)$ is the point reached at time $t$ by the geodesic starting at time $t_0$ at point $p$ with velocity $v$.

We will often reduce the Exponential to an abusive notation $p(v)$ where $t_0 = 0$ and $t = 1$ for simplification. This operation of starting at point $p$ and following the direction $v$ for a unit time is called *shooting*.

We say that a Riemannian space is complete if the map: $p: TM \to \mathcal{M}$ is well-defined for all point $p$ in the manifold $\mathcal{M}$ (the ODE admits solutions defined on the open interval $\mathbb{R}$). Note that in the Euclidean case $p(v) = p + v$.

In order to finally reach the last notion of this section, which is the Exp-parallelization, we require one last tool: parallel transport. It will enable us to define parallel vectors in tangents space along a trajectory according to a given connection $\nabla$.

<footnote n="31">³¹M. Do Carmo, *Riemannian Geometry*, en.</footnote>
<figure><img src="image_56.png" /><figcaption>Figure 1: Left: Sketch of the Exp-parallellisation in an abstract manifold. Right: Exp-parallelisation for the unit square as a manifold with the metric $g_{(p_1,p_2)}((x_1,x_2), (y_1,y_2)) := \frac{x_1y_1}{p_1^2(1-p_1)^2} + \frac{x_2y_2}{p_2^2(1-p_2)^2}$. The average trajectory in red and an individual trajectory in blue or purple. Reproduced from³²</figcaption></figure>

**Definition 9 (Parallel transport).** Let $\gamma: t \mapsto \gamma(t)$ be a smooth curve on $M$ and $X$ a vector field on $\gamma(t)$. The vector field $X$ is said to be parallel along $\gamma(t)$ iff $\nabla_{\gamma(t)}X = 0$.

The parallel transport is unique: given a curve $\gamma$ and a vector $v_0$ tangent to $\gamma(t_0)$, there exists a unique vector field $X$ parallel along $\gamma$ such that $X(\gamma(t_0)) = v_0$. We denote $P_{\gamma,t_0,t}(v_0)$ the isometry mapping $v_0$ to $X(\gamma(t))$.

Remark that in the Euclidean space $X(\gamma(t)) = v_0$ (parallel vector are colinear).

**Definition 10 (Exp-parallelization).** Let $(M, g)$ be a geodesically complete Riemannian manifold equipped with its Levi-Civita connection $\nabla^g$. Let $\gamma: I \to M$ be a curve on $M$, $t_0 \in I$ and $w \in T\gamma(\cup_t), w \neq 0$. The Exp-parallelization of $\gamma$ according to $w$ at $t_0$ is the curve:

$$ \eta_w(\gamma, \dot{\gamma}) : \begin{cases} I \to M \\ t \mapsto \mathrm{Exp}_{\gamma(t_0)}(P_{(\gamma,t_0,t)}) \end{cases} $$

Intuitively, the curve $\eta_w(\gamma, \dot{\gamma})$ is a shift from $\gamma$ in the direction $w$ at $t_0$. At $t_0$, the point constructed is simply by shooting with $w$ starting from $\gamma(t_0)$. For the rest of the curve, the only thing changing is that we don't take directly $w$ (as it lives in $TM_{\gamma(t_0)}$), but the result of the transport of this vector $w$ parallel to the curve $\gamma$.

# Disease Course Mapping formulas

## Introduction

This section is a summary of all the formulas for logistic disease course mapping (DCM) model: basic formulation of the model, parameters introduction, gradients and finally Fisher information
matrix. We provide the formulas as they are implemented, which includes some simplifications over the theoretical model introduced in<fnref n="33" />.

# Model description

## Observations and dataset

The input of the model is a longitudinal dataset, for which we introduce the following notations :

* a number of patients N, that will be indexed with i,

* for each patient a certain number nᵢ of visits, that will be indexed with j,

* for each visit of a patient, a time tᵢⱼ is associated and a set of features yᵢⱼₖ with 1 ≤ k ≤ d being the index of the feature

## Model equation

The framework has been first developed in<fnref n="34" /> in a general formulation.

The age of the patient is mapped onto the common disease timeline through time reparametrization ψᵢ(t):

$$ \psi_i(t) = \alpha_i(t - \tau_i) \qquad (.0.10) $$

where αᵢ is the acceleration factor of patient i and τᵢ is its time-shift. The acceleration is reparametrized as αᵢ = eξᵢ to enforce its positivity. The theoretical formulation is ψᵢ(t) = αᵢ(t - t₀ - τᵢ) + t₀ with τ = 0, but this is equivalent to equation (.0.10) above with t₀ = τ while the +t₀ term is absorbed by other parameters.

The trajectory is then obtained by computing an exp-parallelization of a geodesic in the Riemannian manifold. In the logistic case the geodesic is:

$$ \gamma_k(t) = \left(1 + \left(\frac{1}{p_k} - 1\right) \times \exp\left(-\frac{v_k}{p_k(1-p_k)}(t-t_0)\right)\right)^{-1} \qquad (.0.11) $$

where pₖ and vₖ are the k-th coordinate of the position and velocity (derivative of the curve) at time t₀ = τ of the average trajectory. After exp-parallelization you get:

$$ \forall k \in [|1, d|], \quad \eta_k(\psi_i(t_{ij}); \mathbf{w}_i) = \left(1 + \left(\frac{1}{p_k} - 1\right) \exp\left(-\frac{v_k \psi_i(t_{ij}) + w_{ik}}{p_k(1-p_k)}\right)\right)^{-1} \qquad (.0.12) $$

with wᵢₖ being the k-th coordinate of the space-shift wᵢ. Since we will use this notation very often, we introduce the abbreviation ηₖ(ψᵢ(tᵢⱼ); wᵢ) = η which we will call the "model".

The space-shifts are decomposed along Ns sources (0 ≤ Ns ≤ d - 1):

<footnote n="33">³³J.-B. Schiratti, S. Allassonnière, A. Routier, et al., "A Mixed-Effects Model with Time Reparametrization for Longitudinal Univariate Manifold-Valued Data," en, in *Information Processing in Medical Imaging*, S. Ourselin, D. C. Alexander, C.-F. Westin, et al., Eds., ser. Lecture Notes in Computer Science, Cham: Springer International Publishing, 2015, pp. 564–575.</footnote>

<footnote n="34">³⁴J.-B. Schiratti, S. Allassonnière, A. Routier, et al., "A Mixed-Effects Model with Time Reparametrization for Longitudinal Univariate Manifold-Valued Data," en, in *Information Processing in Medical Imaging*, S. Ourselin, D. C. Alexander, C.-F. Westin, et al., Eds., ser. Lecture Notes in Computer Science, Cham: Springer International Publishing, 2015, pp. 564–575.</footnote>
$$w_i = As_i \qquad (.0.13)$$

where $A \in \mathbb{R}^{d \times N_s}$ is called the mixing matrix and is composed of $N_s$ columns orthogonal to $v_0$ for identifiability. This is done by computing an orthogonal basis $(B_m)_{1 \le m \le d-1}$ of $\text{Span}(v)^\perp$ with the Householder method. Then the columns of $A$ are linear combinations of this basis $B$:

$$A_{:,l} = \sum_{m=1}^{d-1} \beta_{ml} B_m \qquad (.0.14)$$

where the coefficients $(\beta_{ml}, 1 \le m \le d-1, 1 \le l \le N_s)$ are the parameters which are learned in the model.

Some more notations are introduced below:

$$p_k = \frac{1}{1 + g_k} \qquad (.0.15)$$

$$g_k = e^{\tilde{g}_k} \qquad (.0.16)$$

$$v_k = e^{\tilde{v}_k} \qquad (.0.17)$$

A calculation that comes back often is $\frac{1}{p_k(1-p_k)} = \frac{(1+g_k)^2}{g_k}$.

## Bayesian model

The disease course mapping framework adds a probabilistic model over the geometrical model (which is not structurally identifiable as there can be several triplets $(t_0, p_0, v_0)$ describing the same set of trajectories). Some parameters are introduced to have priors over individual parameters.

## Priors

Here is the complete prior list:

[Population parameters (fixed effects)]

$$\tilde{g}_k \sim \mathcal{N}(\bar{g}_k, \sigma_g^2) \qquad (.0.18)$$

$$\tilde{v}_k \sim \mathcal{N}(\bar{v}_k, \sigma_v^2) \qquad (.0.19)$$

$$\beta_{ml} \sim \mathcal{N}(\bar{\beta}_{ml}, \sigma_{\beta}^2) \qquad (.0.20)$$

[Individual parameters (random effects)]

$$\tau_i \sim \mathcal{N}(\bar{\tau}, \sigma_\tau^2) \qquad (.0.21)$$

$$\xi_i \sim \mathcal{N}(0, \sigma_\xi^2) \qquad (.0.22)$$

$$s_{il} \sim \mathcal{N}(0, 1) \qquad (.0.23)$$
## Parameter classes

Let's divide the parameters into categories:

* hyperparameters: these parameters are fixed by the user ahead of estimation. These include $\sigma_g, \sigma_v, \sigma_\beta$ which are all fixed to 0.01 by default. They appear in priors involving the same scheme: $x \sim N(\bar{x}, \sigma_x^2)$. These priors are not really constraining. Intuitively they only provide the exploration strength during training for the population parameters. They do not shape the posterior of related parameters as the $\bar{g}, \bar{v}, \bar{\beta}$ are updated to maximize the likelihood at each step of the EM algorithm which amounts to setting them to $g, v, \beta$ respectively.

* model parameters: $\theta = (\bar{g}, \bar{v}, \bar{\beta}, \bar{\tau}, \sigma_\tau, \sigma_\xi)$ these are parameters which appear as the mean of gaussian priors. Intuitively they are the mean value of the related parameters, which for population parameters $(g, v, \beta)$ is the best value to use in the model formula (0.12). For individual parameters $(\tau, \xi)$ they imply that the distribution of such random effects is expected to be gaussian and $(\bar{\tau}, 0)$ is the mean of this distribution. These parameters are all maximized during the M-step of the EM algorithm used for estimation.

* latent parameters: $z = (z_{pop}, (z_i)_{1 \le i \le N})$ which are sampled during the E-step of the EM. We distinguish between population parameters $z_{pop} = (\tilde{g}, \tilde{v}, \beta)$ which involve the whole dataset (all the population) in order to be updated during estimation, and individual parameters $z_i = (\tau_i, \xi_i, s_i)$ which only require data points of individual $i$ to be updated.

## From model to data

Depending on the type of data of the observations, the link between the model $\eta$ and the data values $y_{ijk}$ is different. Here we provide the expressions for continuous data modelled with Gaussian noise, binary data modelled with Bernoulli noise and ordinal data modelled with a cumulative probit model.

**Gaussian noise** The data points are continuous values, supposedly normalized between 0 and 1, and are supposed to follow a normal law:

$$y_{ijk} \sim N(\eta_k(\psi_i(t_{ij}); \mathbf{w}_i), \sigma_k^2) \qquad (0.24)$$

where $\sigma_k$ is the standard deviation of the noise of item $k$. This parameter is learned and is a model parameter, which means it is updated during the M-step of the EM.

**Bernoulli noise** The data points are binary values, either 0 or 1, and we model them with a Bernoulli:

$$y_{ijk} \sim B(\eta_k(\psi_i(t_{ij}); \mathbf{w}_i)) \qquad (0.25)$$

There is no supplementary parameter to learn for this type of data.
**Ordinal noise** Ordinal data is a bit more complicated: observations $y_{ijk}$ are supposed to be integers between 0 and a given maximum level $L_k$. The model $\eta$ describes the cumulative distribution function instead of the probability distribution function. Due to the logistic curve going up to 1 as time passes, this implies that the ordinal scale should be ordered such that individuals start at 0 and eventually reach the maximum score in the infinite time horizon. Here is how it formulates:

$$ \forall k \in [|1, d|], \forall l \in [|1, L_k|], \mathbb{P}(y_{ijk} \ge l) = \eta_k(\psi_i^{k,l}(t_{ij}); \mathbf{w}_i) \quad (.0.26) $$

$$ \psi_i^{k,l}(t) = e^{\xi_i} (t - \tau_i) - \sum_{m=1}^{l} \delta_k^m \quad (.0.27) $$

with $\delta_k^m > 0$ being the average time delay between levels $m-1$ and $m$ for item $k$. These $(\delta_k^m)_{1 \le k \le d, 1 \le m \le L_k}$ parameters are new population parameters, thus latent, i.e. sampled during the E-step of the algorithm. This adds the following prior:

$$ \delta_k^m \sim \mathcal{N}(\delta_k^{\bar{m}}, \sigma_{\delta}^2) \quad (.0.28) $$

with $(\delta_k^{\bar{m}})_{1 \le k \le d, 1 \le m \le L_k}$ being the corresponding new model parameters and $\sigma_{\delta}$ is a new hyperparameter fixed at 0.01.

Complete equation (0.26) with logistic curves model formulates as:

$$ \eta_k(\psi_i^{k,l}(t_{ij}); \mathbf{w}_i) = \left( 1 + \left( \frac{1}{p_k} - 1 \right) \exp \left( - \frac{v_k \left( e^{\xi_i} (t_{ij} - \tau_i) - \sum_{m=1}^{l} \delta_k^m \right) + w_{ik}}{p_k (1 - p_k)} \right) \right)^{-1} $$

## Likelihood

Estimating the model is done by maximizing the likelihood. Model parameters $\theta$ are not random, which is the meaning of the ";" separator before including $\theta$ in the next formulas. We usually compute the log-likelihood formulas. The total log-likelihood of the model writes:

$$ \log q(\mathbf{y}, \mathbf{z}; \theta) = \log q(\mathbf{y} | \mathbf{z}; \theta) \quad (.0.29) $$

$$ + \log q(\mathbf{z}_{pop}; \theta) \quad (.0.30) $$

$$ + \sum_{i=1}^{N} \log q(\mathbf{z}_i | \mathbf{z}_{pop}; \theta) \quad (.0.31) $$

We will detail each term below.

## Data attachment

The first term, referred to as *attachment*, is measuring how the model fits the data with the parameters $\theta$ and $\mathbf{z}$. This term depends on the current noise model, thus on the data type:
$$ \log(q(\mathbf{y}|\mathbf{z}; \theta)) = - \sum_{k=1}^{d} \left( N_{tot} \log(\sigma_k \sqrt{2\pi}) + \frac{1}{2\sigma_k^2} \sum_{i=1}^{N} \sum_{j=1}^{N_i} (y_{i,j,k} - \eta_k(\psi_i(t_{ij}); \mathbf{w}_i))^2 \right) \quad [\text{Gaussian noise}] $$

$$ \log(q(\mathbf{y}|\mathbf{z}; \theta)) = \sum_{i=1}^{N} \sum_{j=1}^{N_i} \sum_{k=1}^{d} y_{i,j,k} \log(\eta_k(\psi_i(t_{ij}); \mathbf{w}_i)) + (1 - y_{i,j,k}) \log(1 - \eta_k(\psi_i(t_{ij}); \mathbf{w}_i)) \quad [\text{Bernoulli noise}] $$

$$ \log(q(\mathbf{y}|\mathbf{z}; \theta)) = \sum_{i=1}^{N} \sum_{j=1}^{N_i} \sum_{k=1}^{d} \sum_{l=1}^{L_k} \mathbb{1}(y_{i,j,k} = l) \log(\eta_k(\psi_i^{k,l}(t_{ij}); \mathbf{w}_i) - \eta_k(\psi_i^{k,l+1}(t_{ij}); \mathbf{w}_i)) \quad [\text{Ordinal noise}] $$

with $\eta_k(\psi_i^{k,0}(t_{ij}), \mathbf{w}_i) = 1$ and $\eta_k(\psi_i^{k,L_k+1}(t_{ij}), \mathbf{w}_i) = 0$ in the ordinal case (i.e. $\mathbb{P}(y_{ijk} \ge 0) = 1$ and $\mathbb{P}(y_{ijk} \ge L_k) = 0$), since they were not defined by equation (0.26). For the Gaussian noise equation, $N_{tot} = \sum_{i=1}^{N} N_i$ is the total number of observations. Data could be missing, with different features being observed or not at various timepoints so the formula above needs to be adjusted by masking out missing timepoints (in this case the $N_{tot}$ would vary depending on $k$ and would be the total of observations in feature $k$). Notice that we obtain mean squared error (MSE) in the Gaussian noise and a binary crossentropy (BCE) in the Bernoulli noise, which in Machine Learning are the standard losses to deal with continuous and binary data respectively. The next two prior log-likelihood terms correspond to the regularity term:

## Population parameters prior

The second term is the log-likelihood of the population parameters' prior. This includes the Gaussian prior for $(\tilde{g}, \tilde{v}, \tilde{\beta})$:

$$ \log q(\mathbf{z}_{pop}; \theta) = -d \log(\sigma_g \sqrt{2\pi}) - \frac{1}{2\sigma_g^2} \|\bar{\mathbf{g}} - \tilde{\mathbf{g}}\|^2 \quad (.0.32) $$

$$ -d \log(\sigma_v \sqrt{2\pi}) - \frac{1}{2\sigma_v^2} \|\bar{\mathbf{v}} - \tilde{\mathbf{v}}\|^2 \quad (.0.33) $$

$$ -(d-1)N_s \log(\sigma_\beta \sqrt{2\pi}) - \frac{1}{2\sigma_\beta^2} \sum_{m=1}^{d-1} \sum_{l=1}^{N_s} (\bar{\beta}_{ml} - \tilde{\beta}_{ml})^2 \quad (.0.34) $$

## Individual parameters prior

The last term is the log-likelihood of the individual parameters' prior. We have thus:

$$ \log q(\mathbf{z}_i | \mathbf{z}_{pop}; \theta) = - \log(\sigma_\tau \sqrt{2\pi}) - \frac{1}{2\sigma_\tau^2} (\bar{\tau} - \tau_i)^2 \quad (.0.35) $$

$$ - \log(\sigma_\xi \sqrt{2\pi}) - \frac{1}{2\sigma_\xi^2} \xi_i^2 \quad (.0.36) $$

$$ -N_s \log(\sqrt{2\pi}) - \frac{1}{2} \|\mathbf{s}_i\|^2 \quad (.0.37) $$
Gradients

We provide here the formulas for the partial derivatives of the log-likelihood with respect to the dif-
ferent latent parameters. We only compute gradients for latent parameters as the other parameters
are either fixed or maximized during the M-step, thus not requiring gradient information during
optimization. On the other hand, latent parameter sampling can be improved with gradient. We
can also mention the personalization process of the DCM which involves the search of the maximum
a posteriori for individual latent parameters while all other parameters are fixed. A local optimum
can be obtained by an optimization solver, the performance of which is greatly improved with the
gradient being provided.

Noise gradient

We first derive the data attachment term of formula (.0.29). The formulas for Gaussian and Bernoulli noise allow us to chain derivatives, so it is easier to first derive the log-likelihood with respect to the model and then derive the model η with respect to the different parameters (where η is the short notation for η_k(ψ_i(t_{ij}); w_i)). Since we can rewrite log q(y | z; θ) = f(η(x)) for any parameter x, we can chain rule the gradient for any list of parameters X:

$$
\nabla_X \log q(\mathbf{y} | \mathbf{z}; \theta) = \frac{\partial \log q(\mathbf{y} | \mathbf{z}; \theta)}{\partial \eta} \nabla_X \eta \quad (.0.38)
$$

and for each different noise model we have:

$$
\frac{\partial \log q(\mathbf{y} | \mathbf{z}; \theta)}{\partial \eta} = - \sum_k \frac{1}{\sigma_k^2} \sum_{i,j} (\eta - y_{i,j,k}) \quad [\text{Gaussian noise}]
$$

$$
\frac{\partial \log q(\mathbf{y} | \mathbf{z}; \theta)}{\partial \eta} = \sum_{i,j,k} \frac{y_{i,j,k} - \eta}{\eta(1 - \eta)} \quad [\text{Bernoulli noise}]
$$

For the ordinal noise this is a little different due to the likelihood term depending on η_k(ψ_i^{k,l}(t_{ij}); w_i) and η_k(ψ_i^{k,l+1}(t_{ij}); w_i) for which the simplified notations are η_k^l and η_k^{l+1} respectively. The gradient for any list of parameters X thus becomes:

$$
\nabla_X \log q(\mathbf{y} | \mathbf{z}; \theta) = \sum_{i,j,k,l} \mathbb{1}(y_{i,j,k} = l) \frac{\nabla_X \eta_k^l - \nabla_X \eta_k^{l+1}}{\eta_k^l - \eta_k^{l+1}} \quad [\text{Ordinal noise}]
$$

Model gradients

We list here all the derivatives with respect to the different parameters. Let h denote the term
$$ \tilde{g} - \frac{v_k \psi_i(t_{ij}) + w_{ik}}{p_k(1-p_k)} $$ so that equation (.0.12) becomes $\eta = \frac{1}{1+\exp(h)}$. We can chain the gradients for
simpler calculus.

$$
\frac{\partial \eta}{\partial h} = - \exp(h) \eta^2 = \eta(\eta - 1) \quad (.0.39)
$$

and then by computing the derivatives of h w.r.t. the parameters we obtain:
[Individual parameters (random effects)]

$$ \frac{\partial \eta}{\partial \tau_i} = \eta(\eta - 1) \frac{v_k e^{\xi_i}}{p_k(1 - p_k)} \quad (.0.40) $$

$$ \frac{\partial \eta}{\partial \xi_i} = \eta(\eta - 1) \left( - \frac{v_k \psi_i(t_{ij})}{p_k(1 - p_k)} \right) \quad (.0.41) $$

$$ \frac{\partial \eta}{\partial s_{il}} = \eta(\eta - 1) \left( - \frac{\beta_{kl}}{p_k(1 - p_k)} \right) \quad (.0.42) $$

[Population parameters (fixed effects)]

$$ \frac{\partial \eta}{\partial \beta_{ml}} = \eta(\eta - 1) \left( - \frac{s_{il} \mathcal{B}_{mk}}{p_k(1 - p_k)} \right) \quad (.0.43) $$

$$ \frac{\partial \eta_l}{\partial \delta_k^m} = \eta(\eta - 1) \frac{1(l > m)v_k}{p_k(1 - p_k)} \quad [\text{Ordinal added parameters}] \quad (.0.44) $$

For $v$ and $g$ computations are trickier. Indeed, the space-shift vector $w_i = \sum_{l=1}^{N_s} \sum_{m=1}^{d-1} \beta_{ml} s_{il} \mathcal{B}_m$ is a complicated function of $v$ and $p$ since the orthonormal basis $\mathcal{B}$ of $\text{Span}(v)^\perp$ is obtained with a Householder method. The basis is orthonormal with the inner product of the Riemannian manifold, which depends on the metric $G$ of the current point $p$: $\langle x, y \rangle_p = \langle x, G(p)y \rangle$ where $\langle, \rangle$ is the Euclidean inner product. In the logistic curves case, the metric is diagonal due the product of manifolds and the diagonal elements are: $G(p)_{ii} = \frac{1}{p_i^2(1-p_i)^2} = \frac{(1+\exp(\tilde{g}_i))^4}{\exp(\tilde{g}_i)^2}$. The Householder method arbitrarily selects a dimension, which here is selected as the first one without loss of generality, and then computes the orthonormal basis according to the Euclidean inner product with the following steps:

$$ \omega = G(p)v \quad (.0.45) $$

$$ \kappa = \operatorname{sign}(\omega_1) ||\omega|| \quad (.0.46) $$

$$ u = \omega - \kappa e_1 \quad (.0.47) $$

$$ \tilde{u} = \frac{u}{||u||} \quad (.0.48) $$

$$ B = I_d - 2\tilde{u}\tilde{u}^T \quad (.0.49) $$

$$ \mathcal{B} = B_{2:d} \quad (.0.50) $$

where $e$ is the canonical basis, $I_d$ the identity matrix and $\mathcal{B}$ is the matrix $B$ stripped of the first column. This first column is colinear to $G(p)v$, so the other columns are orthogonal to it in the canonical inner product, thus orthogonal to $v$ in the manifold inner product.

Differentiating through this algorithm yields the following gradients:
$$ \frac{\partial \eta}{\partial \tilde{g}_{k'}} = \eta(\eta - 1) \left( 1 - \delta_{kk'} (v_k \psi_i(t_{ij}) + w_{ik})(e^{\tilde{g}_k} - e^{-\tilde{g}_k}) - \frac{1}{p_k(1-p_k)} \sum_{l=1}^{N_s} \sum_{m=1}^{d-1} \beta_{ml} s_{il} \frac{\partial \mathcal{B}_{mk}}{\partial \tilde{g}_{k'}} \right) \quad (.0.51) $$

$$ \frac{\partial \eta}{\partial \tilde{v}_{k'}} = \eta(\eta - 1) \frac{1}{p_k(1-p_k)} \left( \delta_{kk'} v_k \psi_i(t_{ij}) + \sum_{l=1}^{N_s} \sum_{m=1}^{d-1} \beta_{ml} s_{il} \frac{\partial \mathcal{B}_{mk}}{\partial \tilde{v}_{k'}} \right) \quad (.0.52) $$

$$ \frac{\partial \mathcal{B}_{mk}}{\partial \tilde{g}_{k'}} = 2(e^{\tilde{g}_{k'}} + 2 + e^{-\tilde{g}_{k'}})(e^{\tilde{g}_{k'}} - e^{-\tilde{g}_{k'}}) \frac{\partial \mathcal{B}_{mk}}{\partial \tilde{\omega}_{k'}} \quad (.0.53) $$

$$ \frac{\partial \mathcal{B}_{mk}}{\partial \tilde{v}_{k'}} = \frac{1}{p_k^2(1-p_k)^2} \frac{\partial \mathcal{B}_{mk}}{\partial \tilde{\omega}_{k'}} \quad (.0.54) $$

$$ \frac{\partial \mathcal{B}_{mk}}{\partial \tilde{\omega}_{k'}} = -2(\tilde{u}_m \frac{\partial \tilde{u}_k}{\partial \tilde{\omega}_{k'}} + \tilde{u}_k \frac{\partial \tilde{u}_m}{\partial \tilde{\omega}_{k'}}) \quad (.0.55) $$

$$ \frac{\partial \tilde{u}_k}{\partial \omega_{k'}} = \frac{1}{\|\mathbf{u}\|} (\delta_{kk'} - \delta_{k1}\mathrm{sign}(\omega_1) \frac{\omega_k'}{\|\omega\|}) - \frac{u_k}{\|\mathbf{u}\|^3} (u_{k'} + \omega_{k'}(1 - \frac{|\omega_1|}{\|\omega\|})) \quad (.0.56) $$

where $\delta$ is the classical Kronecker symbol. Notice that even though $\eta$ stands for the model value of individual $i$ for dimension $k$, its derivative w.r.t to population parameters $g_{k'}$ and $v_{k'}$ with $k' \neq k$ is not null, whereas for individual parameters such as $\tau_{i'}$, $\xi_{i'}$ the gradient is null if $i' \neq i$.

## Priors gradient

For each of the individual and population parameters we also have a gradient coming from the prior (the regularity term). As we only have Gaussian priors, for any parameter $x$:

$$ \frac{\partial \log q(x; \theta)}{\partial x} = - \frac{x - \bar{x}}{\sigma_x^2} \quad (.0.57) $$

where $\bar{x}$ and $\sigma_x$ are the mean and standard deviation of the prior of $x$.

## Fisher information

The Fisher information is the variance of the partial derivative of the log-likelihood w.r.t. a given parameter $X$. In a multidimensional setting we have the Fisher information matrix: $\mathcal{I}(X) = \mathbb{E}(\nabla_X \log q(y; X) \nabla_X \log q(y; X)^T | X)$. The formula requires a distribution for the observations, however in the DCM model we provided the noise model as a distribution for $y$ given the time $t$. We never assumed a distribution for the times $t_{ij}$. We do not want to specify this distribution as it varies a lot depending on the disease and dataset and integrating over such a distribution is intractable. However we can specify the information given the time of the observation $t$. Then we compute the empirical Fisher information by summing over the observed times $t_{ij}$. We provide below the formulas for all the latent parameters $z = (z_{pop}, z_i)$ of the model:
$$ \mathcal{I}(t, z) = \sum_{k=1}^{d} \left( \frac{1}{\sigma_k^2} \nabla_z \eta \nabla_z \eta^T \right) + R \quad [\text{Gaussian noise}] (.0.58) $$

$$ \mathcal{I}(t, z) = \sum_{k=1}^{d} \left( \frac{1}{\eta(1-\eta)} \nabla_z \eta \nabla_z \eta^T \right) + R \quad [\text{Bernoulli noise}] (.0.59) $$

$$ \mathcal{I}(t, z) = \sum_{k=1}^{d} \sum_{l=1}^{L_k} \left( \frac{(\nabla_z \eta_l - \nabla_z \eta_{l+1})(\nabla_z \eta_l - \nabla_z \eta_{l+1})^T}{\eta_l - \eta_{l+1}} \right) + R \quad (.0.60) $$

[Ordinal noise]

$$ \mathcal{I}(z) = \sum_{i=1}^{N} \sum_{j=1}^{N_i} \mathcal{I}(t_{ij}, z) \quad [\text{Empirical Fisher information}] (.0.61) $$

where R is the information matrix of the priors of z, i.e. $R_{a,b} = \frac{(a-\bar{a})(b-\bar{b})}{\sigma_a^2 \sigma_b^2}$ for the coefficients corresponding to two latent variables a, b with Gaussian priors.

In the ordinal case, for any individual parameter x, we have $\frac{\partial \eta_l}{\partial x} = \eta_l(\eta_l - 1)\frac{\partial h_l}{\partial x} = \eta_l(\eta_l - 1)\frac{\partial h}{\partial x}$ which does not depend on the level l anymore, thus simplifying the formula into: $\mathcal{I}(t, z_i) = \sum_{k=1}^{d} \sum_{l=1}^{L_k} \frac{(\eta_l(\eta_{l-1}-\eta_{l+1})(\eta_{l+1}-\eta_{l-1}))^2}{\eta_l-\eta_{l+1}} \nabla_{z_i} h \nabla_{z_i} h^T + R$.

The Fisher information formulas for model parameters $\theta$ have not been explicitly computed here as we were mostly interested in latent parameters (cf sections 2.4.1). Note that these model parameters do not appear in the model expression $\eta$. In practice this means that the Fisher information for model parameters is not varying with the times of observations $t_{ij}$, so we cannot perform the same analysis as in section 2.4.1. As for the use of the Fisher information matrix for the natural gradient in Riemannian samplers<fnref n="35" />, this is only relevant for sampled parameters, i.e. latent variables. However the Fisher information can still be used to approximate the variance of the estimation of the model parameters, using the fact that $\sqrt{n}(\hat{\theta} - \theta^*)$ converges in distribution towards $N(0, \mathcal{I}(\theta^*)^{-1})$ as the number of observations n grows to infinity, where $\hat{\theta}$ is the estimator for $\theta$ and $\theta^*$ is the true parameter.

<footnote n="35">³⁵M. Girolami and B. Calderhead, “Riemann manifold Langevin and Hamiltonian Monte Carlo methods,” en, *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, vol. 73, no. 2, pp. 123–214, 2011.</footnote>
<figure><img src="image_57.png" /></figure>
Bibliography

[1] C. Abi Nader, N. Ayache, P. Robert, M. Lorenzi, and Alzheimer's Disease Neuroimaging Initiative, "Monotonic Gaussian Process for spatio-temporal disease progression modeling in brain imaging data," eng, *NeuroImage*, vol. 205, p. 116266, Jan. 2020.

[2] J. Abrevaya, Y.-C. Hsu, and R. P. Lieli, "Estimating Conditional Average Treatment Effects," *Journal of Business & Economic Statistics*, vol. 33, no. 4, pp. 485-505, Oct. 2015.

[3] S. Allassonnière and J. Chevallier, "A new class of stochastic EM algorithms. Escaping local maxima and handling intractable sampling," en, *Computational Statistics & Data Analysis*, vol. 159, p. 107159, Jul. 2021.

[4] S. Allassonnière and E. Kuhn, "Stochastic algorithm for Bayesian mixture effect template estimation," en, *ESAIM: Probability and Statistics*, vol. 14, pp. 382-408, Jan. 2010.

[5] S. Allassonnière and E. Kuhn, "Convergent stochastic Expectation Maximization algorithm with efficient sampling in high dimension. Application to deformable template model estimation," en, *Computational Statistics & Data Analysis*, vol. 91, pp. 4-19, Nov. 2015.

[6] L. Arrington, S. Ueckert, M. Ahamadi, S. Macha, and M. O. Karlsson, "Performance of longitudinal item response theory models in shortened or partial assessments," en, *Journal of Pharmacokinetics and Pharmacodynamics*, vol. 47, no. 5, pp. 461-471, Oct. 2020.

[7] D. Arthur and S. Vassilvitskii, "K-means++: The Advantages of Careful Seeding," en,

[8] P. C. Austin and A. Laupacis, "A tutorial on methods to estimating clinically and policy-meaningful measures of treatment effects in prospective observational studies: A review," eng, *The International Journal of Biostatistics*, vol. 7, no. 1, p. 6, 2011.

[9] L. A. Beckett, "Maximum Likelihood Estimation in Mallows's Model Using Partially Ranked Data," en, in *Probability Models and Statistical Analyses for Ranking Data*, M. A. Fligner and J. S. Verducci, Eds., ser. Lecture Notes in Statistics, New York, NY: Springer, 1993, pp. 92-107.

[10] A. L. Benabid, "Deep brain stimulation for Parkinson's disease," en, *Current Opinion in Neurobiology*, vol. 13, no. 6, pp. 696-706, Dec. 2003.

[11] D. Bertolini, K. Arnemann, D. Hall, and J. Walsh, "Machine Learning Enables Smaller ALS Clinical Trials (P1-13.003)," en, *Neurology*, vol. 98, no. 18 Supplement, May 2022.

[12] M. Betancourt, *A Conceptual Introduction to Hamiltonian Monte Carlo*, Jul. 2018.

[13] R. M. A. de Bie, C. E. Clarke, A. J. Espay, S. H. Fox, and A. E. Lang, "Initiation of pharmacological therapy in Parkinson's disease: When, why, and how," eng, *The Lancet. Neurology*, vol. 19, no. 5, pp. 452-461, May 2020.
[14] M. Bilgel, J. L. Prince, D. F. Wong, S. M. Resnick, and B. M. Jedynak, "A multivariate non-linear mixed effects model for longitudinal image analysis: Application to amyloid imaging," en, *NeuroImage*, vol. 134, pp. 658-670, Jul. 2016.

[15] J. F. Bobb, L. Valeri, B. Claus Henn, D. C. Christiani, R. O. Wright, M. Mazumdar, J. J. Godleski, and B. A. Coull, "Bayesian kernel machine regression for estimating the health effects of multi-pollutant mixtures," eng, *Biostatistics (Oxford, England)*, vol. 16, no. 3, pp. 493-508, Jul. 2015.

[16] A. Bône, "Learning adapted coordinate systems for the statistical analysis of anatomical shapes. Applications to Alzheimer's disease progression modeling," en, Ph.D. dissertation, Sorbonne Université, Jan. 2020.

[17] A. Bône, M. Louis, B. Martin, and S. Durrleman, "Deformetrica 4: An Open-Source Software for Statistical Shape Analysis," Sep. 2018.

[18] G. E. P. Box and D. R. Cox, "An Analysis of Transformations," en,

[19] H. Braak, K. D. Tredici, U. Rüb, R. A. I. de Vos, E. N. H. Jansen Steur, and E. Braak, "Staging of brain pathology related to sporadic Parkinson's disease," en, *Neurobiology of Aging*, vol. 24, no. 2, pp. 197-211, Mar. 2003.

[20] P. Chan and N. Holford, "Drug Treatment Effects on Disease Progression," *Annual Review of Pharmacology and Toxicology*, vol. 41, no. 1, pp. 625-659, 2001.

[21] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud, *Neural Ordinary Differential Equations*, Dec. 2019.

[22] S. Chib and E. Greenberg, "Understanding the Metropolis-Hastings Algorithm," *The American Statistician*, vol. 49, no. 4, pp. 327-335, Nov. 1995.

[23] S. Chib and B. H. Hamilton, "Semiparametric Bayes analysis of longitudinal data treatment models," en, *Journal of Econometrics*, vol. 110, no. 1, pp. 67-89, Sep. 2002.

[24] H. A. Chipman, E. I. George, and R. E. McCulloch, "BART: Bayesian additive regression trees," *The Annals of Applied Statistics*, vol. 4, no. 1, pp. 266-298, Mar. 2010.

[25] R. Cilia, E. Cereda, A. Akpalu, F. S. Sarfo, M. Cham, R. Laryea, V. Obese, K. Oppon, F. Del Sorbo, S. Bonvegna, A. L. Zecchinelli, and G. Pezzoli, "Natural history of motor symptoms in Parkinson's disease and the long-duration response to levodopa," *Brain*, vol. 143, no. 8, pp. 2490-2501, Aug. 2020.

[26] O. o. t. Commissioner, *FDA Converts Novel Alzheimer's Disease Treatment to Traditional Approval*, en, Jul. 2023.

[27] S. F. Cook and R. R. Bies, "Disease Progression Modeling: Key Concepts and Recent Developments," en, *Current Pharmacology Reports*, vol. 2, no. 5, pp. 221-230, Oct. 2016.

[28] R. Couronné, "Progression models for Parkinson's Disease," en, Ph.D. dissertation, Sorbonne Université, Sep. 2021.

[29] R. Couronné, A. Valladier, M. Vidhaillet, J. C. Corvol, S. Lehéricy, and S. Durrleman, "Modeling the progression of Parkinson's Disease: Comparison of subjects with and without Sleep Disorders," en, p. 4,

[30] R. Couronné, M. Vidailhet, J.-C. Corvol, S. Lehéricy, and S. Durrleman, "Learning disease progression models with longitudinal data and missing values," in *ISBI 2019 - International Symposium on Biomedical Imaging*, Venice, Italy, Apr. 2019.
[31] A. Cremaschi, M. De Iorio, Y. Seng Chong, B. Broekman, M. J. Meaney, and M. Z. L. Kee, “A Bayesian nonparametric approach to dynamic item-response modeling: An application to the GUSTO cohort study,” en, *Statistics in Medicine*, vol. 40, no. 27, pp. 6021–6037, 2021.

[32] T. Dao and C. D. Sa, “Gaussian Quadrature for Kernel Features,” en,

[33] V. Debavelaere, S. Durrleman, S. Allassonnière, and for the Alzheimer’s Disease Neuroimaging Initiative, “Learning the Clustering of Longitudinal Shape Data Sets into a Mixture of Independent or Branching Trajectories,” en, *International Journal of Computer Vision*, vol. 128, no. 12, pp. 2794–2809, Dec. 2020.

[34] M. Delattre, M. Lavielle, and M.-A. Poursat, “A note on BIC in mixed-effects models,” en, *Electronic Journal of Statistics*, vol. 8, no. 1, Jan. 2014.

[35] B. Delyon, M. Lavielle, and E. Moulines, “Convergence of a Stochastic Approximation Version of the EM Algorithm,” *The Annals of Statistics*, vol. 27, no. 1, pp. 94–128, 1999.

[36] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum Likelihood from Incomplete Data Via the EM Algorithm,” en, *Journal of the Royal Statistical Society: Series B (Methodological)*, vol. 39, no. 1, pp. 1–22, 1977.

[37] T. G. Dietterich, “Ensemble Methods in Machine Learning,” en, in *Multiple Classifier Systems*, ser. Lecture Notes in Computer Science, Berlin, Heidelberg: Springer, 2000, pp. 1–15.

[38] M. Do Carmo, *Riemannian Geometry*, en.

[39] M. C. Donohue, H. Jacqmin-Gadda, M. Le Goff, R. G. Thomas, R. Raman, A. C. Gamst, L. A. Beckett, C. R. Jack, M. W. Weiner, J.-F. Dartigues, and P. S. Aisen, “Estimating long-term multivariate progression from short-term data,” *Alzheimer's & dementia : the journal of the Alzheimer's Association*, vol. 10, no. 0, S400–S410, Oct. 2014.

[40] J. Dumurgier and C. Tzourio, “Epidemiology of neurological diseases in older adults,” en, *Revue Neurologique*, vol. 176, no. 9, pp. 642–648, Nov. 2020.

[41] S. Durrleman, X. Pennec, A. Trouvé, J. Braga, G. Gerig, and N. Ayache, “Toward a Comprehensive Framework for the Spatiotemporal Statistical Analysis of Longitudinal Shape Data,” en, *International Journal of Computer Vision*, vol. 103, no. 1, pp. 22–59, May 2013.

[42] R. S. Eisinger, C. W. Hess, D. Martinez-Ramirez, L. Almeida, K. D. Foote, M. S. Okun, and A. Gunduz, “Motor subtype changes in early Parkinson’s disease,” en, *Parkinsonism & Related Disorders*, vol. 43, pp. 67–72, Oct. 2017.

[43] A. Elbaz, L. Carcaillon, S. Kab, and F. Moisan, “Epidemiology of Parkinson’s disease,” en, *Revue Neurologique*, Neuroepidemiology, vol. 172, no. 1, pp. 14–26, Jan. 2016.

[44] S. E. Embretson and S. P. Reise, *Item response theory for psychologists* (Item response theory for psychologists). Mahwah, NJ, US: Lawrence Erlbaum Associates Publishers, 2000.

[45] A. Espay, A. Fasano, B. van Nuenen, M. Payne, A. Snijders, and B. Bloem, ““On” state freezing of gait in Parkinson disease,” *Neurology*, vol. 78, no. 7, pp. 454–457, Feb. 2012.

[46] D. Ferreira, A. Nordberg, and E. Westman, “Biological subtypes of Alzheimer disease: A systematic review and meta-analysis,” en, *Neurology*, vol. 94, no. 10, pp. 436–448, Mar. 2020.

[47] C. R. Fields, N. Bengoa-Vergniory, and R. Wade-Martins, “Targeting Alpha-Synuclein as a Therapy for Parkinson’s Disease,” *Frontiers in Molecular Neuroscience*, vol. 12, 2019.
[48] B. Fischl, D. H. Salat, E. Busa, M. Albert, M. Dieterich, C. Haselgrove, A. van der Kouwe, R. Killiany, D. Kennedy, S. Klaveness, A. Montillo, N. Makris, B. Rosen, and A. M. Dale, "Whole brain segmentation: Automated labeling of neuroanatomical structures in the human brain," eng, *Neuron*, vol. 33, no. 3, pp. 341-355, Jan. 2002.

[49] R. A. Fisher, "XV.—The Correlation between Relatives on the Supposition of Mendelian Inheritance.," en, *Earth and Environmental Science Transactions of The Royal Society of Edinburgh*, vol. 52, no. 2, pp. 399-433, Jan. 1919.

[50] T. Fletcher, "Geodesic Regression on Riemannian Manifolds," en, *Sep.* 2011, p. 75.

[51] H. M. Fonteijn, M. Modat, M. J. Clarkson, J. Barnes, M. Lehmann, N. Z. Hobbs, R. I. Scahill, S. J. Tabrizi, S. Ourselin, N. C. Fox, and D. C. Alexander, "An event-based model for disease progression and its application in familial Alzheimer's disease and Huntington's disease," eng, *NeuroImage*, vol. 60, no. 3, pp. 1880-1889, Apr. 2012.

[52] N. Fournier and S. Durrleman, "Covariate-Aware Longitudinal Modelling for Neurodegenerative Diseases," en,

[53] A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin, "Bayesian Data Analysis," en, *Chapman and Hall/CRC, Nov. 2013*.

[54] A. Gelman and D. B. Rubin, "Inference from Iterative Simulation Using Multiple Sequences," *Statistical Science*, vol. 7, no. 4, pp. 457-472, Nov. 1992.

[55] M. Gilat, N. D'Cruz, P. Ginis, W. Vandenberghe, and A. Nieuwboer, "Freezing of gait and levodopa," *English, The Lancet Neurology*, vol. 20, no. 7, pp. 505-506, Jul. 2021.

[56] M. Girolami and B. Calderhead, "Riemann manifold Langevin and Hamiltonian Monte Carlo methods," en, *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, vol. 73, no. 2, pp. 123-214, 2011.

[57] C. G. Goetz, B. C. Tilley, S. R. Shaftman, G. T. Stebbins, S. Fahn, P. Martinez-Martin, W. Poewe, C. Sampaio, M. B. Stern, R. Dodel, B. Dubois, R. Holloway, J. Jankovic, J. Kulisevsky, A. E. Lang, A. Lees, S. Leurgans, P. A. LeWitt, D. Nyenhuis, C. W. Olanow, O. Rascol, A. Schrag, J. A. Teresi, J. J. van Hilten, and N. LaPelle, "Movement Disorder Society-sponsored revision of the Unified Parkinson's Disease Rating Scale (MDS-UPDRS): Scale presentation and clinimetric testing results: MDS-UPDRS: Clinimetric Assessment," en, *Movement Disorders*, vol. 23, no. 15, pp. 2129-2170, Nov. 2008.

[58] R. Gorter, J.-P. Fox, and J. W. R. Twisk, "Why item response theory should be used for longitudinal questionnaire data analysis in medical research," *BMC Medical Research Methodology*, vol. 15, no. 1, p. 55, Jul. 2015.

[59] G. Gottipati, A. Berges, S. Yang, C. Chen, M. Karlsson, and E. Plan, "Item Response Model Adaptation for Analyzing Data from Different Versions of Parkinson's Disease Rating Scales," *Pharmaceutical Research*, vol. 36, Jul. 2019.

[60] D. P. Green and H. L. Kern, "Modeling Heterogeneous Treatment Effects in Survey Experiments with Bayesian Additive Regression Trees," *The Public Opinion Quarterly*, vol. 76, no. 3, pp. 491-511, 2012.

[61] S. Gruffaz, P.-E. Poulet, E. Maheux, B. Jedynak, and S. DURRLEMAN, "Learning Riemannian metric for disease progression modeling," in *Advances in Neural Information Processing Systems*, M. Ranzato, A. Beygelzimer, Y. Dauphin, P. S. Liang, and J. W. Vaughan, Eds., vol. 34, Curran Associates, Inc., 2021, pp. 23780-23792.
[62] H. Haario, E. Saksman, and J. Tamminen, "An Adaptive Metropolis Algorithm," *Bernoulli*, vol. 7, no. 2, pp. 223-242, 2001.

[63] O. Hansson, J. Seibyl, E. Stomrud, H. Zetterberg, J. Q. Trojanowski, T. Bittner, V. Lifke, V. Corradini, U. Eichenlaub, R. Batrla, K. Buck, K. Zink, C. Rabe, K. Blennow, L. M. Shaw, Swedish BioFINDER study group, and Alzheimer's Disease Neuroimaging Initiative, "CSF biomarkers of Alzheimer's disease concord with amyloid-PET and predict clinical progression: A study of fully automated immunoassays in BioFINDER and ADNI cohorts," eng, *Alzheimer's & Dementia: The Journal of the Alzheimer's Association*, vol. 14, no. 11, pp. 1470-1481, Nov. 2018.

[64] J. A. Hardy and G. A. Higgins, "Alzheimer's disease: The amyloid cascade hypothesis," eng, *Science* (*New York, N.Y.*), vol. 256, no. 5054, pp. 184-185, Apr. 1992.

[65] R. Harpaz and R. Haralick, "TR-2006001: The EM Algorithm As a Lower Bound Optimization Technique," *Computer Science Technical Reports*, Jan. 2006.

[66] W. K. Hastings, "Monte Carlo sampling methods using Markov chains and their applications," *Biometrika*, vol. 57, no. 1, pp. 97-109, Apr. 1970.

[67] S. r. Hauberg, O. Freifeld, and M. Black, "A Geometric take on Metric Learning," in *Advances in Neural Information Processing Systems*, vol. 25, Curran Associates, Inc., 2012.

[68] T. R. T. Have, A. R. Kunselman, E. P. Pulkstenis, and J. R. Landis, "Mixed Effects Logistic Regression Models for Longitudinal Binary Response Data with Informative Drop-Out," *Biometrics*, vol. 54, no. 1, pp. 367-383, 1998.

[69] B. He and S. Luo, "Joint modeling of multivariate longitudinal measurements and survival data with applications to Parkinson's disease," eng, *Statistical Methods in Medical Research*, vol. 25, no. 4, pp. 1346-1358, Aug. 2016.

[70] M. M. Hoehn and M. D. Yahr, "Parkinsonism: Onset, progression, and mortality," en, *Neurology*, vol. 17, no. 5, pp. 427-427, May 1967.

[71] T. Hofmann, B. Schölkopf, and A. J. Smola, "Kernel methods in machine learning," *The Annals of Statistics*, vol. 36, no. 3, pp. 1171-1220, Jun. 2008.

[72] M. D. Homan and A. Gelman, "The No-U-turn sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo," *The Journal of Machine Learning Research*, vol. 15, no. 1, pp. 1593-1623, Jan. 2014.

[73] J. Horsager, K. B. Andersen, K. Knudsen, C. Skjærbæk, T. D. Fedorova, N. Okkels, E. Schaeffer, S. K. Bonkat, J. Geday, M. Otto, M. Sommerauer, E. H. Danielsen, E. Bech, J. Kraft, O. L. Munk, S. D. Hansen, N. Pavese, R. Göder, D. J. Brooks, D. Berg, and P. Borghammer, "Brain-first versus body-first Parkinson's disease: A multimodal imaging case-control study," *Brain*, vol. 143, no. 10, pp. 3077-3088, Oct. 2020.

[74] A. Hyvärinen and E. Oja, "Independent component analysis: Algorithms and applications," en, *Neural Networks*, vol. 13, no. 4-5, pp. 411-430, Jun. 2000.

[75] J. G. Ibrahim, H. Chu, and L. M. Chen, "Basic Concepts and Methods for Joint Models of Longitudinal and Survival Data," *Journal of Clinical Oncology*, vol. 28, no. 16, pp. 2796-2801, Jun. 2010.

[76] J. C. Immekus, K. E. Snyder, and P. A. Ralston, "Multidimensional Item Response Theory for Factor Structure Assessment in Educational Psychology Research," *Frontiers in Education*, vol. 4, p. 45, 2019.
[77] Y. Iturria-Medina, F. M. Carbonell, R. C. Sotero, F. Chouinard-Decorte, and A. C. Evans, "Multifactorial causal model of brain (dis)organization and therapeutic intervention: Application to Alzheimer's disease," en, *NeuroImage*, vol. 152, pp. 60–77, May 2017.

[78] C. R. Jack, D. A. Bennett, K. Blennow, M. C. Carrillo, H. H. Feldman, G. B. Frisoni, H. Hampel, W. J. Jagust, K. A. Johnson, D. S. Knopman, R. C. Petersen, P. Scheltens, R. A. Sperling, and B. Dubois, "A/T/N: An unbiased descriptive classification scheme for Alzheimer disease biomarkers," eng, *Neurology*, vol. 87, no. 5, pp. 539–547, Aug. 2016.

[79] C. R. Jack, D. S. Knopman, W. J. Jagust, R. C. Petersen, M. W. Weiner, P. S. Aisen, L. M. Shaw, P. Vemuri, H. J. Wiste, S. D. Weigand, T. G. Lesnick, V. S. Pankratz, M. C. Donohue, and J. Q. Trojanowski, "Update on hypothetical model of Alzheimer's disease biomarkers," *Lancet neurology*, vol. 12, no. 2, pp. 207–216, Feb. 2013.

[80] C. R. Jack, D. S. Knopman, W. J. Jagust, L. M. Shaw, P. S. Aisen, M. W. Weiner, R. C. Petersen, and J. Q. Trojanowski, "Hypothetical model of dynamic biomarkers of the Alzheimer's pathological cascade," English, *The Lancet Neurology*, vol. 9, no. 1, pp. 119–128, Jan. 2010.

[81] J.-H. Jacobsen, A. Smeulders, and E. Oyallon, *I-RevNet: Deep Invertible Networks*, Feb. 2018.

[82] J. Jankovic, M. McDermott, J. Carter, S. Gauthier, C. Goetz, L. Golbe, S. Huber, W. Koller, C. Olanow, and I. Shoulson, "Variable expression of Parkinson's disease: A base-line analysis of the DATATOP cohort. The Parkinson Study Group," eng, *Neurology*, vol. 40, no. 10, pp. 1529–1534, Oct. 1990.

[83] B. M. Jedynak, A. Lang, B. Liu, E. Katz, Y. Zhang, B. T. Wyman, D. Raunig, C. P. Jedynak, B. Caffo, and J. L. Prince, "A Computational Neurodegenerative Disease Progression Score: Method and Results with the Alzheimer's Disease Neuroimaging Initiative Cohort," *NeuroImage*, vol. 63, no. 3, pp. 1478–1486, Nov. 2012.

[84] S. Joshi and M. Miller, "Landmark matching via large deformation diffeomorphisms," *IEEE Transactions on Image Processing*, vol. 9, no. 8, pp. 1357–1370, Aug. 2000.

[85] W. Jung, A. W. Mulyadi, and H.-I. Suk, "Unified Modeling of Imputation, Forecasting, and Prediction for AD Progression," en, in *Medical Image Computing and Computer Assisted Intervention – MICCAI 2019*, D. Shen, T. Liu, T. M. Peters, L. H. Staib, C. Essert, S. Zhou, P.-T. Yap, and A. Khan, Eds., ser. Lecture Notes in Computer Science, Cham: Springer International Publishing, 2019, pp. 168–176.

[86] Y. Jung and I. Lee, "Optimal design of experiments for optimization-based model calibration using Fisher information matrix," en, *Reliability Engineering & System Safety*, vol. 216, p. 107968, Dec. 2021.

[87] M. Kanagawa, P. Hennig, D. Sejdinovic, and B. K. Sriperumbudur, *Gaussian Processes and Kernel Methods: A Review on Connections and Equivalences*, Jul. 2018.

[88] C. Kartsonaki, "Survival analysis," en, *Diagnostic Histopathology*, Mini-Symposium: Medical Statistics, vol. 22, no. 7, pp. 263–270, Jul. 2016.

[89] V. Kmetzsch, E. Becker, D. Saracino, D. Rinaldi, A. Camuzat, I. Le Ber, O. Colliot, and f. t. P.-D. s. group for the, "Disease Progression Score Estimation From Multimodal Imaging and MicroRNA Data Using Supervised Variational Autoencoders," *IEEE Journal of Biomedical and Health Informatics*, vol. 26, no. 12, pp. 6024–6035, Dec. 2022.
[90] I. Koval, “Learning Multimodal Digital Models of Disease Progression from Longitudinal Data: Methods & Algorithms for the Description, Prediction and Simulation of Alzheimer’s Disease Progression,” en, Ph.D. dissertation, Institut Polytechnique de Paris, Jan. 2020.

[91] I. Koval, A. Bône, M. Louis, T. Lartigue, S. Bottani, A. Marcoux, J. Samper-González, N. Burgos, B. Charlier, A. Bertrand, S. Epelbaum, O. Colliot, S. Allassonnière, and S. Durrleman, “AD Course Map charts Alzheimer’s disease progression,” eng, *Scientific Reports*, vol. 11, no. 1, p. 8020, Apr. 2021.

[92] I. Koval, T. Dighiero-Brecht, A. J. Tobin, S. J. Tabrizi, R. I. Scahill, S. Tezenas du Montcel, S. Durrleman, and A. Durr, “Forecasting individual progression trajectories in Huntington disease enables more powered clinical trials,” en, *Scientific Reports*, vol. 12, no. 1, p. 18928, Nov. 2022.

[93] J. K. Kueper, M. Speechley, and M. Montero-Odasso, “The Alzheimer’s Disease Assessment Scale-Cognitive Subscale (ADAS-Cog): Modifications and Responsiveness in Pre-Dementia Populations. A Narrative Review,” *Journal of Alzheimer’s Disease*, vol. 63, no. 2, pp. 423–444.

[94] E. Kuhn and M. Lavielle, “Coupling a stochastic approximation version of EM with an MCMC procedure,” fr, *ESAIM: Probability and Statistics*, vol. 8, pp. 115–131, 2004.

[95] K. Lahouel, M. Wells, V. Rielly, E. Lew, D. Lovitz, and B. M. Jedynak, *Learning nonparametric ordinary differential equations from noisy data*, Feb. 2023.

[96] N. M. Laird and J. H. Ware, “Random-effects models for longitudinal data,” eng, *Biometrics*, vol. 38, no. 4, pp. 963–974, Dec. 1982.

[97] M. Lavielle and C. Mbogning, “An improved SAEM algorithm for maximum likelihood estimation in mixtures of non linear mixed effects models,” *Statistics and Computing*, vol. 24, no. 5, pp. 693–707, Sep. 2014.

[98] G. Lebanon, *Learning Riemannian Metrics*, Oct. 2012.

[99] M. J. Lindstrom and D. M. Bates, “Newton—Raphson and EM Algorithms for Linear Mixed-Effects Models for Repeated-Measures Data,” *Journal of the American Statistical Association*, vol. 83, no. 404, pp. 1014–1022, Dec. 1988.

[100] F. Liu, X. Huang, Y. Chen, and J. A. K. Suykens, “Random Features for Kernel Approximation: A Survey on Algorithms, Theory, and Beyond,” *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 44, no. 10, pp. 7128–7148, Oct. 2022.

[101] L. C. Liu and D. Hedeker, “A Mixed-Effects Regression Model for Longitudinal Multivariate Ordinal Data,” en, *Biometrics*, vol. 62, no. 1, pp. 261–268, 2006.

[102] G. Livingston, J. Huntley, A. Sommerlad, D. Ames, C. Ballard, S. Banerjee, C. Brayne, A. Burns, J. Cohen-Mansfield, C. Cooper, S. G. Costafreda, A. Dias, N. Fox, L. N. Gitlin, R. Howard, H. C. Kales, M. Kivimäki, E. B. Larson, A. Ogunniyi, V. Orgeta, K. Ritchie, K. Rockwood, E. L. Sampson, Q. Samus, L. S. Schneider, G. Selbæk, L. Teri, and N. Mukadam, “Dementia prevention, intervention, and care: 2020 report of the Lancet Commission,” English, *The Lancet*, vol. 396, no. 10248, pp. 413–446, Aug. 2020.

[103] F. M. Lord, *Applications of Item Response Theory To Practical Testing Problems*. New York: Routledge, Jul. 1980.
[104] M. Lorenzi, M. Filippone, G. B. Frisoni, D. C. Alexander, and S. Ourselin, "Probabilistic disease progression modeling to characterize diagnostic uncertainty: Application to staging and prediction in Alzheimer's disease," en, *NeuroImage*, Mapping diseased brains, vol. 190, pp. 56-68, Apr. 2019.

[105] M. Lorenzi, X. Pennec, G. B. Frisoni, and N. Ayache, "Disentangling normal aging from Alzheimer's disease in structural magnetic resonance images," en, *Neurobiology of Aging*, Novel Imaging Biomarkers for Alzheimer's Disease and Related Disorders (NIBAD), vol. 36, S42-S52, Jan. 2015.

[106] M. Louis, "Computational and statistical methods for trajectory analysis in a Riemannian geometry setting," en, Ph.D. dissertation, Sorbonne Université, Oct. 2019.

[107] C. Lungu, J. M. Cedarbaum, T. M. Dawson, E. R. Dorsey, C. Faraco, H. J. Federoff, B. Fiske, R. Fox, A. M. Goldfine, K. Kieburtz, E. A. Macklin, H. Matthews, G. Rafaloff, R. Saunders-Pullman, N. F. Schor, M. A. Schwarzschild, B.-A. Sieber, T. Simuni, D. J. Surmeier, A. Tamiz, M. H. Werner, C. B. Wright, and R. Wyse, "Seeking progress in disease modification in Parkinson disease," en, *Parkinsonism & Related Disorders*, vol. 90, pp. 134-141, Sep. 2021.

[108] S. Luo, H. Zou, G. T. Stebbins, M. A. Schwarzschild, E. A. Macklin, J. Chan, D. Oakes, T. Simuni, C. G. Goetz, and P. S. G. S.-P. Investigators, "Dissecting the Domains of Parkinson's Disease: Insights from Longitudinal Item Response Theory Modeling," en, *Movement Disorders*, vol. 37, no. 9, pp. 1904-1914, 2022.

[109] A. Ly, M. Marsman, J. Verhagen, R. P. P. P. Grasman, and E.-J. Wagenmakers, "A Tutorial on Fisher information," en, *Journal of Mathematical Psychology*, vol. 80, pp. 40-55, Oct. 2017.

[110] L. v. d. Maaten and G. Hinton, "Visualizing Data using t-SNE," *Journal of Machine Learning Research*, vol. 9, no. 86, pp. 2579-2605, 2008.

[111] E. Maheux, J. Ortholand, C. Birkenbihl, E. Thibeau-Sutre, M. Sood, D. Archetti, V. Bouteloup, I. Koval, and S. Durrleman, "Forecast Alzheimer's disease progression to better select patients for clinical trials," en, Jul. 2021.

[112] K. Marek, D. Jennings, S. Lasch, et al., "The Parkinson Progression Marker Initiative (PPMI)," en, *Progress in Neurobiology*, Biological Markers for Neurodegenerative Diseases, vol. 95, no. 4, pp. 629-635, Dec. 2011.

[113] R. V. Marinescu, A. Eshaghi, M. Lorenzi, A. L. Young, N. P. Oxtoby, S. Garbarino, S. J. Crutch, and D. C. Alexander, "DIVE: A spatiotemporal progression model of brain pathology in neurodegenerative disorders," en, *NeuroImage*, vol. 192, pp. 166-177, May 2019.

[114] R. V. Marinescu, N. P. Oxtoby, A. L. Young, E. E. Bron, A. W. Toga, M. W. Weiner, F. Barkhof, N. C. Fox, P. Golland, S. Klein, and D. C. Alexander, "TADPOLE Challenge: Accurate Alzheimer's disease prediction through crowdsourced forecasting of future data," eng, *PRedictive Intelligence in MEdicine. PRIME (Workshop)*, vol. 11843, pp. 1-10, Oct. 2019.

[115] C. Marras, K. R. Chaudhuri, N. Titova, and T. A. Mestre, "Therapy of Parkinson's Disease Subtypes," en, *Neurotherapeutics*, vol. 17, no. 4, pp. 1366-1377, Oct. 2020.
[116] P. Martinez-Martin, C. Rodriguez-Blazquez, M. Alvarez-Sanchez, T. Arakaki, A. Bergareche-Yarza, A. Chade, N. Garretto, O. Gershanik, M. M. Kurtis, J. C. Martinez-Castrillo, A. Mendoza-Rodriguez, H. P. Moore, M. Rodriguez-Violante, C. Singer, B. C. Tilley, J. Huang, G. T. Stebbins, and C. G. Goetz, "Expanded and independent validation of the Movement Disorder Society-Unified Parkinson's Disease Rating Scale (MDS-UPDRS)," en, *Journal of Neurology*, vol. 260, no. 1, pp. 228-236, Jan. 2013.

[117] J. M. te Marvelde, C. A. W. Glas, G. Van Landeghem, and J. Van Damme, "Application of Multidimensional Item Response Theory Models to Longitudinal Data," en, *Educational and Psychological Measurement*, vol. 66, no. 1, pp. 5-34, Feb. 2006.

[118] P. McCullagh, *Generalized Linear Models*, en. Routledge, Oct. 2018.

[119] I. G. McKeith, "Spectrum of Parkinson's disease, Parkinson's dementia, and Lewy body dementia," eng, *Neurologic clinics*, vol. 18, no. 4, pp. 865-902, Nov. 2000.

[120] P. J. Moore, T. J. Lyons, J. Gallacher, and f. t. A. D. N. Initiative, "Random forest prediction of Alzheimer's disease using pairwise selection from time series data," en, *PLOS ONE*, vol. 14, no. 2, e0211558, Feb. 2019.

[121] P. Moulaire, P. E. Poulet, E. Petit, T. Klockgether, A. Durr, T. Ashisawa, S. Tezenas du Montcel, and f. t. R. Consortium, "Temporal Dynamics of the Scale for the Assessment and Rating of Ataxia in Spinocerebellar Ataxias," en, *Movement Disorders*, vol. 38, no. 1, pp. 35-44, 2023.

[122] C. Murphy, "Olfactory and other sensory impairments in Alzheimer disease," en, *Nature Reviews Neurology*, vol. 15, no. 1, pp. 11-24, Jan. 2019.

[123] M. E. Murray, N. R. Graff-Radford, O. A. Ross, R. C. Petersen, R. Duara, and D. W. Dickson, "Neuropathologically defined subtypes of Alzheimer's disease with distinct clinical characteristics: A retrospective study," en, *The Lancet Neurology*, vol. 10, no. 9, pp. 785-796, Sep. 2011.

[124] T. Nedelec, B. Couvy-Duchesne, F. Monnet, T. Daly, M. Ansart, L. Gantzer, B. Lekens, S. Epelbaum, C. Dufouil, and S. Durrleman, "Identifying health conditions associated with Alzheimer's disease up to 15 years before diagnosis: An agnostic study of French and British health records," English, *The Lancet Digital Health*, vol. 4, no. 3, e169-e178, Mar. 2022.

[125] M. Nguyen, T. He, L. An, D. C. Alexander, J. Feng, and B. T. T. Yeo, "Predicting Alzheimer's disease progression using deep recurrent neural networks," en, *NeuroImage*, vol. 222, p. 117-203, Nov. 2020.

[126] S. S. O'Sullivan, A. H. Evans, and A. J. Lees, "Dopamine dysregulation syndrome: An overview of its epidemiology, mechanisms and management," eng, *CNS drugs*, vol. 23, no. 2, pp. 157-170, 2009.

[127] N. P. Oxtoby, A. L. Young, D. M. Cash, T. L. S. Benzinger, A. M. Fagan, J. C. Morris, R. J. Bateman, N. C. Fox, J. M. Schott, and D. C. Alexander, "Data-driven models of dominantly-inherited Alzheimer's disease progression," *Brain*, vol. 141, no. 5, pp. 1529-1544, May 2018.

[128] K. V. Papp, D. M. Rentz, I. Orlovsky, R. A. Sperling, and E. C. Mormino, "Optimizing the preclinical Alzheimer's cognitive composite with semantic processing: The PACC5," *Alzheimer's & Dementia: Translational Research & Clinical Interventions*, vol. 3, no. 4, pp. 668-677, Nov. 2017.
[129] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, and D. Cournapeau, "Scikit-learn: Machine Learning in Python," en, *MACHINE LEARNING IN PYTHON*,

[130] R. C. Petersen, P. S. Aisen, L. A. Beckett, M. C. Donohue, A. C. Gamst, D. J. Harvey, C. R. Jack, W. J. Jagust, L. M. Shaw, A. W. Toga, J. Q. Trojanowski, and M. W. Weiner, "Alzheimer's Disease Neuroimaging Initiative (ADNI): Clinical characterization," en, *Neurology*, vol. 74, no. 3, pp. 201-209, Jan. 2010.

[131] R. Pinder, "The Honeymoon Period—And After," en, in *The Management of Chronic Illness: Patient and Doctor Perspectives on Parkinson's Disease*, R. Pinder, Ed., London: Macmillan Education UK, 1990, pp. 54–67.

[132] R. B. Postuma, A. E. Lang, J. F. Gagnon, A. Pelletier, and J. Y. Montplaisir, "How does parkinsonism start? Prodromal parkinsonism motor changes in idiopathic REM sleep behaviour disorder," *Brain*, vol. 135, no. 6, pp. 1860–1870, Jun. 2012.

[133] P.-E. Poulet and S. Durrleman, "Mixture Modeling for Identifying Subtypes in Disease Course Mapping," en, in *Information Processing in Medical Imaging*, A. Feragen, S. Sommer, J. Schnabel, and M. Nielsen, Eds., ser. Lecture Notes in Computer Science, Cham: Springer International Publishing, 2021, pp. 571–582.

[134] P.-E. Poulet and S. Durrleman, "Multivariate disease progression modeling with longitudinal ordinal data," en, *Statistics in Medicine*, vol. n/a, no. n/a,

[135] *Projections de population à l'horizon 2060 - Insee Première - 1320.*

[136] C. Proust-Lima, V. Philipps, and B. Liquet, "Estimation of Extended Mixed Models Using Latent Classes and Latent Processes: The R Package lcmm," en, *Journal of Statistical Software*, vol. 78, pp. 1–56, Jun. 2017.

[137] C. Proust-Lima, V. Philipps, B. Perrot, M. Blanchin, and V. Sébille, "Modeling repeated self-reported outcome data: A continuous-time longitudinal Item Response Theory model," en, *Methods*, vol. 204, pp. 386–395, Aug. 2022.

[138] M. Quiroga, P. G. Garay, J. M. Alonso, J. M. Loyola, and O. A. Martin, *Bayesian additive regression trees for probabilistic programming*, Oct. 2022.

[139] L. L. Raket, "Statistical Disease Progression Modeling in Alzheimer Disease," *Frontiers in Big Data*, vol. 3, 2020.

[140] L. L. Raket, "Progression models for repeated measures: Estimating novel treatment effects in progressive diseases," en, *Statistics in Medicine*, vol. 41, no. 28, pp. 5537–5557, 2022.

[141] M. A. Rana, A. Li, D. Fox, B. Boots, F. Ramos, and N. Ratliff, "Euclideanizing Flows: Diffeomorphic Reduction for Learning Stable Dynamical Systems," en, in *Proceedings of the 2nd Conference on Learning for Dynamics and Control*, PMLR, Jul. 2020, pp. 630–639.

[142] D. Ravi, D. C. Alexander, and N. P. Oxtoby, "Degenerative Adversarial NeuroImage Nets: Generating Images that Mimic Disease Progression," en, in *Medical Image Computing and Computer Assisted Intervention - MICCAI 2019*, D. Shen, T. Liu, T. M. Peters, L. H. Staib, C. Essert, S. Zhou, P.-T. Yap, and A. Khan, Eds., ser. Lecture Notes in Computer Science, Cham: Springer International Publishing, 2019, pp. 164–172.
[143] S. L. Risacher, W. H. Anderson, A. Charil, P. F. Castelluccio, S. Shcherbinin, A. J. Saykin, A. J. Schwarz, and For the Alzheimer's Disease Neuroimaging Initiative, "Alzheimer disease brain atrophy subtypes are associated with cognition and rate of decline," en, *Neurology*, vol. 89, no. 21, pp. 2176-2186, Nov. 2017.

[144] H. Robbins and S. Monro, "A Stochastic Approximation Method," *The Annals of Mathematical Statistics*, vol. 22, no. 3, pp. 400-407, Sep. 1951.

[145] J. Salvatier, T. V. Wiecki, and C. Fonnesbeck, "Probabilistic programming in Python using PyMC3," en, *PeerJ Computer Science*, vol. 2, e55, Apr. 2016.

[146] F. Samejima, "Estimation of latent ability using a response pattern of graded scores," en, *Psychometrika*, vol. 34, no. 1, pp. 1-97, Mar. 1969.

[147] T. Saulnier, V. Philipps, W. G. Meissner, O. Rascol, A. Pavy-Le Traon, A. Foubert-Samier, and C. Proust-Lima, "Joint models for the longitudinal analysis of measurement scales in the presence of informative dropout," en, *Methods*, vol. 203, pp. 142-151, Jul. 2022.

[148] B. Sauty and S. Durrleman, "Progression Models for Imaging Data with Longitudinal Variational Auto Encoders," en, in *Medical Image Computing and Computer Assisted Intervention - MICCAI 2022*, L. Wang, Q. Dou, P. T. Fletcher, S. Speidel, and S. Li, Eds., ser. Lecture Notes in Computer Science, Cham: Springer Nature Switzerland, 2022, pp. 3-13.

[149] B. Sauty and S. Durrleman, "Riemannian Metric Learning for Progression Modeling of Longitudinal Datasets," in *2022 IEEE 19th International Symposium on Biomedical Imaging (ISBI)*, Mar. 2022, pp. 1-5.

[150] B. Sauty and S. Durrleman, "Impact of sex and APOE-4 genotype on patterns of regional brain atrophy in Alzheimer's disease and healthy aging," *Frontiers in Neurology*, vol. 14, 2023.

[151] J.-B. Schiratti, S. Allassonnière, A. Routier, O. Colliot, and S. Durrleman, "A Mixed-Effects Model with Time Reparametrization for Longitudinal Univariate Manifold-Valued Data," en, in *Information Processing in Medical Imaging*, S. Ourselin, D. C. Alexander, C.-F. Westin, and M. J. Cardoso, Eds., ser. Lecture Notes in Computer Science, Cham: Springer International Publishing, 2015, pp. 564-575.

[152] J.-B. Schiratti, "Methods and algorithms to learn spatio-temporal changes from longitudinal manifold-valued observations," en, Ph.D. dissertation, Université Paris Saclay (COmUE), Jan. 2017.

[153] J.-B. Schiratti, S. Allassonnière, O. Colliot, and S. Durrleman, "A Bayesian mixed-effects model to learn trajectories of changes from repeated manifold-valued observations," *The Journal of Machine Learning Research*, vol. 18, no. 1, pp. 4840-4872, Jan. 2017.

[154] J. Serroyen, G. Molenberghs, G. Verbeke, and M. Davidian, "Nonlinear Models for Longitudinal Data," *The American Statistician*, vol. 63, no. 4, pp. 378-388, Nov. 2009.

[155] C. Shand, P. J. Markiewicz, D. M. Cash, D. C. Alexander, M. C. Donohue, F. Barkhof, N. P. Oxtoby, and Alzheimer's Disease Neuroimaging Initiative, "Heterogeneity in Preclinical Alzheimer's Disease Trial Cohort Identified by Image-based Data-Driven Disease Progression Modelling," eng, *medRxiv: The Preprint Server for Health Sciences*, p. 2023.02.07.23285572, Feb. 2023.
[156] A. Sharma, V. Goyal, M. Behari, A. Srivastva, G. Shukla, and D. Vibha, "Impulse control disorders and related behaviours (ICD-RBs) in Parkinson's disease patients: Assessment using "Questionnaire for impulsive-compulsive disorders in Parkinson's disease" (QUIP)," *Annals of Indian Academy of Neurology*, 2015.

[157] S. L. Silverman, "From Randomized Controlled Trials to Observational Studies," en, *The American Journal of Medicine*, vol. 122, no. 2, pp. 114-120, Feb. 2009.

[158] T. Simuni, C. Caspell-Garcia, C. Coffey, S. Lasch, C. Tanner, and K. Marek, "How stable are Parkinson's disease subtypes in de novo patients: Analysis of the PPMI cohort?" en, *Parkinsonism & Related Disorders*, vol. 28, pp. 62-67, Jul. 2016.

[159] G. T. Stebbins, C. G. Goetz, D. J. Burn, J. Jankovic, T. K. Khoo, and B. C. Tilley, "How to identify tremor dominant and postural instability/gait difficulty groups with the movement disorder society unified Parkinson's disease rating scale: Comparison with the unified Parkinson's disease rating scale," en, *Movement Disorders*, vol. 28, no. 5, pp. 668-670, 2013.

[160] R. J. Swartz and S. W. Choi, "A Burdened CAT: Incorporating Response Burden With Maximum Fisher's Information for Item Selection," en,

[161] B. O. Taddé, H. Jacqmin-Gadda, J.-F. Dartigues, D. Commenges, and C. Proust-Lima, "Dynamic modeling of multivariate dimensions and their temporal relationships using latent processes: Application to Alzheimer's disease," en, *Biometrics*, vol. 76, no. 3, pp. 886-899, 2020.

[162] A. Tam, C. Laurent, S. Gauthier, and C. Dansereau, "Prediction of Cognitive Decline for Enrichment of Alzheimer's Disease Clinical Trials," *The Journal Of Prevention of Alzheimer's Disease*, 2022.

[163] C. E. Teunissen, I. M. W. Verberk, E. H. Thijssen, L. Vermunt, O. Hansson, H. Zetterberg, W. M. van der Flier, M. M. Mielke, and M. del Campo, "Blood-based biomarkers for Alzheimer's disease: Towards clinical implementation," en, *The Lancet Neurology*, vol. 21, no. 1, pp. 66-77, Jan. 2022.

[164] M. H. Tosin, G. T. Stebbins, C. Comella, C. G. Patterson, D. A. Hall, and t. S. S. Group, "Does MDS-UPDRS Provide Greater Sensitivity to Mild Disease than UPDRS in De Novo Parkinson's Disease?" en, *Movement Disorders Clinical Practice*, vol. 8, no. 7, pp. 1092-1099, 2021.

[165] M. Vandemeulebroecke, B. Bornkamp, T. Krahnke, J. Mielke, A. Monsch, and P. Quarg, "A Longitudinal Item Response Theory Model to Characterize Cognition Over Time in Elderly Subjects," *CPT: Pharmacometrics & Systems Pharmacology*, vol. 6, no. 9, pp. 635-641, Sep. 2017.

[166] C. S. Venuto, N. B. Potter, E. Ray Dorsey, and K. Kieburtz, "A review of disease progression models of Parkinson's disease and applications in clinical trials," en, *Movement Disorders*, vol. 31, no. 7, pp. 947-956, 2016.

[167] G. Verbeke and E. Lesaffre, "A Linear Mixed-Effects Model With Heterogeneity in the Random-Effects Population," *Journal of the American Statistical Association*, vol. 91, no. 433, pp. 217-221, 1996.
[168] J. Weickenmeier, M. Jucker, A. Goriely, and E. Kuhl, "A physics-based model explains the prion-like features of neurodegeneration in Alzheimer's disease, Parkinson's disease, and amyotrophic lateral sclerosis," en, *Journal of the Mechanics and Physics of Solids*, vol. 124, pp. 264-281, Mar. 2019.

[169] D. Weintraub, J. Koester, M. N. Potenza, A. D. Siderowf, M. Stacy, V. Voon, J. Whetteckey, G. R. Wunderlich, and A. E. Lang, "Impulse Control Disorders in Parkinson Disease: A Cross-Sectional Study of 3090 Patients," *Archives of Neurology*, vol. 67, no. 5, pp. 589-595, May 2010.

[170] E. Yang, M. Farnum, V. Lobanov, T. Schultz, R. Verbeeck, N. Raghavan, M. N. Samtani, G. Novak, V. Narayan, A. DiBernardo, and Alzheimer's Disease Neuroimaging Initiative, "Quantifying the pathophysiological timeline of Alzheimer's disease," eng, *Journal of Alzheimer's disease: JAD*, vol. 26, no. 4, pp. 745-753, 2011.

[171] L. Younes, "Diffeomorphic Learning," *Journal of Machine Learning Research*, vol. 21, no. 220, pp. 1-28, 2020.

[172] A. L. Young, R. V. Marinescu, N. P. Oxtoby, M. Bocchetta, K. Yong, N. C. Firth, D. M. Cash, D. L. Thomas, K. M. Dick, J. Cardoso, J. van Swieten, B. Borroni, D. Galimberti, M. Masellis, M. C. Tartaglia, J. B. Rowe, C. Graff, F. Tagliavini, G. B. Frisoni, R. Laforce, E. Finger, A. de Mendonça, S. Sorbi, J. D. Warren, S. Crutch, N. C. Fox, S. Ourselin, J. M. Schott, J. D. Rohrer, and D. C. Alexander, "Uncovering the heterogeneity and temporal complexity of neurodegenerative diseases with Subtype and Stage Inference," en, *Nature Communications*, vol. 9, no. 1, pp. 1-16, Oct. 2018.

[173] A. L. Young, J. W. Vogel, L. M. Aksman, P. A. Wijeratne, A. Eshaghi, N. P. Oxtoby, S. C. R. Williams, D. C. Alexander, and, "Ordinal SuStaIn: Subtype and Stage Inference for Clinical Scores, Visual Ratings, and Other Ordinal Data," *Frontiers in Artificial Intelligence*, vol. 4, 2021.

[174] W. J. Zetusky, J. Jankovic, and F. J. Pirozzolo, "The heterogeneity of Parkinson's disease: Clinical and prognostic implications," eng, *Neurology*, vol. 35, no. 4, pp. 522-526, Apr. 1985.

[175] X. Zhang, E. C. Mormino, N. Sun, R. A. Sperling, M. R. Sabuncu, B. T. T. Yeo, and t. A. D. N. Initiative, "Bayesian model reveals latent atrophy factors with dissociable cognitive trajectories in Alzheimer's disease," en, *Proceedings of the National Academy of Sciences*, vol. 113, no. 42, E6535-E6544, Oct. 2016.