models are usually built on continuous variables such as fluid biomarkers, imaging or aggregated cognitive scores. Their extension to discrete data involves using a logit<fnref n="13" /> or a probit model<fnref n="14" />. For ordinal data the cumulative probit model remains one of the gold standard for longitudinal modelling of discrete markers. One recent example of such model combine with survival analysis (joint modelling) is<fnref n="15" />.

Dealing with items has been at the heart of Item Response Theory (IRT), from which we will borrow tools to provide an extension of the disease course mapping framework. IRT has long been favoured to deal with this type of data since the first developments in the second half of the 20th century with the works of Lord and his collaborators<fnref n="16" />. IRT has mostly been used for cross-sectional data (even in Parkinson’s disease<fnref n="17" />), and extensions to longitudinal data<fnref n="18" />,<fnref n="19" />,<fnref n="20" /> are less common, even though neurodegeneration is best highlighted by the collection of repeated measurements. In IRT each individual *i* is assumed to have a latent trait $\theta_i$, which becomes a function of time (usually linear) in the case of longitudinal IRT. In disease modelling this latent trait stands for the disease progression score. Longitudinal IRT models have been used on medical cohorts to design new composite scores<fnref n="21" />, however the model has too few latent parameters (only one per individual) and thus does little for the understanding of heterogeneity. Multidimensional IRT models<fnref n="22" /> use a multidimensional latent parameter $\theta_i$ for each individual, thus allowing for more heterogeneity in the model, but these have not been applied to disease progression modelling to our knowledge.

A recent attempt<fnref n="23" /> at bridging IRT with other frameworks has been successful, where a Bayesian non-parametric approach was used to handle the dynamics of longitudinal observations. However the approach taken to model the longitudinal aspect was to discretize time and describe the evolution of the latent trait as unidimensional with an autoregressive scheme. We want to pursue this trend with the introduction of IRT to geometric mixed-effect models.

<footnote>pp. 963–974, Dec. 1982.</footnote>

<footnote n="13">¹³T. R. T. Have, A. R. Kunselman, E. P. Pulkstenis, et al., “Mixed Effects Logistic Regression Models for Longitudinal Binary Response Data with Informative Drop-Out,” *Biometrics*, vol. 54, no. 1, pp. 367–383, 1998.</footnote>

<footnote n="14">¹⁴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.</footnote>

<footnote n="15">¹⁵T. Saulnier, V. Philipps, W. G. Meissner, et al., “Joint models for the longitudinal analysis of measurement scales in the presence of informative dropout,” en, *Methods*, vol. 203, pp. 142–151, Jul. 2022.</footnote>

<footnote n="16">¹⁶F. M. Lord, *Applications of Item Response Theory To Practical Testing Problems*. New York: Routledge, Jul. 1980.</footnote>

<footnote n="17">¹⁷G. Gottipati, A. Berges, S. Yang, et al., “Item Response Model Adaptation for Analyzing Data from Different Versions of Parkinson’s Disease Rating Scales,” *Pharmaceutical Research*, vol. 36, Jul. 2019.</footnote>

<footnote n="18">¹⁸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.</footnote>

<footnote n="19">¹⁹M. Vandemeulebroecke, B. Bornkamp, T. Krahnke, et al., “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.</footnote>

<footnote n="20">²⁰C. Proust-Lima, V. Philipps, B. Perrot, et al., “Modeling repeated self-reported outcome data: A continuous-time longitudinal Item Response Theory model,” en, *Methods*, vol. 204, pp. 386–395, Aug. 2022.</footnote>

<footnote n="21">²¹L. Arrington, S. Ueckert, M. Ahamadi, et al., “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.</footnote>

<footnote n="22">²²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.</footnote>

<footnote n="23">²³A. Cremaschi, M. De Iorio, Y. Seng Chong, et al., “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.</footnote>
## 2.2 Binary model

### 2.2.1 Equation

Based on the logistic curves model we will now apply it to non-continuous data, starting from simple categorical data. The basic framework is when observations $y_{ijk}$ are binary: either 0 or 1. The model adaptation is rather simple: the observations $y_{ijk}$ are not treated as a noisy measure of the true trajectory $\eta_{ik}(\psi_i(t_{ij}))$ but as a realization of a random Bernoulli variable with probability $\eta_{ik}(\psi_i(t_{ij}))$.

$$y_{ijk} \sim \mathcal{B}(\eta_{ik}(\psi_i(t_{i,j}))) \qquad (2.2.1)$$

As in the framework of IRT, if we consider the value to be the result of one item in a test as right or wrong (0 or 1), we can evaluate the observed result as a probabilistic outcome relying on a latent disease state passing through a logit model.

### 2.2.2 Application to symptom's prediction

The binary version of the DCM model was motivated by a specific cohort for Parkinson's disease (PD) coming from the NS-Park consortium. This dataset was recording only symptoms, most of which were encoded as 0 or 1. The data at the time of study came from 25 PD expert centers in France and gathered 18018 patients with a total of 33798 visits. Patients were included around 68 years old and followed for a mean duration of 1.8 years. The inclusion criterion was a PD diagnosis, for which the mean age was 61. It is notable that patients were included quite late in the study: about 7 years after the diagnosis. However it is known that the dopaminergic treatment in PD is a very efficient symptomatic treatment, which leads to a singular period during the first years of the disease often referred to as the "honeymoon"<fnref n="24" />. During this period, patients almost return to a normal life thanks to the treatment. Only a few years later does the treatment start to show signs of limited efficiency. The end of the honeymoon period after which patients start to suffer from the disease is when various symptoms start to appear, which is what the NS-Park consortium aimed at studying. Figure 2.1 summarizes the frequency of symptoms among the population.

We did not provide the raw frequency per visit since some symptoms are repeated a lot for one patient across its multiple visits while others are appearing only once and are treated afterwards (case of hallucinations for instance). This last possibility is annoying for our model since the logistic is increasing, meaning that the frequency of symptom occurrence is assumed to increase with time. As many symptoms can be treated, it would require a model being able to take treatment into account. This will be the focus of the dedicated chapter on treatment models.

Our first attempt consisted in using the full range of symptoms with a Bernoulli noise model in a prediction setting. We separated the dataset into a training set of individuals and a test set of individuals so that there is no data leakage between the two. In the test set, for each individual we divide the visits of one individual in two by setting a virtual time separation: one set constitutes the "past" visits, upon which we can perform a personalization (i.e. estimation of individual parameters), the other set constitutes the "future" visits that we try to predict.

As shown in Figure 2.2, in the second row, the predictions seem to be very good. ROC AUC range from 0.66 to 0.97 depending on the symptoms. However this is quite deceptive as in most

<footnote n="24">²⁴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.</footnote>
<figure><img src="image_17.png" /><figcaption>Figure 2.1: Proportion of selected individuals presenting the symptom at least once. Hoehn-Yahr is a PD severity rating scale ranging from 0 to 5 but we binarized it as ≥ 2. RBD: rapid eye movement sleep behavior disorder. DDS: dopamine dysregulation syndrome.</figcaption></figure>

cases the predictions are quite easy: for most untreated symptoms, once there has been at least one occurrence before, they tend to be observed every time afterwards. This means that once the personalization is done on an individual where one symptom has already reached 1, all the values of this symptom after are very likely to be 1 also. This explains why the AUC looks so good.

In a more realistic setting, we changed the prediction task. For each symptom, we selected patients in the test set for which we could split visits between past and future so that the symptom of interest had *never* been observed in the past. In this setting we are using the model to try and predict the first occurrence of the symptom based on all the other symptoms. This task is much more interesting for clinical practitioners as it allows them to anticipate symptom appearance and therefore adapt treatment if needed.
<figure><img src="image_18.png" /><figcaption>Figure 2.2: Prediction of symptoms. First row: example of prediction for one individual; green dots are past visits used to estimate individual parameters (estimation is performed with the 42 symptoms, not only the one showed) and red ones are future visits that we try to predict using the personalized DCM trajectory (blue curve). Middle row: ROC curves for this prediction task. Bottom row: ROC AUC when predicting the first occurrence of a symptom. Performances are way lower for this task.</figcaption></figure>
The results are much worse as shown in the last row of Figure 2.2. The mean AUC over all symptoms is 0.64. This is to be compared to a very naive approach that would predict symptom occurrence with the median time of onset after diagnosis. Such an approach would barely have 0.54 mean AUC. We hope that these results showcase how hard this task is. We can distinguish two groups of symptoms:

* ROC AUC $\ge$ 0.6: including apathy, weariness or falls. These symptoms are partially predicted, which means that part of their variance has been captured by the model

* ROC AUC $\le$ 0.6: including dyskinesia, RBD or postural deformation. These symptoms are un predictable using our binary DCM model without additional information

As we start to lean towards trying to detect the first occurrence of a symptom, we need to mention that this is more suited to survival analysis (see<fnref n="25" /> for an overview of survival analysis) as we are actually dealing with events. Using the other symptoms for this task could lead to using joint modelling<fnref n="26" />. But survival analysis and joint modelling is outside the scope of our work in this chapter.

We propose an interesting practical approach to still use the DCM with logistic curves to predict first symptom occurrence. We first select symptoms for which this first occurrence is crucial to predict, and which is not likely to have occurred prior to study inclusion. This selection has been done after discussing the list of symptoms in depth with a Parkinson’s disease expert Dr. J-C. Corvol. We retained 26 symptoms out of the 42, and then modified the data observations by setting all observations to 1 after the first occurrence. In this setting the model curves $\eta$ can be interpreted as the probability of the symptom not having occurred yet. The results of the fit in this case are shown in Figure 2.3.

The average trajectory shows a clear inclusion bias (which is not enforced as showcased by the nociceptive pain curve), which is why we hypothesised that symptoms did not occur before the inclusion in the study. Otherwise we would have left censoring, i.e. we would need to take into account that the patients could have had the symptom earlier without the information whether the event occurred. These curves are similar to survival curves for each symptom. The model in this case is equivalent to fitting a logistic (but reversed) to the Kaplan-Meier curve with patients having been aligned temporally on a common disease timeline. The interpretation of the different slopes is showcased in the bottom plot of the figure. The slope defines the transition window between during which the symptom is most likely to occur for the first time, which is very different from the interpretation of $v$ in the continuous DCM case.

To conclude on binary DCM models, we see that the application to symptom data is mitigated. Due to the variable nature of the occurrence of symptoms combined with the possibility of treatment, the DCM model has a hard time predicting the symptoms, even with a multivariate setting. Indeed there is too few information coming from binary data, which undermines the personalization step and worsens predictions. One way to improve would be to include continuous measurements along the binary ones, allowing biomarkers to guide the forecast of symptoms’ first occurrence.

We need to avoid being lured into two pitfalls. The first being using a Bernoulli noise model when we are only interested in the first occurrence of a symptom. This situation warrants instead

<footnote n="25">²⁵C. Kartsonaki, “Survival analysis,” en, *Diagnostic Histopathology*, Mini-Symposium: Medical Statistics, vol. 22, no. 7, pp. 263–270, Jul. 2016.</footnote>

<footnote n="26">²⁶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.</footnote>
<figure><img src="image_19.png" /><figcaption>Figure 2.3: Average trajectory of the first occurrence model. Top plot: average trajectory shown for 10 of the 26 symptoms, with the inclusion date shown in red. Bottom row: average trajectory for 4 symptoms shown separately, with the "transition window" highlighted in red.</figcaption></figure>

a survival model, or a joint model. The second pitfall is a computational one. When trying to optimize a model using continuous and discrete measurements at the same time, maximizing the likelihood leads to an unbalanced fit on the two types of data. Indeed, the log-likelihood of the Bernoulli leads to a binary cross-entropy while Gaussian noise leads to a mean squared error loss. These two losses are not on the same scale (one log loss and one square loss), and when looking at the gradient of those two we have the formulas .0.38 .0.39 (in appendix), where the gradient of the cross-entropy is scaled by $\frac{y-\eta}{\eta(1-\eta)}$ while the gradient of the mean squared error is scaled by $\frac{y-\eta}{\sigma^2}$.
However as we mentioned, symptoms tend to be either present or absent, quickly forcing $\eta$ to values very close to 0 or 1. The scaling factor thus quickly turns in favor of the binary items, meaning that the model tends to fit better the binary items. This is counter-intuitive as continuous items are more reliable and should provide more information, and thus drive the fit. This multi-loss issue currently prevents us to have decent DCM models using different data types as input. One solution would be to implement an adaptive loss scheme, using profiling weight functions (PWF) as in$^{27}$.

<footnote>$^{27}$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 – MIC-CAI 2019*, D. Shen, T. Liu, T. M. Peters, et al., Eds., ser. Lecture Notes in Computer Science, Cham: Springer International Publishing, 2019, pp. 164–172.</footnote>
## 2.3 Ordinal model

### 2.3.1 Equation

Now let's assume the measurements lie on a discontinuous space with ordered values such as for a cognitive score. Without loss of generality we assume these values are integers between 0 and L included. As for the binary case, we will again consider that the model values belonging to the manifold (which is continuous) describe a form of probability (here through the cumulative distribution function). We now specify the logistic curves model for ordinal data similar to Samejima's model (also known as cumulative logit model) from IRT<fnref n="28" />:

$$ 
\begin{align} 
\forall l \in [|1, L|], \mathbb{P}(y_{ijk} \ge l) &= \eta_{ik}(\psi_{ik}^l(t_{i,j})) \nonumber \\
&= \left(1 + \left(\frac{1}{p_k} - 1\right) \exp\left(-\frac{v_k \psi_{ik}^l(t_{ij}) + w_{ik}}{p_k(1-p_k)}\right)\right)^{-1} \tag{2.3.1} \\
\psi_{ik}^l(t_{ij}) &= e^{\xi_i}(t_{ij} - t_0 - \tau_i) + t_0 - \sum_{m=1}^{l} \delta_k^m \nonumber \\
\forall l \in [|1, L|], \delta_k^l > 0 \nonumber
\end{align} 
$$

where $\psi_{ik}^l$ is the time reparametrization of individual i for the item k at level l. As in the DCM model, the population parameters ($p_k, v_k, t_0$) define the base logistic curve (i.e. the curve for level 1 of the item) as the initial position, initial velocity and initial time respectively. $\tau_i$ and $\xi_i$ are the time-related parameters, namely the time-shift and the log-acceleration. $w_{ik}$ is the space-shift. $\delta_k^l$ are new parameters introduced to deal with the different levels of the item, as explained thereafter.

The idea is to model the cumulative distribution function with the logistic curves instead of the probabilities of each level directly. Each level is modelled with a logistic curve. However we do not want to overparametrize the model so we choose to enforce parallelism by setting the same velocity v for each level of a same item. This introduces only one parameter for each added level. The new parameter is $\delta_k^l$, the delay in time between levels $l-1$ and $l$. This delay is a fixed effect, so it is shared by the whole population. Since it is added to the time reparametrization function $\psi_{ik}$ after the affine transformation of time by $\tau$ and $e^\xi$, this time delay is a duration in the common disease timeline. This means that if one patient has an acceleration factor of 2, the expected time to jump from one level to the next is divided by 2. With the $\delta_k^l$ being positive, we also ensure that the condition $\mathbb{P}(y_{ijk} \ge l) < \mathbb{P}(y_{ijk} \ge l-1)$ is verified. From an interpretability standpoint, these delays are helpful in understanding which levels of an ordinal scale last longer than the others, or on the contrary which levels represent only a short transition.

### Bayesian modelling

We will specify the probabilistic model in order to estimate the parameters by maximizing the likelihood. We set the priors as in the DCM for $\tau_i, \xi_i, s_i, v, g, \beta$. In the ordinal model we added a new parameter: the $\delta$ parameters, which must be positive, so we set their prior distribution as a log-normal.

<footnote n="28">²⁸F. Samejima, “Estimation of latent ability using a response pattern of graded scores,” en, *Psychometrika*, vol. 34, no. 1, pp. 1–97, Mar. 1969.</footnote>
As mentioned before, the probabilistic model allows us to deal with missing observations without needing to impute them as the likelihood is only computed for the points we observe. The robustness to the missing data as been evaluated in$^{29}$.

In terms of optimization of the log-likelihood, changing the model from a Gaussian noise to a Bernoulli variable amounts to optimizing a crossentropy between the latent model and the trajectory instead of a mean squared error:

$$ \hat{y}_{i,j} = \eta_i(\psi_i(t_{i,j})) \qquad (2.3.2) $$

[Simplified notation]

$$ \log(p(\mathbf{y}|\mathbf{z}; \theta)) = -N_{tot} \log(\sigma\sqrt{2\pi}) - \frac{1}{2\sigma^2} \sum_{i=1}^{n} \sum_{j=1}^{N_i} \|\mathbf{y}_{i,j} - \hat{\mathbf{y}}_{i,j}\|^2 \qquad (2.3.3) $$

[Continuous model with Gaussian noise]

$$ \log(p(\mathbf{y}|\mathbf{z}; \theta)) = \sum_{i=1}^{n} \sum_{j=1}^{N_i} \sum_{k=1}^{d} y_{i,j,k} \log(\hat{y}_{i,j,k}) + (1 - y_{i,j,k}) \log(1 - \hat{y}_{i,j,k}) \qquad (2.3.4) $$

[Categorical model with Bernoulli realization]

For the ordinal model, maximizing the log-likelihood is straightforward, as we can simply come back to $\mathbb{P}(y_{ijk} = l)$ by a difference $\mathbb{P}(y_{ijk} \ge l) - \mathbb{P}(y_{ijk} \ge l + 1)$. One thing to note here is that the log-likelihood cannot be expressed in the curved exponential family and we are outside the scope of the theoretical guarantees for the MCMC-SAEM convergence. However, convergence has always been reached in practice.

### Initialization method

Since the estimation algorithm will only converge towards a local optimum of the likelihood, we need to ensure that either we replicate the estimation with several initial points or we provide a starting point which is already a "good guess". For the population parameters of the logistic curves model, we compute the initial point for $(\mathbf{p}, \mathbf{v})$ by computing a regression on the individual progression above the first level, with $\mathbf{p}$ being deduced from the intercept and $\mathbf{v}$ from the slope. For the $(\delta_k^l)_{k,l}$, we computed them as the mean time to reach level $l$ from level $l-1$. Indeed each $\delta_k^l$ can be understood as the time between the progression of the probabilities $\mathbb{P}(y_{ijk} \ge l)$ and $\mathbb{P}(y_{ijk} \ge l-1)$.

### Implementation details

All the code for the ordinal disease course mapping model is already available on Gitlab: https://gitlab.com/icm-institute/aramislab/leaspy

### 2.3.2 Experiments: Simulation study

**First part:** In this first experience we generate synthetic ordinal data in order to highlight two benefits of our modelling approach. On the one hand, we show how modelling items leads to a

<footnote>$^{29}$R. Couronné, M. Vidailhet, J.-C. Corvol, et al., "Learning disease progression models with longitudinal data and missing values," in ISBI 2019 - International Symposium on Biomedical Imaging, Venice, Italy, Apr. 2019.</footnote>
finer-grained description of the disease progression, as compared to modelling aggregated scores. In particular, we show that a model based on item response exhibits patients clusters that were not visible with a continuous model. On the other hand, we compare how the ordinal model fares as opposed to the logistic curves model. We show that the flexibility introduced with the $\delta$ parameters allows to learn a step function which is less rigid than an imposed template such as a logistic curve.

<figure><img src="image_20.png" /><figcaption>Figure 2.4: Average trajectory $\gamma_0$ for the ordinal model. Top row: for each item the cumulative distribution functions are shown. Note that $\delta_k^l$ correspond to the shift on the x-axis between the logistic curves. Middle row: raw probabilities are shown, directly obtained by difference between the logistic curves of top row. Last row: items are shown together by summarizing the multinomial distributions of items. On the left the value at each timepoint is the maximum likelihood; on the right the value is the expectation. We also plot the aggregated score computed as the mean of the items.</figcaption></figure>

*Generation process:* We use our ordinal model in order to simulate synthetic patients and obtain their observations at randomly chosen timepoints. For the initialization, we chose the statistical model parameters $\theta$ and the population parameters $z_{pop}$. The values are chosen arbitrarily in the range of plausible values for a neurodegenerative-like disease progression model. In this simulation part, we kept a low dimensionality for the data with only $d = 3$ dimensions, which is the lowest number of dimensions that allows a meaningful ICA decomposition with at least 2 sources. We then adjusted the $\delta_k^l$ so that each of the three items has an atypical progression, by which we mean it can not be easily interpolated by a linear curve or a logistic curve. Each item has $L = 4$ levels, which sums into an aggregated score ranging from 0 to 12. Figure 2.4 shows the average trajectory $\gamma_0$ of
the model, as well as the average aggregated trajectory. Once the population parameters are fixed, we randomly generate new individuals: for each individual parameter, we sample the distribution mentioned in section 2.1.2. However to generate two different clusters, we choose to randomly attribute a class 1 or 2 to a new individual. Then for the individual parameters simulation, the only difference between individuals of the two clusters is that the first source in cluster 1 follows a Gaussian distribution centered around -4 instead of 0 with a standard deviation of 1 and the source 1 in cluster 2 is oppositely centered around +4. This source acts as an inverter of the order of the curves. The mean trajectory in cluster 1 has the curve Y0 starting first, then Y1 follows and finally Y2, whereas for cluster 2 the order is switched: Y1 first, then Y2 and finally Y0. However the aggregated score is roughly the same in the two clusters. With an equal repartition of individuals between the two clusters, we therefore obtain a balanced dataset of individuals. Once this is done, we need to generate random timepoints for the visits of each individual. For this we decided to select the number of visits with a binomial distribution with probability p = 0.5 and n = 10 to which we add 2 to guarantee at least two visits per individual. Then each of the timepoint is sampled with a Gaussian distribution around the middle of the disease progression curve, with a time window calibrated to span about 5 years of follow-up for each subject.

We estimated four models on this synthetic dataset:

* a continuous model (equation 1.1.4) on items: to highlight the purpose of using a dedicated ordinal model for such scales, we fit a disease mapping course model with logistic curves as presented in the method section with observations considered as continuous with a Gaussian noise. The model is three dimensional with two sources

* an ordinal model (equation 2.3.1) on items: we fit the model presented in the method section, here the observations are treated as ordinal. The model is also three dimensional with two sources

* a univariate continuous model on the aggregated score: to show how the loss of item information impacts the ability of a model to produce relevant information about patients, we fit a single logistic curve model on the sum of the three items. The model is thus one-dimensional without sources

* a univariate ordinal model on the aggregated score

The first point of our synthetic experiment was to show how the gain of finer-grain information allows to learn meaningful data patterns. We extracted the individual parameters of the four models. Each individual $i$ thus has a set of parameters $(\tau_i, \xi_i, s_i^1, s_i^2)$ in the multivariate models and $(\tau'_i, \xi'_i)$ in the univariate models. For visualizing the points in Figure 2.5, we decided to embed the individual in a two-dimensional space using a t-distributed stochastic neighbor embedding<fnref n="30" />. We see that the two clusters are clearly separated in the multivariate models, whereas the aggregated model is not able to disentangle them. A quantitative classification with a logistic regression taking in input the set of individual parameters and predicting the cluster of origin confirms this conclusion: with $(\tau_i, \xi_i, s_i^1, s_i^2)$ from either multivariate model the classifier gets close to 100% accuracy while with $(\tau'_i, \xi'_i)$ from the best univariate it only has 65%, which is much worse.

While the multivariate continuous model also captures individual heterogeneity by identifying the two clusters, it is not as performant in terms of pure data reconstruction. Figure 2.6 illustrates

<footnote n="30">³⁰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.</footnote>
<figure><img src="image_21.png" /><figcaption>Figure 2.5: t-SNE plots for individual parameters. On the left side the parameters in the multi-variate model are shown, on the right side the parameters in the univariate models. Top row shows the two ordinal models while bottom row shows the continuous models. Each color corresponds to one of the clusters (1 in blue, 2 in red).</figcaption></figure>

the difference between the continuous and the ordinal model in terms of reconstruction errors on the items, with two different methods for predicting the outcome of the ordinal model. It is clear that the continuous model is worse, with a first burden being the lower bound of the noise due to the discretized scale (every 0.25 in this case). This can be mitigated in part by comparing the rounded predictions of the continuous model to the maximum likelihood predictions of the ordinal model, but even in this case we average 0.38 points of absolute error for continuous rounded predictions versus 0.34 for the ordinal maximum likelihood predictions. The second burden is the prior choice of the form of the curve (here the logistic). It is important to note that for small scale items in general there is hardly a perfect continuous curve due to the variance of item levels. For instance, a standard Likert scale for a symptom would be a 0-4 scale with each level corresponding to a severity measure: null, mild, moderate, severe, very severe. The difference between two successive levels may vary, which highlights the role of the $\delta_k^l$ parameters.

These preliminary results encourage the use of dedicated ordinal models on ordinal scales, and this proves true even for items with more levels (especially when levels are unevenly distributed).

**Second part:** In the first part of the synthetic experiment we generated data using an ordinal
<figure><img src="image_22.png" /><figcaption>Figure 2.6: Error of reconstruction scatter plot. For each individual, we computed the mean absolute error of reconstruction over all visits and all items between the observations and the model prediction for the multivariate models. In the case of the continuous model the prediction is directly the outcome of the model, whereas for the ordinal model we computed the maximum likelihood for each value in the left plot, and the expectation in the right plot. The plots show that the ordinal model outperforms the continuous model for almost all patients in the reconstruction of the items.</figcaption></figure>

model and showed that the ordinal model was better at reconstructing such data than a continuous model. It is important to check that this also holds for data not generated by the ordinal model.

We used a continuous model similar to the one used in the first part, but without the $\delta_k^l$ parameters. For each individual $i$ generated, at each visit $j$ for each item $k$ we compute a noisy observation $y_{ijk} \in [0, 1]$. We then transform this data into ordinal data. We choose a Likert scale (0-4), so we multiply the observation by 4 and round it. Both multivariate models, ordinal and continuous, are calibrated on the data generated and we compute the reconstruction errors as in Figure 2.6. Results are presented in Figure 2.7 below.

We provided experimental results with data generated with a continuous model rather than the ordinal model. These results show that the ordinal model is also able to outperform a continuous model even in this setting.

### 2.3.3 Application to Parkinson's disease data

#### Data set

**Data availability statement:** Data used in the preparation of this article were obtained from the Parkinson's Progression Markers Initiative (PPMI) database<fnref n="31" /> (www.ppmi-info.org/access-data-specimens/download-data). For up-to-date information on the study, visit ppmi-info.org. <fnref n="32" />

The cohort includes mainly Parkinson's disease (PD) patients within two years of their diagnosis. There are also prodromal individuals, but we excluded those as they are a minority and some of

<footnote n="31">³¹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.</footnote>

<footnote n="32">³²PPMI – a public-private partnership – is funded by the Michael J. Fox Foundation for Parkinson's Research and funding partners, including [list the full names of all of the PPMI funding partners found at www.ppmi-info.org/about-ppmi/who-we-are/study-sponsors].</footnote>
<figure><img src="image_23.png" /><figcaption>Figure 2.7: Error of reconstruction scatter plot. For each individual, we computed the mean absolute error of reconstruction over all visits and all items between the observations and the model prediction for the multivariate models. In the left plot the predictions of the ordinal model are the maximum likelihood estimation while the predictions of the continuous model are rounded. In the right plot for the continuous model the prediction is directly the outcome of the model, whereas for the ordinal model we computed the expectation. The right plots show that the ordinal model outperforms the continuous model for almost all patients in the reconstruction of the items.</figcaption></figure>

them do not convert to PD. We selected all PD diagnosed patients with at least 2 visits in the study in order to have a longitudinal history for our training set, which left 900 subjects for a total of 7918 visits. From the wide range of biomarkers and clinical assessments we chose to focus on the MDS-UPDRS score. It is composed of 65 items divided in 4 sections: non-motor aspects of experiences of daily living (I) which is filled half by the clinician and half by the patient in a self-report questionnaire; motor aspects of experiences of daily living (II) which is a form only filled by the patient; motor examination (III) by the clinician; motor complications (IV) which is assessed by the clinician.

Each item is rated on a Likert scale (0: Normal, 1: Slight, 2: Mild, 3: Moderate, 4: Severe). For the third part, two versions of the test are recorded depending on the treatment effect. A version "OFF" is assessed when the patient hasn't taken medication in the last 6 hours, which is a state when the treatment should have little effect and motor symptoms will usually be observed. A version "ON" corresponds to the same motor examination shortly after the patient has taken medication, which is a state where motor complications are alleviated. Therefore we used "OFF" measures when selecting features in order to build the disease natural history with as few mitigation of treatment effect as possible. We also removed the fourth part of the MDS-UPDRS since it focuses on treatment secondary effects. In the end we obtained 59 items rated on scales from 0 to 4.

Most of the items have a low mean: 44 items have a mean between 0.4 and 1.25. Most of the observed values are 0 (66%), with 22% of 1 values, and only 12% of values above 2. For some items, the higher levels are very seldom seen so their progression cannot be modelled fully due to lack of data. The average total score on the parts I, II and III is 28.7 ± 22.9 with 90% of the values within the [1, 60] range.
Qualitative experiment: a description of the MDS-UPDRS

As shown in the previous summary, patients in PPMI do not display all possible symptoms, and most of their items remain at 0 throughout their follow-up. However, items with non-zero values often occur in groups, suggesting that several items capture a similar biological phenomenon. Our goal is to unravel these trends using Independent Component Analysis (ICA) on the space-shifts.

We calibrated an ordinal DCM model on the 900 subjects using their complete set of 59 items. Missing values do not need to be imputed since the model only computes the likelihood for observed values, which is one of the advantages of using DCM. The estimation was performed with the MCMC-SAEM algorithm, using a block Gibbs sampler. The sampling was grouped for all 59 dimensions of the population parameters **v**, **g**, **δ** during the first 16,000 iterations. For the final 4,000 iterations, a more refined Gibbs sampler was used, sampling each scalar coordinate separately. The samplers employed an oscillating tempered scheme 19. The overall calibration time was 56 hours.

Due to the high number of dimensions, we chose 8 as the number of sources to reduce the dimension of individual parameters and prevent overfitting. Previous trials with fewer sources resulted in significantly lower log-likelihood at the end of the calibration, while increasing the number of sources became too computationally expensive. ICA enforces correlations within the space-shifts due to the relation $w_i = As_i$, with **A** being a rectangular matrix. We can analyze the correlations across individuals resulting from these space-shift values, which can be seen as a time delay in the progression of each item compared to the population's average trajectory. Figure 2.8 shows the Pearson correlation matrix between the space-shifts corresponding to each item. We performed hierarchical clustering to isolate groups of highly correlated items.
<figure><img src="image_24.png" /><figcaption>Figure 2.8: Pearson correlation matrix of space-shifts at the item level. Hierarchical clustering has been performed to order the items of the different parts of the MDS-UPDRS.</figcaption></figure>
As we can see in Figure 2.8, the first group of items is comprised of tremor related items. The next two clusters gather items of the motor part (III) of the MDS-UPDRS, with one cluster containing items for the left part of the body and the other being dedicated to the right part. This is explained by the fact that Parkinson’s disease very often declares itself in one side of the body before the other. Then we have a small trio of gait, freezing of gait and posture, which can be included in a larger cluster of items concerning motor aspects of daily life (essentially items of the part II related to movement and balance). Then the last large cluster includes all the elements of part I (non-motor aspects of daily life) and the last items of part II. This large cluster can be separated into two subclusters, one being the part I of the MDS-UPDRS and the second being the leftover items of part II which are more related to mouth tasks (eating, drooling...) and communication. Overall we can notice a strong correlation between items of part I and part II. These parts are both questionnaires and mostly filled by the patient, which can explain the correlation bias due to the self-assessment of one’s abilities.

Parkinson’s patients are commonly separated into two or three subtypes<fnref n="33" />: tremor dominant (TD), postural instability and gait dominant (PIGD) and mixed patients (at the limit between TD and PIGD). Here we see that the TD and PIGD tendencies are represented in the correlations of items of the first and fourth clusters.

Some items were not regrouped with others or were associated in different ways. For instance we can note that neck rigidity and facial expression were at the junction between the left and right motor clusters, suggesting that higher body motor impairment usually follows motor symptoms in one side of the body. However these two items are less correlated to motor impairment in the legs and the feet. This is very interesting as we can hypothesise that the disease implies a localised neurodegeneration of motor neurons, and then symptoms spread from their original onset, explaining the anticorrelation between top and bottom motor impairment. Speech (part III question 1) has been isolated, while spontaneity of movement (part III question 14) seems to correlate with all the motor assessments of part III, as a global appreciation of the motor state of the patient by the clinician.

We gave a specific look at the item called DDS (dopamine dysregulation syndrome) which is a disturbance of dopamine therapy with addictive patterns such as impulse control disorders (ICD) which include gambling, hypersexuality or compulsive eating<fnref n="34" />. This syndrome is a consequence of too much dopamine uptake (sometimes voluntarily on the patient’s side). What appears is that the DDS is anticorrelated to all motor items of part III *except* for the freezing of gait (part III question 11). This backs the conclusion that freezing of gait is also a secondary effect of dopaminergic treatment rather than a disease symptom<fnref n="35" />.

We then examined the sources individually to identify which statistically independent linear combinations of the items were found, and tried to understand if they were linked to a pathological trend. We provide a description of the 8 sources:

* source 1 (variance explained 11%): this source induces a positive delay in all motor items of the left part of the body while it pushes for a negative delay in almost all of the other items.

<footnote n="33">³³G. T. Stebbins, C. G. Goetz, D. J. Burn, et al., "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.</footnote>

<footnote n="34">³⁴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.</footnote>

<footnote n="35">³⁵A. Espay, A. Fasano, B. van Nuenen, et al., ""On" state freezing of gait in Parkinson disease," *Neurology*, vol. 78, no. 7, pp. 454-457, Feb. 2012.</footnote>
We can interpret this source as a left body delay

* source 2 (variance explained 10%): this source is almost the same as source 1 but for right body items. We call this one right body delay. These two first sources with equivalent prevalence show that Parkinson's disease is as likely to start in the left part of the body as the right

* source 3 (variance explained 12%): this source has a positive delay on all motor items linked to the limbs' movement as opposed to the other items. We name this one motor source

* source 4 (variance explained 15%): this source has a positive delay on the left part of the body and the rigidity items. This source embeds a left body and rigidity combined progression

* source 5 (variance explained 12%): this source has a positive delay on all the questions that are self-evaluated as opposed to all the motor evaluations and the non-motor assessment by the clinician. We can therefore say that this source is the patient's self report effect

* source 6 (variance explained 13%): this source regroup items of the part I (non-motor aspects of daily life) with tremor manifestations as opposed to motor issues. This source is seen as the brain-first source

* source 7 (variance explained 14%): this source is driven mostly by the motor items and more especially the gait, posture and balance items. We call this one the gait and balance source

* source 8 (variance explained 13%): this source gathers items from the part I with a strong emphasis on anxiety, depression and apathy, along with speech issues and opposed to tremor items. We coin this one the cognitive source

We observe that some sources are in practice not totally independent, for instance source 5 (self report effect) and source 4 (left body rigidity) are correlated. This suggest that the number of sources is not fully adequate. We see that sometimes a source seem to include several groups of items at the same time (source 4 is a mix of left body items and rigidity), which could be the consequence of too few dimensions to disentangle these effects.

If we come back to the classification of patients between tremor dominant and postural insta-
bility/gait dominant, we see that PIGD is solely driven by source 7, whereas tremor information is
spread across sources 1 to 4.

We won't dwell on more clinical subtleties but hope this showcased the power of the ICA on
the spatial effects unfolding from the Riemannian framework of the DCM. The results here could
only be obtained thanks to the modelling of items with an ordinal method. The global model
calibrated on the 59 items showed that some items had a very low progression profile (or could
not be completely learned due to the few data available on the last levels), especially in the motor
examination part. However the second part of the MDS-UPDRS scale was more reliable since it
rarely stayed at zero and the progression was faster and steadier. We thus used the 13 items of the
second part in the next experiment to show quantitatively how well the ordinal model performs.

Quantitative results: prediction task

In order to understand the value brought by our new take on the model, we compare the ordinal
DCM to the logistic version of the DCM taking the cognitive scores as continuous values. We also
compare to a linear mixed-effect model and a no-change prediction model as our baseline. We will have two different options in order to show what the model can bring to the table: one model for the total score of one MDS part (the second part as mentioned previously) which is thus an integer between 0 and 42; and one model taking the items into account, as 42 is the sum of 13 questions with a scale from 0 to 4. In this last case we will show that taking the score as a continuous value is not valid since the noise here is lower than the step between possible values. We present the several models used in this experiment:

### Continuous univariate model

The MDS UPDRS part II (/42) is normalized between 0 and 1 in order for the logistic model to work. We then maximize the log-likelihood with Gaussian noise, thus implying a fit in the least squares sense.

### Continuous multivariate model

Each item (/4) of the MDS UPDRS part II is normalized between 0 and 1. As for the univariate maximize the log-likelihood with Gaussian noise, where the standard deviation of the noise is feature-dependent. For the prediction of the MDS UPDRS II total we sum the predictions of the items.

### Ordinal models

In order to operate a fair comparison we use the same hyperparameters as for the continuous version, so that the ordinal model shares the same parameters as the continuous one except for the $\delta$. We fit by maximizing the log-likelihood, and we compute the predictions using the expectation of the ordinal model: $E(\eta_{ij}) = \sum_l P(\hat{y}_{ij} \ge l)$. As for the continuous case, we have a univariate model built with only the MDS UPDRS II total score and a multivariate model using the 13 items. For the multivariate case, we compute the prediction as the expectation of the sum of the items.

### Baseline models

We also included a baseline constant model that predicts the last known value for the patient (hypothesising no change in the future) and a linear mixed-effect model for comparison. Note that they perform very well: it is hard to beat these models when predicting in a very short future since the state of the patient evolves very slowly.

To sum up, we provide below the formulas of the models used.

The continuous univariate model computes the MDS-UPDRS part II total as:

$$y_{ij} = \left( 1 + \left( \frac{1}{p_0} - 1 \right) \exp \left( - \frac{v_0 \left( e^{\xi_i} (t_{ij} - t_0 - \tau_i) + t_0 \right)}{p_0 (1 - p_0)} \right) \right)^{-1} + \epsilon_{ij}$$

with $\epsilon_{ij} \sim N(0, \sigma)$.

The multivariate continuous model computes the 13 items of the MDS-UPDRS part II as follows:

$$\forall k \in [1, 13], y_{ijk} = \left( 1 + \left( \frac{1}{p_k} - 1 \right) \exp \left( - \frac{v_k \left( e^{\xi_i} (t_{ij} - t_0 - \tau_i) + t_0 \right) + w_{ik}}{p_k (1 - p_k)} \right) \right)^{-1} + \epsilon_{ijk}$$

with $\epsilon_{ijk} \sim N(0, \sigma)$.

The ordinal univariate model computes the MDS-UPDRS part II total as:

$$\forall l \in [1, 52], P(y_{ij} \ge l) = \left( 1 + \left( \frac{1}{p_0} - 1 \right) \exp \left( - \frac{v_0 \left( e^{\xi_i} (t_{ij} - t_0 - \tau_i) + t_0 - \sum_{m=1}^{l} \delta_0^m \right)}{p_0 (1 - p_0)} \right) \right)^{-1}$$
The multivariate ordinal model computes the 13 items of the MDS-UPDRS part II as follows:

$$ \forall k \in [1, 13], \forall l \in [1, 4], \mathbb{F}(y_{ijk} \ge l) = \left( 1 + \left( \frac{1}{p_k} - 1 \right) \exp \left( - \frac{v_k(e^{\xi_i}(t_{ij} - t_0 - \tau_i) + t_0 - \sum_{m=1}^l \delta_k^m) + w_{ik}}{p_k(1 - p_k)} \right) \right)^{-1} $$

The constant prediction model for the MDS-UPDRS part II total formulates as $\hat{y}_{ij} = y_{ij,final}$ where $j_{final} < j$ is the last visit known for subject i in the training data.

The linear mixed-effect model for the MDS-UPDRS part II total formulates as :

$$ y_{ij} = \beta + \alpha \times t_{ij} + \beta_i + \alpha_i \times t_{ij} + \epsilon_{ij} $$

with $\beta$ being the fixed intercept, $\alpha$ the fixed slope, $\beta_i \sim N(0, \sigma_\beta)$ the random intercept, $\alpha_i \sim N(0, \sigma_\alpha)$ the random slope and $\epsilon_{ij} \sim N(0, \sigma)$ the Gaussian noise.

<figure><img src="image_25.png" /><figcaption>Figure 2.9: Mean absolute prediction error as a function of the current patient state. Bin ranges correspond to the MDSII value of the patient. Error bars for mean model error are estimated from cross-validation. LMM: linear mixed-effect model; continuous: disease course mapping model; ordinal: disease course mapping model, ordinal version; univariate: model learned of MDSII total directly; multivariate: model learned on MDSII items and then aggregated.</figcaption></figure>

The prediction framework we use consists in two steps. First a training set is used to calibrate a model, i.e. learning the fixed effects for the mixed-effect models (linear and DCM). In a second step, we use the patients of the test set. Out of each patient's set of visits we randomly select some past visits and a visit in the future for which we want to predict. We learn the random effects of the models on the past visits and use the personalized model to compute the value at the desired future visit's timepoint. We use a repeated 8-fold cross-validation to measure the models errors. The Figure 2.9 shows the mean absolute prediction error.
As mentioned previously, we observe that constant prediction and linear mixed model perform best when the patient score is low, as they progress slowly at the beginning of the disease. However we see that the prediction from those simple models worsens as the score increases, meaning that linear and constant models are less able to capture disease progression. As a matter of fact, they largely under estimate the change for patients later in the disease.

On the ordinal side we observe a similar trend, although both ordinal models are more reliable at higher disease stages. They outperform continuous models on the low stages, due to the flexibility of the step parameters $(\delta_l)_l$. However when it comes to later disease stages continuous models catch up with ordinal models in performance. This is mainly due to the data being much sparser at these stages, leaving a larger role to the shape priors. In the long term setting it confirms that the logistic prior improves on linear and constant curve priors. Ordinal models show that their low granularity allows for better prediction at low levels when changes are small, while we lack data for higher disease stages thus leading to a mitigated performance.

Comparing univariate versus multivariate models is very interesting. At early disease stages the univariate models seem to better capture the essence of the progression. On the other hand, the higher the disease stage the better the multivariate models perform, showing that multidimensional models are better suited for disease progression prediction when the rate of change increases.

Figure 2.11 shows the average disease progression learned by univariate mixed-effect models. Spaghetti plots with wrapped individual trajectories are provided below. The Figure 2.11 shows that the logistic and the ordinal provide a somehow similar description of the aggregate score up until age 85. After the logistic model is bound to the assumption that the score will evolve to 1, which in practice is highly unlikely as it is a complete state of degeneracy and patients usually do not reach this stage. This is why the ordinal curve slows compared to the logistic at the end. The logistic model assumes that the disease progresses faster after 85 years old. On the other hand the linear mixed-effect model tends to underestimate the rate of change at this stage. Note that the shape of the step function is mostly driven by the $\delta$ values learned, and thus is more flexible. This explains why it is more accurate for the values of the MDSII total between 0 and 15 where most of the observations lie according to the histograms.

We provide a visualization of the fit of univariate models with Figure 2.10. The spaghetti plots show the individual trajectories mapped onto the average trajectories. We do this by removing the random effects, hence only leaving the noise of individual observations as a deviation from the mean trajectory, shown in black, which accounts for the fixed effects. Note that the score we are modelling is quite noisy, with some individuals experiencing very sharp and sudden surges or collapses. This is also one of the reasons for the use of an ordinal model: the noise is more flexibly modelled as a spread of probability weights across several levels rather than a noise distribution chosen in model design (here a simple Gaussian noise).

The experiment was performed again on parts I and III of the MDS-UPDRS with continuous and ordinal models, and the results were essentially the same:

* For the MDS-UPDRS part I (/42), the ordinal univariate had a mean absolute error of 4.64 while the multivariate ordinal had 4.54; continuous univariate was 4.76 and continuous multivariate was 5.41

* For the MDS-UPDRS part III (/132), the ordinal univariate had a mean absolute error of 12.8 while the multivariate ordinal had 12.3; continuous univariate was 15.8 and continuous multivariate was 13.9
<figure><img src="image_26.png" /><figcaption>Figure 2.10: Spaghetti plots with individuals mapped onto the mean trajectory. Mapped trajectories are obtained by removing the random effects learned by each model. Top left: raw data; bottom left: linear mixed-effect model; top right: ordinal univariate model; bottom right: continuous univariate DCM model with logistic curve.</figcaption></figure>

What we seem to notice though is that the more items in a score, the better the multivariate approach performs, although it comes at a higher computational price.
<figure><img src="image_27.png" /><figcaption>Figure 2.11: Top figure: average disease progression of the univariate mixed-effect models as a function of the age. Bottom figure: histograms of patients mapped onto the corresponding timeline of their models. Aligned with the top figure it shows where data points are distributed according to each model. Parkinson age corresponds to that mapping on the time axis, as a result of the individual random effects linked to time ($$(\xi_i, \tau_i)$$ for ordinal and continuous DCM, random slope and intercept for LMM). LMM: linear mixed-effect model; continuous: disease course mapping model; ordinal: disease course mapping model, ordinal version.</figcaption></figure>
### 2.3.4 Discussion

Our work shows the potential of leveraging items with a small ordinal scale over working with an aggregated score. We built upon a non-linear mixed-effect model used mainly for continuous markers and extended it to ordinal data. This allowed to address cognitive decline by processing items of the test instead of their aggregated value. We showed how this could help understand disease trends in Parkinson's disease within the MDS-UPDRS. We believe this could be extended to other cognitive scales in other neurodegenerative scales with the aim to better understand items dynamics.

In the prediction task, which is especially hard, ordinal models had a finer prediction for early disease stages, when disease evolution is slower. Our results highlight the need to rely on multi-variate data to increase prediction performance at later disease stages, when rating scores changes are of higher magnitude.

Future applications could use this model as a tool to extract information from specific items after a multivariate analysis. We believe this can be the first step to building new composite scores such as the PACC5 in Alzheimer's disease<fnref n="36" />. One idea would be to use the Fisher information as the metric of choice to select meaningful items in a composite score, since in IRT it is common to compute Fisher information of items for the individual parameters.

On the downside we see the limitations of the ordinal method when used on data with few values for certain levels (in our case these were the high values of the MDS UPDRS II), as the added parameters of the ordinal DCM require more data in order to be properly learned. These extra parameters also come at a computational cost within the estimation algorithm, therefore we deem the method useful when the items are rated on short scales like the Likert one.

## 2.4 Parallel with Item Response Theory

As mentioned in the introduction of this chapter, the study of discrete data is the focus of the Item Response Theory (IRT) field. This has been developed mostly in psychology<fnref n="37" />, and was aimed at understanding how to best design tests to evaluate a latent trait of the patient. It makes sense that the same theory applies to cognitive tests. The other field of application of IRT has been the design of educational tests. In this context, the latent trait refers to the students' level of knowledge, which is unknown but needs to be evaluated by the test designer. By utilizing the prior knowledge of the estimated mean latent trait, the test difficulty can be tailored to minimize variance in the estimation of the students' latent trait. This theory has found extensive applications in binary and ordinal tests. The models extend to several items for multivariate tests, and can also include a multidimensional latent trait for more complex application cases. The IRT has also been adapted to longitudinal cases<fnref n="38" /> and recently new applications to Parkinson's disease have emerged<fnref n="39" />, but

<footnote n="36">³⁶K. V. Papp, D. M. Rentz, I. Orlovsky, et al., "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.</footnote>

<footnote n="37">³⁷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.</footnote>

<footnote n="38">³⁸C. Proust-Lima, V. Philipps, B. Perrot, et al., "Modeling repeated self-reported outcome data: A continuous-time longitudinal Item Response Theory model," en, *Methods*, vol. 204, pp. 386-395, Aug. 2022.</footnote>

<footnote n="39">³⁹S. Luo, H. Zou, G. T. Stebbins, et al., "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.</footnote>
there are very few attempts at combining both into a multidimensional longitudinal latent trait<fnref n="40" />.

Let us introduce some notations. We will keep the index *i* for the individual, *k* for the feature which is called item in this section, and *j* for the visit index in the longitudinal framework. The observations $y_{ijk}$ are the results of the answers to item *k* for individual *i* at visit *j*, and in the case of pure cross-sectional data we drop the visit index $y_{ik}$. The most basic formula for an IRT model with binary items and no longitudinal modelling is as follows:

$$ \mathbb{P}(y_{ik} = 1) = \eta_k(\theta_i) \qquad (2.4.1) $$

$$ \eta_k(\theta_i) = \frac{1}{1 + \exp(-a_k(\theta_i - b_k))} \qquad (2.4.2) $$

where $(a_k, b_k)$ are model parameters and $\theta_i$ is the latent trait of individual *i*. $b_k$ refers to the level of the test, i.e. the higher $b_k$ the lower the probability of having a right answer to the item. Ideally the designer of the test should aim for a $b_k$ close to the mean of the individual latent traits, or in the case of an educational test, the objective should be to have $b_k$ equal to estimated level required to pass the test. $a_k$ is the discriminating power of the item. If $b_k$ is required level to pass the test, then a strong discriminating power $a_k$ will imply that individuals with sufficient knowledge should almost always pass while individuals below the threshold have almost no chance. Sometimes in the case of multiple answer questionnaire, the model takes into account a part of randomness inherent to the possibility of answering correctly by chance:

$$ \mathbb{P}(y_{ik} = 1) = c_k + (1 - c_k)\eta_k(\theta_i) \qquad (2.4.3) $$

Now let's complicate the model by introducing the longitudinal model. The latent trait $\theta_i$ becomes a function of time. The most simple choice for modelling it is to have a linear mixed-effect structure for the time:

$$ \theta_i(t) = \alpha_i t + \beta_i \qquad (2.4.4) $$

with the time slope $\alpha_i \sim N(\bar{\alpha}, \sigma_{\alpha}^2)$ and the intercept $\beta_i \sim N(\bar{\beta}, \sigma_{\beta}^2)$. Then if we add a multi-dimensional structure for $\theta_i$ we obtain: $\theta_i \in \mathbb{R}^{N_s}$, with $N_s$ the dimension of the latent trait, and the formula for each item now includes a scalar product $d_k^T \theta_i$ instead of $\theta_i$. All in all we have the formula below:

$$ \mathbb{P}(y_{ijk} = 1) = \frac{1}{1 + \exp(-a_k(d_k^T(\alpha_i t_{ij} + \beta_i) - b_k))} \qquad (2.4.5) $$

which is to be compared with the binary DCM model formula:

$$ \mathbb{P}(y_{ijk} = 1) = \frac{1}{1 + \exp\left(\tilde{g}_k - \frac{v_k(e^{\xi_i(t_{ij} - t_0 - \tau_i) + t_0}) + w_{ik}}{p_k(1 - p_k)}\right)} $$

<footnote n="40">⁴⁰J. M. te Marcelde, C. A. W. Glas, G. Van Landeghem, et al., “Application of Multidimensional Item Response Theory Models to Longitudinal Data,” en, *Educational and Psychological Measurement*, vol. 66, no. 1, pp. 5–34, Feb. 2006.</footnote>
We observe that even though there are subtleties in the formulas, with different identifiability conditions and distribution priors over the random effects, some parameters have essentially the same role. For instance, we see that the DCM model only has a one-dimensional acceleration, so it can be directly compared to a specific case of the latent trait formula where $\alpha_i = \lambda_i \bar{\alpha}$ with $\lambda_i$ being an uni-dimensional factor. This boils down to the hypothesis of parallel trajectories in the DCM that does not allow a patient's acceleration to vary depending on the item $k$. In this case we see that $e^{\xi_i}$ and $\lambda_i$ play the same role, while $a_k d_k^T \bar{\alpha} \approx \frac{v_k}{p_k(1-p_k)}$ reflects the discriminating power of the item, $b_k \approx \tilde{g}_k p_k (1-p_k)$ is the difficulty of the item, and finally the space-shift can be decomposed over the sources with the $k$-th column of the mixing matrix $A$ as $\frac{A_k s_i - v_k e^{\xi_i} \tau_i}{p_k(1-p_k)} \approx a_k d_k \beta_i$. The last formula shows that in multidimensional IRT $\beta_i$ plays a similar role to the sources in the DCM model, and we could also enforce the ICA as in the DCM by specifying the following prior: $\beta_i \sim N(\bar{\beta}, \sigma_{\bar{\beta}}^2 I_{N_s})$.

In the ordinal case, the comparison with Samejima's model based on the cumulative logit$^{41}$ is ever so simpler:

$$P(y_{ijk} \geq l) = \eta_{kl}(\theta_i(t_{ij})) \qquad (2.4.6)$$

$$\eta_{kl}(\theta_i) = \frac{1}{1 + \exp(-a_k(\theta_i - b_{kl}))} \qquad (2.4.7)$$

where the difficulty of the item $b_{kl}$ now depends on the level $l$ of the item, with increasing difficulty as levels increase. These parameters can be directly compared to the $\sum_{m=1}^l \delta_k^m$, which encodes this increasing difficulty of the item in the ordinal DCM model.

Overall the discrete versions of the DCM model are very similar to IRT models. The specificity lies in the hypotheses of the DCM model, namely the uni-dimensional latent acceleration random effect and the independence prior of the sources ($\approx \beta_i$), constituting the space-shift ($\approx$ intercept).

### 2.4.1 Fisher information

Computing Fisher information in IRT is usually the way to estimate the variance of the parameters$^{42}$ and also a way to understand which values in the ordinal scale are bringing the most information. We transpose these ideas to the longitudinal case. For the logistic curves model, we computed explicitly the Fisher information formula for both the binary and the ordinal models. These formulas can be found in the appendix 15. In IRT, computing the Fisher information for item $k$ with the standard model 2.4 gives:

$$\mathcal{I}_k(\theta) = \mathbb{E}_{X \sim f(\cdot; \theta)} \left( \left( \left( \frac{\partial}{\partial \theta} \log f(X; \theta) \right)^2 \right) \middle| \theta \right) \qquad (2.4.8)$$

$$= a_k^2 \eta_k(\theta) (1 - \eta_k(\theta)) \qquad (2.4.9)$$

The main use case of this formula relies on the fact that in the asymptotic case we can have the approximation $V(\hat{\theta}) = \frac{1}{\bar{I}(\theta)}$ where $\hat{\theta}$ denotes the maximum likelihood estimation (MLE) of $\theta$.

<footnote>$^{41}$F. Samejima, “Estimation of latent ability using a response pattern of graded scores,” en, *Psychometrika*, vol. 34, no. 1, pp. 1–97, Mar. 1969.</footnote>

<footnote>$^{42}$A. Ly, M. Marsman, J. Verhagen, et al., “A Tutorial on Fisher information,” en, *Journal of Mathematical Psychology*, vol. 80, pp. 40–55, Oct. 2017.</footnote>
Indeed, if we consider that data was generated according to the model with a true parameter $\theta^*$, then the MLE will converge towards the true parameter $\theta^*$ as the number of samples $n$ goes to infinity according to the following formula:

$$ \sqrt{n}(\hat{\theta} - \theta^*) \xrightarrow[n \to \infty]{D} \mathcal{N}(0, \mathcal{I}(\theta^*)^{-1}) \qquad (2.4.10) $$

where $D$ means that the random variable on the left converges in distribution towards the distribution on the right. This formula still holds for the multidimensional $\theta$ case, where $\mathcal{I}(\theta)$ is the Fisher information matrix.

Of course, since our model uses a Bayesian framework, computing confidence intervals is not an issue as it will be a byproduct of our posterior estimation for latent parameters. The formula can still be applied to the model parameters which are maximized, and therefore do not have a posterior distribution. However we usually are more interested in individual parameters, as they allow us to personalize the model and run various analysis, including correlations with covariates, clustering, prediction, patient's stratification... For those parameters, the use of Fisher information as a mean to compute a confidence interval is of no use after the calibration by the MCMC-SAEM algorithm. The only potential application is for the personalization step: when we add new patients and try to compute their optimal individual trajectory by maximizing the likelihood over the individual parameters, we obtain only point estimates for those parameters and therefore the Fisher information formula for confidence intervals comes in handy.

But Fisher information can also be of interest in order perform item selection$^{43}$. In our case, it could be used to evaluate the influence of items over time: we can for instance use it as a criterion to discriminate "early" versus "late" items, as well as understanding which items are most informative for each specific parameter of the model. In the case of the ordinal model we can compute the amount of information per level of each item. As the information is computed as a function of the time, we obtain profiles describing the informativeness of items at different stages of the disease.

Another interesting application is the use of this information over time to build new composite scores with a focus on specific aspects of the disease. For instance if we need to spot early motor patients in Parkinson's disease, we can look into which items provide the most early information for motor-related model parameters. We provide an example below for the Parkinson's disease ordinal model presented in the previous section.

As we see in the formula 0.58 or in the formula 2.4.8, the Fisher information of the chosen parameter depends on the model value for this parameter. However the model is function of time, so the Fisher information is varying over time. If we included a prior over the time of the visits $t_{ij} \sim D$ then we could have integrated the information over time, but this is less interesting for us since we can leverage this time dependence to extract crucial insight from the items.

To illustrate our assertion, we selected one parameter of interest: $\tau_i$ the time-shift, which is a proxy for the onset age of the disease. Figure 2.12 shows the evolution in time of the Fisher information for $\tau$ for four selected items. We computed these curves as the mean over all the individual trajectories parametrized by $z_i$. The patterns for the different items are very different, for instance we see that item 8 of part I has very low information, but also seems to be less informative as time passes. On the contrary items such as MDS-UPDRS II-1 or III-1 bring more

<footnote>$^{43}$R. J. Swartz and S. W. Choi, "A Burdened CAT: Incorporating Response Burden With Maximum Fisher's Information for Item Selection," en,</footnote>
<figure><img src="image_28.png" /><figcaption>Figure 2.12: Tau information as a function of time for four selected items. The curve was obtained as the mean information over all individuals.</figcaption></figure>

information as the disease progresses. Other items such as MDS-UPDRS II-10 have fluctuating information.

The information profiles can be used to design specific sub-scores targeting early disease stage for instance. We can select items that provide enough information in the early stages to form a new composite score such that the estimation of the individual parameters in a DCM model is as precise as possible.

As mentioned previously, the Fisher information can be used to approximate the variance of individual parameters estimation during personalization. We used this approximation to compute the standard deviation of $\tau$ in Figure 2.13 depending on the number of visits used for personalization. This is very useful in order to understand how reliable the estimates obtained by personalization are, and therefore whether the predictions we extrapolate based on these individualized trajectories are precise. The violin plot shows that with increasing number of visits the estimation gets better.

More importantly we observe that with less than three visits the estimated variance of $\tau_i$ is quite high, which means that the model requires multiple visits to perform well. This is paramount in the context of applications to clinical trials where a disease progression model can provide a prognosis score to be included in a PROCOVA<sup>44</sup> for instance, or perform patient selection based on some criterion such as $\xi_i > 1$ to only retain fast-progressors<sup>45</sup>. Other applications of cognitive decline predictions have been shown to allow clinical trial enrichment<sup>46</sup>. In this context we see that using only a baseline visit to perform the personalization is rather limiting compared to using a few follow-ups visits. The question that hinders this larger visit sampling is whether the pay-off for following patients before starting the trial is worth it given that medical examinations might be costly and that there is a very limited legal time-window between patient enrollment and clinical trial debut. The optimal application case however is when patients are enrolled from an observational cohort, where they tend to have several years of prior follow-up.

Finally the Fisher information can be leveraged to improve the design of the test<sup>47</sup>. We computed

<footnote n="44"><sup>44</sup>D. Bertolini, K. Arnemann, D. Hall, et al., "Machine Learning Enables Smaller ALS Clinical Trials (P1-13.003)," en, *Neurology*, vol. 98, no. 18 Supplement, May 2022.</footnote>

<footnote n="45"><sup>45</sup>E. Maheux, J. Ortholand, C. Birkenbihl, et al., "Forecast Alzheimer's disease progression to better select patients for clinical trials," en, Jul. 2021.</footnote>

<footnote n="46"><sup>46</sup>A. Tam, C. Laurent, S. Gauthier, et al., "Prediction of Cognitive Decline for Enrichment of Alzheimer's Disease Clinical Trials," *The Journal Of Prevention of Alzheimer's Disease*, 2022.</footnote>

<footnote n="47"><sup>47</sup>Y. Jung and I. Lee, "Optimal design of experiments for optimization-based model calibration using Fisher</footnote>
<figure><img src="image_29.png" /><figcaption>Figure 2.13: Violin plot of the approximated standard deviation of $\tau_i$ depending on the number of visits used for personalization.</figcaption></figure>

the information for each individual parameter, for each item, level per level. The results are shown in Figure 2.14. The total information is obtained by summing over all visits, rather than averaging, so that it also reflects the imbalance of item level representation. The goal is to spot whether some items levels are not useful and whether the item should be reworked. We see that some items prove to be more useful for $\tau$ such as MDS-UPDRS III-2 or III-14, whereas some items bring little to nothing such as I-6 or III-11. Overall the information is mainly dominated by item level 1 with a bit level 2. Indeed the MDS-UPDRS scale is known to lean towards late disease stages sensitivity to the detriment of early sensitivity<sup>48</sup>, which is one of its main drawbacks when used as an endpoint for clinical trials which are targeting early stages. Items with level 3 and 4 information are very seldom, meaning that the item scales are not really adapted to the distribution of patients in this cohort. However PPMI is not skewed towards early PD, so this indicates that the MDS-UPDRS scale suffers from an uninformativeness of the higher level values of its items, which are probably too severe to be observed.

information matrix,” en, *Reliability Engineering & System Safety*, vol. 216, p. 107968, Dec. 2021.

<sup>48</sup>M. H. Tosin, G. T. Stebbins, C. Comella, et al., “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.
<figure><img src="image_30.png" /><figcaption>Figure 2.14: Fisher information for $\tau_i$ per item, per level. The information was summed over all visits, therefore there is an imbalance because of the item levels not being equally represented.</figcaption></figure>
## 2.5 Conclusion

This chapter was centered around the extension of the DCM model to discrete data, i.e. binary and ordinal data. These types of data are ubiquitous in neurodegenerative cohorts where symptoms and cognitive scales are at the core of disease progression. We showed how the application of the DCM model to this type of data enlarged the horizon of potential analysis. From symptom onset prediction to fine-grained analysis of a cognitive scale, we showcased a wide array of interesting applications.

We also encountered the limitations due to modelling hypotheses and observations distribution in the cohort. Symptoms might be treated, which is in contradiction with the model's mono-tonic increase of the logistic. In the ordinal case, estimation of the model is very sensitive to the representation of the levels in the cohort and requires much more computational time.

The analysis of discrete data paved the way for our next chapter. Indeed, due to the nature of binary and ordinal observations, patients tend to be easily separable into groups, for instance those with a specific symptom and those without. Therefore we observed a lot of disparity in the cohorts based on those discrete features, which reflected in the estimated individual parameters of the DCM models that we estimated. The simulated data in the synthetic experiment section was purposefully designed to have the clusters in Figure 2.5. This illustration is an exaggeration of real-life data, but it shows that we now have the tools to uncover new clusters of patients. Disparities that were previously merged in aggregated scores can be disentangled with an item-level modelling. Techniques to extract those clusters are the focus of the next chapter of this thesis.
# Chapter 3

## Subtyping

This chapter addresses a fundamental issue of disease progression modelling. Theoretical hypotheses underlying the use of mixed-effect models with unimodal priors apply when the dataset is an homogeneous population of individuals. However neurodegenerative diseases most often describe a spectrum, including more specific diseases, for instance Parkinson's disease and Lewy-body disease<fnref n="1" />. We therefore introduced the mixture models for Disease Course Mapping (DCM), allowing to take disease subtypes into account.

The method presented in this chapter led to the paper<fnref n="2" />. My contributions include the extension of the model to include latent classes, declined in two possible variants. The main novelty in this work however lies in a new heuristic to improve the estimation of the mixture, which is known to be the challenging part. This heuristic is described in the estimation chapter 5.2.1. The chapter is composed of two main parts, describing the two variants of the mixture models for DCM.

### Contents of the chapter

<table><tr><td>3.1 Context</td><td>86</td></tr><tr><td>3.2 Mixture of DCM models</td><td>88</td></tr><tr><td>3.3 Mixture of individual parameters</td><td>89</td></tr><tr><td>3.4 Experiments</td><td>90</td></tr><tr><td>3.4.1 Simulated data: Univariate model</td><td>90</td></tr><tr><td>3.4.2 Applications to Alzheimer's disease data</td><td>91</td></tr><tr><td>3.4.3 Experiment with mixture of individual parameters</td><td>95</td></tr><tr><td>3.4.4 Combination with DCM models on discrete data</td><td>97</td></tr><tr><td>3.5 Discussion</td><td>99</td></tr><tr><td>3.6 Conclusion</td><td>99</td></tr></table>

<footnote n="1">¹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.</footnote>

<footnote n="2">²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, et al., Eds., ser. Lecture Notes in Computer Science, Cham: Springer International Publishing, 2021, pp. 571–582.</footnote>
## 3.1 Context

Large medical cohorts with numerous biomarkers and cognitive scores are particularly well suited to analyse disease at a population level. Understanding how the disease progresses and what is the expected variability amongst individual is key. This variability can be driven by certain covariates, as showcased in Alzheimer's disease (AD). Indeed the genetic factor APOE-ε4 has major influence over the speed of the disease, and sex is also implicated in neurodegeneration, as shown in<fnref n="3" />. But when no prior knowledge and no specific covariate are available, the variability of individual disease patterns which is learned by the model is constrained by the specifications of the said model.

For the DCM model, as well as for standard mixed-effect models, the main limitation is the hypothesis of a unimodal distribution for random effects. The gaussian prior used in the DCM implicitly describes the population as homogeneous, and regularizes around it. Another constraint on the variability specific to the DCM model lies in the uniform speed hypothesis. Indeed the only random effect influencing the speed of the disease is the acceleration factor $\alpha_i$, and it accelerates the disease timeline $\psi_i$ which then accelerates uniformly all the markers progression $\eta_k$. This choice of modelling does not allow a patient to have a cognitive decline twice as fast as the normal progression while showing standard MRI brain atrophy for instance. Figure 3.1 illustrates how the value of $v$ shapes the possible trajectories in the manifold $(0, 1)^2$.

These limitations are not adapted to the reality of neurodegenerative diseases, where different subtypes are often observed. For instance in AD, biological subtypes have been identified based on the tau pathology and brain atrophy<fnref n="4" />. In Parkinson's disease (PD), there is a commonly used classification between tremor dominant (TD) patients and postural gait instability disorder (PIGD) patients<fnref n="5" />. However these two groups do not clearly separate the whole population, and the question remains whether these subtypes are not actually two extremes of a continuum of variability. Moreover recent studies have shown that the subtyping is not stable in time<fnref n="6" />, thereby questioning the classification system. We will discuss this particular case in the experiments section.

The heterogeneity in disease presentation is therefore structured around the subtype in which the individual belongs. Thus, subtyping is an essential task in disease progression modelling. However uncovering unknown clusters of individuals is a complex challenge on top of the estimation of the model parameters. It might thus be tempting to separate the clustering and the model estimation, for instance by performing a *posteriori* clustering. However this naive approach does not allow an unbiased estimation of model parameters. We will show with the DCM model to what extent ignoring subtypes in the model estimation leads to estimation errors.

The problem of subtyping has been well studied and we will focus on the literature of disease progression models. For mixed-effect models subtyping implies multimodal distributions, leading to mixture models. The mixture can be introduced at different levels. The most common is the mixture of random effects, as is implemented in longitudinal mixed-effect statistical frameworks

<footnote n="3">³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.</footnote>

<footnote n="4">⁴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.</footnote>

<footnote n="5">⁵G. T. Stebbins, C. G. Goetz, D. J. Burn, et al., “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.</footnote>

<footnote n="6">⁶T. Simuni, C. Caspell-Garcia, C. Coffey, et al., “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.</footnote>
<figure><img src="image_31.png" /><figcaption>Figure 3.1: Disease Course Mapping geometric variability depending on $v$. The bundle of trajectories is obtained by varying the values of the space shift $w_i$ in a logistic curves DCM model (variability due to temporal random effects is not seen in the spatial representation). Top plot: $v = (e^{-1}, e^{-1.5})$. Bottom plot: $v = (e^{-1.5}, e^{-1})$.</figcaption></figure>
such as latent classes mixed models (LCMM)<fnref n="7" /> or MONOLIX<fnref n="8" />. Another approach is to introduce a mixture at the residual noise level, which means that we can have several structural models describing the population subtypes. The challenge behind mixture models is the estimation, and approaches vary depending on the mixed-effect model chosen. The choices of algorithms will be further discussed in the dedicated chapter on optimization.

For disease subtyping many other solutions have been proposed outside of the scope of mixed-effects models. The Subtyping and Stage Inference algorithm (SuStaIn)<fnref n="9" /> is an extension of an event-based model with a mixture and has shown very specific patterns of disease progression in AD<fnref n="10" />. Also in the context of AD and identification of disease subtypes in the Alzheimer's Disease Neuroimaging Initiative (ADNI) cohort, *ad hoc* methods have been developed for the classification of subtypes. The most famous one is the Murray-Dickson algorithm<fnref n="11" />, and its adaptations<fnref n="12" />. More traditional models such as latent Dirichlet allocation have also been applied to the cohort<fnref n="13" />.

Our choice to deal with this heterogeneity is thus situational, as it depends on the specification of the DCM model. The main concern for the ability of the model to accurately describe subtypes lies behind the limitation shown in figure 3.1, as we hypothesize that subtypes might lead to different *geometric* patterns rather than *temporal* patterns. This is why we first constructed the mixture model at the residual noise level, meaning that each cluster of individuals has a different DCM model, thus a different *v*. In a second time we also describe a model with the mixture of random effects, which is more adapted to the identification of temporal clusters such as fast progressors versus slow progressors.

The work in this chapter presents materials from the published paper<fnref n="14" /> which focused on the mixture at the residual noise level. We extend it with additional experiments as well as a comparison with a second mixture model at the individual parameters level.

## 3.2 Mixture of DCM models

The mixture of DCM models is a new layer atop of the hierarchical structure already built. If we write $q(y | z; \theta)$ the residual noise likelihood of the model parametrized by $\theta$, with latent variables

<footnote n="7">⁷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.</footnote>

<footnote n="8">⁸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="9">⁹A. L. Young, R. V. Marinescu, N. P. Oxtoby, et al., "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.</footnote>

<footnote n="10">¹⁰C. Shand, P. J. Markiewicz, D. M. Cash, et al., "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.</footnote>

<footnote n="11">¹¹M. E. Murray, N. R. Graff-Radford, O. A. Ross, et al., "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.</footnote>

<footnote n="12">¹²S. L. Risacher, W. H. Anderson, A. Charil, et al., "Alzheimer disease brain atrophy subtypes are associated with cognition and rate of decline," en, *Neurology*, vol. 89, no. 21, pp. 2176–2186, Nov. 2017.</footnote>

<footnote n="13">¹³X. Zhang, E. C. Mormino, N. Sun, et al., "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.</footnote>

<footnote n="14">¹⁴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, et al., Eds., ser. Lecture Notes in Computer Science, Cham: Springer International Publishing, 2021, pp. 571–582.</footnote>
z and observations y, then the residual noise likelihood of a mixture model with L clusters is:

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

where $\pi^c$ denotes the probability of each cluster c with its own set of parameters $\theta^c = (\bar{\tau}^c, \bar{\xi}^c, \sigma_\tau^c, \sigma_\xi^c)$ and fixed effects $z_{pop}^c = (g^c, v^c, \beta^c)$. We assume that each individual has a true latent class $c_i$ such

that $\pi^c = \frac{\sum_{i=1}^{N} 1_{c_i=c}}{N}$ and a set of individual parameters $z_i^{c_i} = (\tau_i^{c_i}, \xi_i^{c_i}, s_i^{c_i})$ corresponding to the model attached to the cluster $c_i$. The total log-likelihood of the mixture model then writes:

$$\log Q(\mathbf{y}, \mathbf{z}, \theta^1, \dots, \theta^L) = \sum_{i=1}^{N} \log q(\mathbf{y} | \mathbf{z}_i^{c_i}, \mathbf{z}_{pop}^{c_i}; \theta^{c_i}) \qquad (3.2.2)$$

$$+ \sum_{i=1}^{N} \log q(\mathbf{z}_i^{c_i} | \mathbf{z}_{pop}^{c_i}; \theta^{c_i}) \qquad (3.2.3)$$

$$+ \sum_{c=1}^{L} \log q(\mathbf{z}_{pop}^{c}; \theta^{c}) \qquad (3.2.4)$$

The estimation algorithm is described in the optimization chapter. It is a mixture version of the Monte Carlo Markov Chain Stochastic Approximation Expectation Maximization (mMCMC-SAEM) algorithm, where we add the latent variables $\pi_i^c$ for the probability of each individual i to be in each cluster c, and $(z_i^c)_{i}$ are the individual parameters for individual i in cluster c. This allows for soft clustering throughout the estimation, avoiding one of the pitfalls of optimization with mixtures where sometimes individuals are assigned to a given cluster with hard labels.

## 3.3 Mixture of individual parameters

The second approach is to add the mixture at the individual parameters level. We still keep the notations for $\pi^c$ the probability of cluster c and $c_i$ the true class of individual i. The individual parameters priors become:

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

$$q(s_{ik} | \theta) = \sum_{c=1}^{L} 1_{c_i=c} \phi(s_{ik}; \bar{s}_k^c, \sigma_{s_k}^c) \qquad (3.3.2)$$

where $\phi$ is the probability density function of a normal distribution (multivariate in the case of $\tau, \xi$), and $(\bar{\tau}^c, \bar{\xi}^c, \Sigma_{\tau, \xi}^c)$ are model parameters corresponding to the mean and covariance of the $(\tau_i, \xi_i)$ in cluster c. We also introduce a mean and a standard deviation for the sources. Parameters $(\beta, g, v)$ are shared by the clusters.

**Remark:** We removed the independence of individual parameters in the priors, by adding a covariance matrix for temporal parameters, while still keeping the sources independent (for the
Independent Component Analysis on the space-shifts to be relevant). In the experiment section we show why this was necessary.

Some conditions are required for the model to still be identifiable. Namely:

$$ \mathbb{E}(\xi_i|\theta) = \sum_{c=1}^{L} \pi^c \bar{\xi}^c = 0 \qquad (3.3.3) $$

$$ \mathbb{E}(\boldsymbol{s}_i|\theta) = \sum_{c=1}^{L} \pi^c \bar{\boldsymbol{s}}^c = \mathbf{0} \qquad (3.3.4) $$

$$ \mathbb{V}(s_{ik}|\theta) = \sum_{c=1}^{L} \pi^c ((\sigma_{s_k}^c)^2 + (\bar{s}_k^c)^2) = 1 \qquad (3.3.5) $$

Due to the maximization step in the EM algorithm not being explicit with the last two constraints on the sources 3.3.4 3.3.5, we only implemented a mixture for the time-related individual parameters $(\tau_i, \xi_i)$.

## 3.4 Experiments

### 3.4.1 Simulated data: Univariate model

We first show that the mixture of DCM models accurately recovers the ground truth parameters on simulated data. We generated data by initializing two unimodal disease course mapping models with chosen population parameters. We use the Kullback-Leibler divergence to determine the difference between the two model distributions. We then assume the two clusters have equal prevalence and we create 512 individuals, which consists in randomly attributing the individual to a cluster and sample individual parameters according to the distribution of the cluster. Next we arbitrarily decide the number and time of the "visits" for each individual, and we compute the values associated to these timepoints. Finally we add a small Gaussian noise to the output.

A simplified version of the model consists in using only one dimension $d = 1$, leading to an univariate model with no space-shifts. In this case, the only individual parameters left are the time-related ones: $(\xi_i)_i$ and $(\tau_i)_i$. This special case of the model converges much faster (about one minute for 2000 iterations), which allows us to initialize the model without having to use the initialization method described in the optimization section 19. We evaluated the class estimation with the area under the ROC curve (ROC AUC) metric. For the evaluation of population and individual parameters, we computed the absolute error between the ground truth parameters of each cluster and the parameters of the closest estimated cluster, this error being normalized by the standard deviation of the parameter. Figure 3.2 shows the reconstruction metrics on simulated data as a function of the Kullback-Leibler divergence between clusters.

With the ROC AUC, we are able to see that the mixture perfectly separates two clusters, once the two clusters are different enough. The estimation errors for population parameters increases slightly as the clusters get farther from each other, which is understandable since the algorithm needs to do more exploration. The mean error on individual parameters shows a very stable curve. The reconstruction is limited by the noise of the generated data. In all the univariate experiments,
<figure><img src="image_32.png" /><figcaption>Figure 3.2: Evaluation metrics for the reconstruction of ground truth parameters. Left: ROC AUC for the estimated cluster probabilities of individuals $\pi_i^c$. Right: Absolute error between estimated parameters and ground truth, normalized by parameter standard deviation.</figcaption></figure>

the estimated parameters allowed the model to reach a L2 loss close to the noise level. The consistent errors for individual parameters seemingly correspond to the range in which individual parameters differences can be mistaken for noise. For the population parameters, we see an increase in error as the clusters get farther from each other. This is likely due to the fact that the parameters need to drift farther from their initialization which hints towards an unfinished convergence. The DCM model implementation lacks a convergence criterion, a topic which will be discussed in the optimization chapter.

### 3.4.2 Applications to Alzheimer's disease data

We apply the method to the ADNI dataset. The data used in this chapter were obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu), launched in 2003 as a public-private partnership, led by Principal Investigator Michael W. Weiner, MD. The primary goal of ADNI has been to test whether serial magnetic resonance imaging (MRI), positron emission tomography (PET), other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of mild cognitive impairment (MCI) and early Alzheimer's disease (AD). The versions of the dataset used for this experiment were ADNI 1, 2, 3 and ADNI GO.

This application is a simple case to verify that our model works in practice, as we will later focus on Parkinson's disease which is the goal of our subtyping study. We first build a model considering only stable Alzheimer's disease (AD) patients and stable controls. This means that across all of their visits the stable controls only have a healthy control diagnosis, while stable AD must have at least two AD confirmed diagnosis. This leaves us with 618 controls and 730 AD patients, with an average of 6.5 visits per patient. The rest of the patients are considered as part of the Mild Cognitive Impairment (MCI) category.

We select 8 features that are most relevant to the progression of AD, which included cognitive scores, MRI-derived regional volumes and biomarkers level. Table 3.1 summarizes the features we used. MRI measures are intracranial volume measures. For MRI and CSF values, the data were
<table><caption>Table 3.1: ADNI demographics (includes controls, AD patients and MCI). MMSE: Mini-Mental State Examination, a cognitive score ranging from 0 (worst state) to 30 (perfectly healthy). ADAS-cog 13: Alzheimer’s Disease Assessment Scale cognitive evaluation is a s-cognitive score ranging from 0 (healthy) to 85 (worst). RAVLT: Rey Auditory Verbal Learning Test, scaled between 0 (best) and 1 (worst). MRI measures for hyppocampus, ventricles and entorhinal are known to be sensitive to AD progression. CSF: Cerebrospinal fluid measures for proteins p-Tau and amyloid β42, which are at the heart of the biological process of neurodegeneration in AD, namely the amyloïd cascade.</caption><thead><tr><th></th><th>MMSE</th><th>ADAS-cog 13</th><th>AVLT</th><th>MRI Hippocampus</th></tr></thead><tbody><tr><th>count</th><td>10246</td><td>9912</td><td>10131</td><td>9008</td></tr><tr><th>mean</th><td>26.726430</td><td>17.458136</td><td>0.529044</td><td>0.429449</td></tr><tr><th>std</th><td>3.931331</td><td>11.409325</td><td>0.183635</td><td>0.174020</td></tr><tr><th>min</th><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><th>25%</th><td>25</td><td>9</td><td>0.40</td><td>0.30</td></tr><tr><th>median</th><td>28</td><td>15</td><td>0.55</td><td>0.42</td></tr><tr><th>75%</th><td>29</td><td>23</td><td>0.66</td><td>0.55</td></tr><tr><th>max</th><td>30</td><td>85</td><td>1</td><td>1</td></tr><tr><th></th><th>MRI Ventricles</th><th>MRI Entorhinal</th><th>CSF p-Tau</th><th>CSF Abeta42</th></tr><tr><th>count</th><td>9000</td><td>9006</td><td>2701</td><td>2702</td></tr><tr><th>mean</th><td>0.482749</td><td>0.490103</td><td>0.446442</td><td>0.397885</td></tr><tr><th>std</th><td>0.180927</td><td>0.177736</td><td>0.185596</td><td>0.199178</td></tr><tr><th>min</th><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><th>25%</th><td>0.35</td><td>0.37</td><td>0.32</td><td>0.23</td></tr><tr><th>median</th><td>0.48</td><td>0.47</td><td>0.44</td><td>0.40</td></tr><tr><th>75%</th><td>0.60</td><td>0.60</td><td>0.57</td><td>0.53</td></tr><tr><th>max</th><td>1</td><td>1</td><td>1</td><td>1</td></tr></tbody></table>

re-scaled between 0 and 1 so that the progression can be modelled by a logistic curve. If needed the score is reversed to be increasing instead of decreasing. We kept cognitive scores on their original scale in table 3.1, but they were applied the same re-scaling before being fed to the model, as we consider them as continuous measures due to the fine-grained scale.

We first fit a single DCM model on it. The posterior analysis of individual parameters highlights the need to take heterogeneity into account. The posterior distribution of the couple of temporal parameters ($\tau_i, \xi_i$) is shown in Figure 3.3. The two modes correspond to the two classes of patients selected for this dataset. The point of this experiment is to emphasize the difference between the mixture of DCM models that we chose and a mixture of individual parameters model, which would only separate AD and controls in this case.

We then used the mixture of DCM. Optimal number of clusters was decided based on the Bayesian Information Criterion (BIC). The results of the mixture model after 4000 iterations of mMCMC-SAEM are shown in Figure 3.4. Figure A shows the probability of individuals to belong to cluster 1 vs cluster 2 (i.e. individuals on the left of the histogram belong to cluster 2). Overall cluster 1 gathers about 70% of the individuals. Notice that almost all controls belong to cluster 1. Figure B shows the geometric trajectories of the two clusters shown on two dimensions (the cognitive score MMSE and the MRI volume measure of the hippocampus). We chose these two dimensions as
<figure><img src="image_33.png" /><figcaption>Figure 3.3: Scatter plot of individual parameters, with kernel density estimation (KDE). Left: KDE on all AD and control patients. Right: KDE estimated separately for AD and controls.</figcaption></figure>

they are probably the most important AD-related features amongst the 8. We wanted to showcase
the differences in the geometric progression of biological measures against cognitive examinations.
The center of the bundle is the average population trajectory $\gamma_0$, the variability shown corresponds
to the curves obtained by varying the source 3 from -5 to +5. Notice that trajectories in cluster
2 show a decline in cognitive ability (MMSE) at a lower atrophy rate compared to cluster 1.
<figure><img src="image_34.png" /><figcaption>Figure 3.4: **A** Probability of individuals to belong to cluster 1 vs cluster 2. Left: all individuals. Right: controls only. **B** Geometric trajectories of the two clusters shown on two dimensions. **C** Cluster average trajectories in time.</figcaption></figure>
Finally figure C shows the cluster's average trajectories in time. For each cluster we plotted the average trajectory for all patients and also for AD patients as the average of cluster 1 is skewed due to most of the controls belonging to this cluster. AD patients in both clusters have a similar cognitive evolution as shown by the cognitive scores (MMSE, RAVLT and ADAS-COG) whereas AD patients of cluster 1 have a faster progression on MRI-based and CSF biomarkers. This is a result of the different geometric patterns learned (as showcased by figure B), which is allowed by the mixture of DCM models and would not have been possible with a mixture of individual parameters only.

Several studies report an association between the atrophy rates of particular brain regions and cognitive decline<sup>15,16</sup>, which lead to the identification of disease subtypes based on the differences in regional atrophies. Our analysis suggest that such associations are not systematic: similar pathological processes may lead to distinct pattern of cognitive decline, like similar cognitive decline may have distinct pathological processes.

Interestingly, the clustering does not correspond to AD and controls classes. Cluster 1 contains 88% of all controls and 65% of AD cases, meaning cluster 2 accounts for a minority of the data and contains mostly AD patients. To understand why the mixture model did not separate AD patients from controls, we estimated one model on AD patients only and one model on controls only. We then compute the log-likelihood of our models. Results are shown in table 3.2. They confirm that the mixture model explains the variability seen in the data better than prior categorisation based on diagnosis. We further estimated a mixture with the initial clusters being AD and controls to try and enforce the convergence towards a separation of AD and controls, and the algorithm still converged towards a similar version to the one presented in Figure 3.4.

<table><caption>Table 3.2: Log-likelihood of models</caption><thead><tr><td>Model</td><td>Fit (log q)</td><td>Regularization (log p)</td><td>Total log-likelihood</td></tr></thead><tbody><tr><td>Mixture</td><td>38100</td><td>-10391</td><td>27708</td></tr><tr><td>AD + Controls</td><td>37226</td><td>-9772</td><td>27454</td></tr></tbody></table>

As a final experiment, we performed individual personalization of MCI in the previous model, i.e. we estimated individual parameters including likelihood to belong to each cluster while keeping population parameters fixed. The distribution of MCI patients in the two clusters confirmed our hypothesis that cluster 2 is a specific subtype of AD: MCI associated to cluster 2 are mostly converters (66%) while 80% of non-converters are in cluster 1.

### 3.4.3 Experiment with mixture of individual parameters

We now show the difference with a mixture at the level of random effects. We first illustrate why we introduced the covariance matrix between $\tau$ and $\xi$ in the mixture 3.3.1 with a practical experiment. We used the exact same dataset as in the previous section, i.e. ADNI with stable AD and controls. As showcased by Figure 3.3, the temporal individual parameters seem to have a bi-modal distribution, inciting us to fit a model with a mixture of individual parameters. We used

<footnote n="15"><sup>15</sup>X. Zhang, E. C. Mormino, N. Sun, et al., “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.</footnote>

<footnote n="16"><sup>16</sup>S. L. Risacher, W. H. Anderson, A. Charil, et al., “Alzheimer disease brain atrophy subtypes are associated with cognition and rate of decline,” en, *Neurology*, vol. 89, no. 21, pp. 2176–2186, Nov. 2017.</footnote>
the model described by equation 3.3.1, with no mixture on the sources since the implementation
would be painstakingly difficult due to the constraints on the sources global distribution to keep the
ICA. If we remove the covariance matrix, meaning that we keep τ and ξ independent, the model
has a hard time finding two modes where the two coordinates are independent. We showcase the
standard fit of two clusters in figure 3.5.

<figure><img src="image_35.png" /><figcaption>Figure 3.5: Distribution of individual parameters during the estimation (with mMCMC-SAEM). The two ellipsis correspond to the two Gaussians learned in the mixture. Each point is an individual, colored by its probability to belong to the green cluster. The red cluster progressively shrinks and eventually reaches zero probability.</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>Figure 3.6: Distribution of time-related individual parameters. Left: individual parameters ($\tau_i, \xi_i$) colored by their probability of belonging to cluster 1. Ellipsis correspond to the covariance matrix of each cluster centered around its mean. Right: individual parameters colored by their true label (AD or control), with KDE estimation.</figcaption></figure>

It might be understandable that fitting two Gaussians with independent coordinates in the mixture is prone to fail, so we added the covariance matrix for more flexibility. With this model we obtain the distribution in Figure 3.6. First there is a clear distinction between this figure and Figure 3.3. The most striking part is that here we see that $\xi$ and $\tau$ are far from independent. More specifically they tend to be anticorrelated, i.e. early onset also means faster progression. Then we notice that the two identified clusters match the KDE estimations of AD and controls, which is what was expected from the mixture of individual parameters. Finally we observe obvious facts about the clusters characteristics: AD patients (red cluster in both plots) have a lower $\tau$ and higher $\xi$, meaning that they start earlier and progress faster than controls.

### 3.4.4 Combination with DCM models on discrete data

Our goal was to be able to identify subtypes for the Parkinson's disease (PD) progression, but most of the meaningful features in PD are non-continuous variables (cognitive scores). To this end, one can combine the mixture models with the Bernoulli and ordinal noise models described in the previous chapter. For the DCM mixture model, we highlighted that one of the main perks of this mixture is to have different $v$ in the clusters of the mixture. However in Bernoulli models, which are typically used for symptoms occurrence, the model parameter $v$ is not as important as for continuous models. This due to the fact that in neurodegenerative disease once a symptom starts to appear, it is almost always present until it is treated. This implies that transitions from 0 to 1 tend to be sharp, and $v_k$ becomes the parameter determining the time window of appearance of the symptom. Therefore trying to fit mixture DCM models with Bernoulli noise is not working well. On the other hand, mixture on individual parameters is a better match for Bernoulli noise models.

For ordinal noise models, the addition of the $\delta$ parameters as population parameters increases
the potential variety of models obtained with a mixture DCM models. However the increase in
parameters due to the ordinal noise multiplied by the number of clusters makes the estimation
almost impossible for large models, i.e. with *d* large (e.g. ≥ 10) or large number of levels per item.
We estimated such as mixture on the PPMI ordinal data used to conduct experiments in chapter
2.3.3.

The fit ended (after about two weeks of estimation) with two clusters with probabilities 36% and 64%. The overall negative log-likelihood of the ordinal DCM model (without mixture) was 5.0 × 10⁵ whereas the mixture with two clusters reached 4.8 × 10⁵. This improvement is also validated but not really significant when using the Bayesian Information Criterion. We computed the BIC adpated to mixed-effects models<fnref n="17" />, i.e. *BIC* = −2 log *q* + *n*<sub>*R*</sub> log(*N*) + *n*<sub>*F*</sub> log(*N*<sub>*tot*</sub>) where *q* is the likelihood, *n*<sub>*R*</sub> the dimension of random parameters and *n*<sub>*F*</sub> the dimension of fixed parameters, *N* the number of individuals and *N*<sub>*tot*</sub> the total number of observations. The Figure 3.7 shows the mean trajectory in each cluster on the 15 first items of the MDS-UPDRS score. We see that cluster 1 (in blue, accounts for 36%) is mostly in advance of cluster 2, but the trajectories are not very different. The clustering does not correlate with known PD subtypes such as Tremor Dominant (TD) and Postural Instability and Gait Disorder (PIGD). It still remains to be proven that Parkinson’s disease is composed of distinct subtypes or whether they are only extremes of a continuum of disease trajectories.

<figure><img src="image_37.png" /><figcaption>Figure 3.7: Maximum likelihood of the average trajectory in each cluster of the ordinal DCM mixture model. In blue, cluster 1. in orange, cluster 2. MDSX_Y is the item Y of part X of the MDS-UPDRS, described in section 2.3.3.</figcaption></figure>

<footnote n="17">¹⁷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.</footnote>
**3.5 Discussion**

The experiments here show the interest of the mixture of DCM models in deciphering different *geometric* patterns, but this is often not adapted to the obvious clusters such as controls vs patients. For this a mixture of individual parameters is more indicated. The main issue with mixture models lies in the estimation, which is tackled in the dedicated section 5. This is also where the main limitations of the method are highlighted, since the convergence is very complicated in practice despite the theoretical guarantees.

When it comes to potential applications, the two mixture possibilities have very different prospects. For the mixture of individual parameters:

* The possibility of identifying new subtypes is limited, since the clusters are designed on temporal markers only. This means that subtypes are characterized by their differences in onset age and speed of progression. These two features of disease progression do not require such a complex model as a DCM to be extracted from a longitudinal cohort, and clustering can be inferred without burdening the estimation of the mixed-effect model. Therefore this mixture is more adapted if subtypes are already known to affect the temporality of the disease. In this case the mixture is used to improve the model parameter estimation and to confirm the disease subtypes

* The model scales way better with the number of clusters than the mixture of DCM models, making it easier to fit in more complex settings than just a handful of clusters

In the case of the mixture of DCM models:

* This version is more indicated when looking for subtypes with specific disease patterns. This is very useful when trying to find new trends in a longitudinal dataset that has often been through plenty of clustering methods, as the clusters found are not based on easily observed features such as speed of progression. The different clusters found rely on differences in geo-metric trajectories, therefore separating different disease courses or pathogenicity mechanisms.

* However it will not separate temporal clusters or simpler variability patterns that can be handled by individual parameters (acceleration, time-shift and space-shift). The results are thus more "cryptic".

* The mixture has a very hard time scaling with the number of clusters as the number of parameters grows linearly and estimation is already very complicated with only 2 or 3 clusters

The code and implementation of these two mixture models has been done with the open-source library Leaspy https://gitlab.com/icm-institute/aramislab/leaspy, and will be released in the near future.

**3.6 Conclusion**

Neurodegenerative diseases are known to be very heterogeneous, with subtypes having very distinct disease manifestations. One of the challenges is to find methods able to separate the individuals in classes corresponding to *biologically* meaningful subtypes. One pitfall among many is to mistake a
continuum of variability for clusters, leading to characterization of extreme patients when the "av-
erage" patient does not actually belong to any cluster. In mixed-effect model terms, we should not
try to enforce clusters with mixture when a random effect already captures the needed variability.

The need to subtype and find biologically meaningful markers of variability amongst patients
is fueled by the hope of treatment. Pharmaceutical companies and laboratories have faced many
failures in trying to develop neurodegenerative cures, especially for Alzheimer's disease. Their
interest is going now towards more personalized approaches, by trying to target a specific population
("find the right treatment for the right person at the right time" - the mantra of personalized
medicine). Hence the need to cluster patients. Treatment development are motivating most of the
research in neurodegenerative diseases, and this conveniently allows us to transition to the next
chapter of this thesis centered around treatment modelling.
<table><tr><td>4.1.2</td><td>Disease-modifying treatment effect</td><td>105</td></tr><tr><td>4.2</td><td>Additive treatment effect model</td><td>106</td></tr><tr><td>4.2.1</td><td>Latent disease progression model</td><td>109</td></tr><tr><td>4.2.2</td><td>Treatment effect model</td><td>109</td></tr><tr><td>4.2.3</td><td>Estimation method</td><td>112</td></tr><tr><td>4.2.4</td><td>Experiments</td><td>113</td></tr><tr><td>4.2.5</td><td>Discussion</td><td>128</td></tr><tr><td>4.3</td><td>Piecewise-geodesic model</td><td>129</td></tr><tr><td>4.3.1</td><td>Method</td><td>129</td></tr><tr><td>4.3.2</td><td>Experiments</td><td>132</td></tr><tr><td>4.3.3</td><td>Discussion</td><td>134</td></tr><tr><td>4.4</td><td>Conclusion</td><td>137</td></tr></table>
4.1 Introduction to treatment modelling

Modelling long term evolution of a process such as a neurodegenerative disease is a complex problem,
which is the stem challenge of disease progression models. Different approaches have been developed
to tackle this issue, either using mixed-effect models, event-based models or dynamic models based
on ODE. All these models have proven to be quite potent proponents in the task of disease history
description. However their potential in generalisation tasks such as prediction often relates to the
(explicit or implicit) priors : the "template" family of curves for mixed-effect models or the class
of functions modelling the dynamics for the ODE.

These approaches all rely on the fact that the disease progression is a common pattern, shared
amongst the population. This is quite convenient in some observational studies, but this doesn't
match the clinical trial framework. Because the treatment impacts the disease progression directly
or through the mitigation of symptoms and other manifestations, the observations are impacted
in a way that causes a rupture with the previous models : the future observations are not a
simple "mechanical" consequence of the current state of the patient, but they take into account
the treatment effect as well. Even though disease progression models perform very well when
describing the disease's evolution without external interference, which we will call the natural disease
progression, they do not adequately deal with a treatment acting over the disease progression. A
successful curative treatment can even revert the progression or decrease some biomarkers' observed
levels, making it difficult to use smooth models. Thus very few attempts have been made to include
treatment in longitudinal models describing disease progression<fnref n="3" />.

Two kinds of treatment are commonly acknowledged : symptomatic treatments and disease-
modifying ones. Symptomatic treatments only affect the consequences of the disease by alleviating
the burden of symptoms (hence the name). However they do not cure the root of the dysfunction-
ment, whereas disease-modifying drugs impact the disease much earlier in the pathology pathways.
These last ones are thus more desirable, provided that it is not too late to cure the patient (or even
better, that the neurodegenerative process can be reverted). Those treatment effects are depicted
in Figure 4.1.

The methodology to identify if a treatment (drug or therapy) is disease-modifying in clinical
trials is debated, but it stands out that a simple placebo arm vs treatment arm is not sufficient.
One has to either add delayed treatment start, or use "staggered withdrawal" strategies, in order to
determine whether patients who followed treatment during a certain period of time catch up with
placebo patients in disease state by the end of the study.

The question now is about modelling, or more especially what should we include in the models
to allow treatment modification on the outcomes, whether it be symptomatic or disease-modifying.
We will now present several ideas, within different models and contexts.

4.1.1 Modelling treatment effect

Traditionally a treatment is evaluated through a randomized controlled trial, which often has a
small size. On the other hand, one can exploit the high amount of data in observational cohorts
to better understand the treatment effect on disease progression<fnref n="4" />. One of the known issues when

<footnote n="3">³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.</footnote>

<footnote n="4">⁴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.</footnote>
<figure><img src="image_38.png" /><figcaption>Figure 4.1: Scheme of treatment effect possible types. The Y-axis represents an observable marker for the cognition (could be another observable disease progression marker). The blue curve represents the natural progression of the disease without treatment. Three types of treatment effect are shown, where a cure is a specific disease-modifying treatment. The main difference to be noted is how the progression is modified after stopping the treatment.</figcaption></figure>

using observational studies is that the data is not subject to randomization<fnref n="5" />. Furthermore, when a treatment has already been proven efficient in a disease, observational cohorts are only composed of treated subjects as everyone is prescribed the treatment. Our application focuses on Parkinson's disease, for which dopaminergic therapy is a known symptomatic treatment<fnref n="6" />,<fnref n="7" />. In the event that a cohort is only composed of patients without controls or placebo, we still would like to assess the treatment influence over the disease progression. This problem has never been addressed to our knowledge yet we believe it has many implications. We aim to learn meaningful information about the treatment without controls using a disease progression model.

But how do we model treatment effect traditionally? This field is well documented in bio-
statistics<fnref n="8" />, however, this is often done outside of the longitudinal scope. Treatment effect models
often apply to clinical trials, where the treatment effect is assessed as the mean difference in the
designated endpoint of the trial between the placebo arm and the treated arm. This difference
is known as the average treatment effect (ATE), and is measured at a precise time-point, e.g. a
fixed number of weeks from trial start. To add some useful information about variability in treat-
ment effect, one can assess the treatment effect within certain groups defined by a covariate. This
quantity is known as conditional average treatment effect. The problem of measuring treatment

<footnote n="5">⁵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.</footnote>

<footnote n="6">⁶R. M. A. de Bie, C. E. Clarke, A. J. Espay, et al., “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.</footnote>

<footnote n="7">⁷C. Lungu, J. M. Cedarbaum, T. M. Dawson, et al., “Seeking progress in disease modification in Parkinson disease,” en, *Parkinsonism & Related Disorders*, vol. 90, pp. 134–141, Sep. 2021.</footnote>

<footnote n="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.</footnote>
effect conditionally on some covariates is a standard regression problem$^9$, for which the current state-of-the-art is Bayesian Additive Regression Trees (BART)$^{10}$. This method provides a flexible and agnostic framework without requiring too much data. However the longitudinal aspect of the treatment effect is often not accounted for, and treatment effect models with progression are still not standard$^{11}$. One can still add time as a covariate for the CATE regression model, but this would leave the regression model to learn both the disease progression dynamics as well as the treatment effect. This is not efficient as we could use disease progression models to alleviate the learning task. One could add the treatment in a mixed-effect model as time-dependent covariate, but this would only estimate the treatment effect without its interaction with other covariates. Having to specify the potential interactions of the treatment with other variables is possible, but tedious and will introduce more complexity in the model.

We will show in the additive treatment model section 4.2 how we can do this. In this section, we propose a two-step approach to address these challenges by combining disease progression models with treatment effects. Our method is fully Bayesian. A longitudinal model describes the latent disease progression and a regression model learns the treatment effect as the difference between the latent disease state and the observed disease state. We will tailor our method for an application to Parkinson’s disease, but it can be generalized fairly straightforwardly. Based on the Parkinson’s Progression Markers Initiative (PPMI) observational cohort$^{12}$, we will use the model to describe the dopaminergic therapy’s effect on clinical assessments, including the motor evaluation, quality of life and imaging data.

4.1.2 Disease-modifying treatment effect

Another drawback of ATE and CATE models is that they often ignore the difference between
disease-modifying and symptomatic effects. As they are defined as a mean difference between a
treated arm and a placebo arm, they do not account for the longitudinal aspect. As noted in<fnref n="13" />,
the disease-modifying treatment effect is poorly accounted for in most models. Often it is assessed
as a change in the rate of disease progression, but this is limited to linear models. Most often
disease progression models are used separately on the treated arm and the placebo arm, and the
disease-modifying effect is assessed as the difference between the two progression models<fnref n="14" />,<fnref n="15" />.

Some classes of models are better suited to model disease-modifying mechanisms, namely dy-
namical models<fnref n="16" /> or generative models where we could simulate treatment effect as an intervention

<footnote n="9">⁹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.</footnote>

<footnote n="10">¹⁰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.</footnote>

<footnote n="11">¹¹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.</footnote>

<footnote n="12">¹²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.</footnote>

<footnote n="13">¹³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.</footnote>

<footnote n="14">¹⁴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.</footnote>

<footnote n="15">¹⁵C. S. Venuto, N. B. Potter, E. Ray Dorsey, et al., “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.</footnote>

<footnote n="16">¹⁶B. O. Taddé, H. Jacqmin-Gadda, J.-F. Dartigues, et al., “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.</footnote>
on the system<fnref n="17" />. On the other hand, non-linear mixed-effect model lack a general framework for such effect, and we will discuss about a possible choice of model in the disease-modifying model section.

## 4.2 Additive treatment effect model

Would it be possible to measure the effect of a treatment using a longitudinal observational cohort that lacks a control arm? To answer this, we propose a novel application of disease progression to assess the treatment effect in neurodegenerative diseases. In this section we provide a method to infer treatment effects by leveraging patient data in the absence of a placebo/control arm.

This model was motivated by the PPMI cohort, and more especially the way L-dopa treatment impacts Parkinson's disease. PD is characterized by the neurodegeneration of dopaminergic neurons, starting from the substantia nigra, and the subsequent loss of dopamine production in the brain. Dopamine is essential to the motor neurons, explaining the symptoms related to movement and tremor. The widely adopted treatment is the levodopa or L-dopa, the precursor of the dopamine. It provides the material needed to the brain of the patient in order to produce the lack-ing dopamine. This essentially results in a disappearing of motor issues, at least during a certain phase of the disease.

This treatment is purely symptomatic, as it does not help prevent further neurodegeneration. However it is highly efficient during the first years after diagnosis. Still, this yields a more simple task of modelling, since we do not have to take disease-modifying effects into account. Therefore the treatment effect can be directly measured at each timepoint, without requiring any longitudinal modelling.

Furthermore, PPMI provides us the very data we need for it: "OFF" and "ON" measures. Baseline visits are required to be clean of treatment, and enrolled patients should not have followed a long treatment prior to the baseline. Then patients can undergo dopaminergic therapy, the choice being to the discretion of the clinician and the patient. Values of dopaminergic treatment dose (L-dopa, dopamine agonists and others) are reported as L-dopa equivalent. For each subsequent visit, patient undergoing therapy will have two measures: one measure "OFF" deprived of treatment effect (in practice, patients are required to be clean of treatment for more than 6 hours) and a measure "ON" when the drug is effective (approximately one hour after taking medication). Even though the OFF state is not exactly the same as untreated PD, studies use it as a proxy to the state of a patient without any treatment<fnref n="18" />. Treatment effect can then be approximated as the raw difference between the ON state and the OFF state.

For this, we propose to leverage disease progression models. Our proposed approach combines them with a treatment model. The disease progression model computes the latent natural disease trajectory under the no-treatment hypothesis. The treatment model is a regression on the difference between the observed trajectory and the predicted one. One can interpret the natural disease progression predictions as a digital control twin of the patient who took the treatment.

Let's introduce notations and the mathematical formulation for this. Let $y_{ij} \in \mathbb{R}^d$ be the d-dimensional observation of patient $i$ at age $t_{ij}$, accompanied by a set of covariates $\mathbf{X}_{ij}$. Each patient $i$ has $n_i$ observations. We use the index $k$ for the dimension coordinates $y_{ijk}$. We introduce the

<footnote n="17">¹⁷Y. Iturria-Medina, F. M. Carbonell, R. C. Sotero, et al., “Multifactorial causal model of brain (dis)organization and therapeutic intervention: Application to Alzheimer’s disease,” en, *NeuroImage*, vol. 152, pp. 60–77, May 2017.</footnote>

<footnote n="18">¹⁸R. Cilia, E. Cereda, A. Akpalu, et al., “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.</footnote>
random variable $Y_{ijk}$ for the observations, whose distribution generated the observations. We will use the terms "ON state" for the trajectory of the patient under treatment, while "OFF state" will be acknowledged as the natural disease history state (which is not exactly true). The ON and OFF states are indicated by the treatment variable $T$: $T_{ij} = 1$ for the ON state or $T_{ij} = 0$ for the OFF state, and the realizations of the variable $Y_{ijk}$ will be superscripted with "ON" and "OFF" accordingly. Conditional average treatment effect in the most general case can therefore be written as the following:

$$CATE_k(\mathbf{x}) = \mathbb{E}(Y|T=1, \mathbf{X}=\mathbf{x}) - \mathbb{E}(Y|T=0, \mathbf{X}=\mathbf{x}) \quad (4.2.1)$$

If we had infinite samples for ON and OFF measures, we could simply estimate the CATE via the difference between the means of samples matching the covariate $x$. But this is not always possible, for instance if $x$ is a continuous variable. Instead, we will use a regression of $CATE_k(x)$ as a function of $x$, which corresponds to the covariates $\mathbf{X}_{ij}$. These covariates can include all types of internal variables, such as the OFF state for instance, and external cofactors such as sex or genetic risk factors. The covariates can be time-dependent as well.

In our observational cohort, we have visits for $T_{ij} = 1$ but rarely with $T_{ij} = 0$ as measure as the disease worsens. In the case where we have a model for disease progression, the OFF state is not a data, but the output of the model:

$$y_{ij}^{OFF} = \eta(\mathbf{z}_{ij}; \psi) + \epsilon_{ij} \quad (4.2.2)$$

where $\eta$ is the disease progression model, parametrized by $\psi$ and depending on individual time-dependent covariates $\mathbf{z}_{ij}$. We call $\eta$ the latent disease progression model. The noise structure $\epsilon_{ij}$ can be chosen depending on the structure of the data. In the continuous case we choose:

$$\epsilon_{ij} \sim \mathcal{N}(\mathbf{0}_d, \Sigma_d)$$

with $\Sigma$ diagonal for the sake of simplicity, but more complex structures can be introduced: time-dependent noise with Gaussian processes, a different noise covariance matrix... We then can infer the approximated treatment effect values:

$$\Delta_{ijk} = y_{ijk}^{ON} - \eta_k(t_{ij}, x; \mathbf{z}_i, \theta) \quad (4.2.3)$$

We will now regress the CATE using these approximated values. What we call the *treatment effect model* writes:

$$\phi(\mathbf{X}_{ij}; \omega) = \Delta_{ij} + \epsilon'_{ij} \quad (4.2.4)$$

where $\omega$ are the parameters of the regression and $\epsilon'$ is the residual noise. $\phi$ is the difference between the value of the marker observed under treatment and the *latent disease progression model* which is the expected value under no treatment. Our method is pictured in Figure 4.2.

Plugging (4.2.2) and (4.2.4) together we obtain:

$$\mathbf{y}_{ij}^{ON} = \eta(\mathbf{z}_{ij}; \psi) + \phi(\mathbf{x}; \omega) + \epsilon_{ij} - \epsilon'_{ij} \quad (4.2.5)$$

Since we consider $\psi$ and $\epsilon_{ij}$ to be learned already, we can assume that $\epsilon'_{ij}$ will summarize the residual noise of both the disease progression model $\eta$ and the treatment effect model $f$. In the
<figure><img src="image_39.png" /><figcaption>Figure 4.2: Model scheme. The plain curve is the progression of one individual, in blue before treatment starts and in orange after. The dashed blue curve is the predicted progression of the disease without treatment. Treatment effect $\phi$ can be estimated as the difference between the dashed blue line and the orange line, as a function of variables $x$.</figcaption></figure>

case where $\epsilon_{ij}$ has a specific bias (for instance the model has a hard time fitting young individuals), we can decompose this noise as the sum of a structured bias and pure noise. In this case the bias will be included in the learning of $\phi$. The condition for $\phi$ to be an unbiased estimator of CATE is then:

$$ \forall x, \mathbb{E}(Y_{ijk} | T_{ij} = 0, \mathbf{X}_{ij} = \mathbf{x}) - \eta_k(\mathbf{z}_{ij}; \psi) = 0 \quad (4.2.6) $$

In practice, this means that the predictions of our progression model should not be biased when conditioned on covariates. Therefore we need to verify after the calibration of the $\eta$ model that the residuals are independent of the covariates used in the treatment effect model $\phi$. It is often easy to validate this condition for covariates that are external to the progression model, such as sex or ethnicity. However the model predictions are usually very dependent on other covariates, such as age or time since diagnosis. Indeed it is well known that prediction models have a hard time predicting a long time ahead. DCM for instance is known to make long-term predictions that are prone to under-estimation.

However we can also understand condition 4.2.6 locally: for all $x$ for which the disease progression models are unbiased, then $\phi(\mathbf{x})$ can be considered as an unbiased estimate of CATE($\mathbf{x}$). This is less constraining as we can empirically select the range of values for $x$ for which we make accurate predictions.

We now present the models used in our application for $\eta$ and $\phi$.
## 4.2.1 Latent disease progression model

Any longitudinal model is suitable in this application, this is not specific to the DCM. Indeed the latent disease progression model can be any disease progression model as long as it provides an unbiased estimation of $E(Y_{ijk}|T_{ij} = 0, x_{ij} = x)$, so that the inferred values $\tilde{y}_{ijk}$ can be unbiased estimators of CATE. Therefore we will make several applications: two with the DCM model, and one with another disease progression model.

For the second disease progression model, we will use a dynamical model fitted using the penalty method. This application was done in collaboration with Pr. Bruno Jedynak and PhD Michael Wells from Portland State University who developed this method. The penalty method is an iterative method for fitting a system of Ordinary Differential Equations (ODE) to data. Specifically, we fit an ODE of the form:

$$ \frac{\partial}{\partial t} \eta_k(t) = f(\eta_k(t)) \qquad (4.2.7) $$

It minimizes a cost over a set of discretized trajectories and functions $f$ in a space of smooth functions, specifically a reproducing kernel Hilbert space. The cost function has three terms: a data term, a penalty term, and a regularization term. The data term keeps the trajectories close to the data. The penalty term ensures the trajectories follow the ODE (4.2.7). The regularization term corresponds to the norm of the function $f$, called the field, in the RKHS and it is meant to prevent overfitting. The penalty term is a loose constraint on the ODE, namely $||\frac{\partial \eta_k}{\partial t} - f(\eta_k)|| \le \frac{1}{\lambda}$, where $\lambda$ is increased at each iteration. The optimal fit is found after several iterations when the algorithm converges. Details and validation of the algorithm are provided in<fnref n="19" />.

## 4.2.2 Treatment effect model

We chose the treatment effect model $\phi: \mathbf{x} \mapsto \phi(\mathbf{x})$ to be Bayesian. Indeed we are interested in a descriptive model for the treatment effect, and a Bayesian framework provides reliable estimates with credibility intervals. We used three treatment effect models:

* Bayesian linear regression (BLR): a model able to capture linear relationships between predictors and outcome, but which does not capture interactions between two input variables $x_l$ and $x_{l'}$, thus limiting the descriptive power of the model

* Bayesian additive regression trees (BART): the state-of-the-art model for CATE regressions, which is flexible enough to capture any relationship between the covariates $\mathbf{x}$ and the outcome. However, in practice the trees did not embed significant interactions between input variables so we added a third model

* Bayesian kernel machine regression (BKMR): a multivariate kernel regression model allows for a lot of flexibility due to the choice of the kernel and the hyperparameters associated with it. This model, although harder to fit and more dependent on variable selection, has the advantage of modeling interactions between predictors which is very valuable from a descriptive point of view

<footnote n="19">¹⁹K. Lahouel, M. Wells, V. Rielly, et al., Learning nonparametric ordinary differential equations from noisy data, Feb. 2023.</footnote>
BLR and BART allow the computation of variable importance as a percentage of the explained outcome variance. Since BKMR is less performant when the dimension of x increases, we used BLR and BART for variable selection to reduce dimensions to a meaningful set of covariates x|S. More details about each method are provided next.

## Bayesian linear regression

A deeper introduction to BLR models can be found in<fnref n="20" />. With our notations the model takes the following form (in a single dimension for the sake of simplicity, extension to multidimensional output is straightforward):

$$ \phi(\mathbf{x}; \alpha, \beta) = \alpha^T \mathbf{x} + \beta \qquad (4.2.8) $$

with $\alpha \in \mathbb{R}^m$ with $m$ being the dimension of the covariates $\mathbf{x}_{ij}$ and $\beta \in \mathbb{R}$. The BLR model describes $\alpha$ and $\beta$ as random variables for which we provide normal priors $\alpha_l \sim \mathcal{N}(0, \sigma_l)$, $\forall l \in [|1, m|]$ and $\beta \sim \mathcal{N}(0, \sigma_\beta)$ where $(\sigma_l)_{l}$ and $\sigma_\beta$ are hyperparameters. We decided to fix those hyperparameters to values depending on the training data: $\sigma_\beta = \text{Var}(y_{ij})$, $\sigma_l = \frac{\text{Cov}(x_{ijl}, y_{ijl})}{\text{Var}(x_{ijl})}$ (which would be the slope coefficient in a non-Bayesian setting).

## Bayesian kernel machine regression

First we will describe standard kernel regression. Here we use a reproducing kernel Hilbert space $\mathcal{H}$, $f \in \mathcal{H}$ with a kernel $k$. In the applications we chose the standard gaussian kernel $k(\mathbf{x}, \mathbf{x}') = \exp(-\frac{\|\mathbf{x}-\mathbf{x}'\|^2}{\sigma^2})$. With the representer theorem, $\phi$ can be expressed as a function of the observations:

$$ \phi(\mathbf{x}) = \sum_{ij} w_{ij} k(\mathbf{x}, \mathbf{x}_{ij}) \qquad (4.2.9) $$

In order to keep the dimension of the problem within a reasonable range and to avoid ill-conditioned kernel matrices, we select only a few control points $(\mathbf{x}_i^c)_{i=1...n_c}$. These points can be placed on a grid in $\mathbb{R}^m$ if $m$ is the dimension of $\mathbf{x}$ for a regular spacing, or we can choose to subsample the data $(\mathbf{x}_{ij})_{i=1...n, j=1...n_i}$ using the Kmeans++ initialization method<fnref n="21" />. Hence:

$$ \phi(\mathbf{x}) = \sum_{i} w_i k(\mathbf{x}, \mathbf{x}_i^c) \qquad (4.2.10) $$

The least squares problem can be formulated as an optimization problem (every uppercase corresponds to the vectorization of the lowercase variable):

$$ \min_{i=1...n_c, w_i \in \mathbb{R}^d} \sum_{ij} \|\mathbf{y}_{ij}^{ON} - \eta(\mathbf{z}_{ij}; \psi) - \phi(\mathbf{x}_{ij})\|^2 + \lambda \|f\|_{\mathcal{H}}^2 \qquad (4.2.11) $$

$$ \min_{W \in \mathbb{R}^{n_c \times d}} \|Y^{ON} - \eta - KW\|^2 + \lambda \sum_{i} W_{:,i}^T K^{c} W_{:,i} \qquad (4.2.12) $$

<footnote n="20">²⁰A. Gelman, J. B. Carlin, H. S. Stern, et al., “Bayesian Data Analysis,” en, Chapman and Hall/CRC, Nov. 2013.</footnote>

<footnote n="21">²¹D. Arthur and S. Vassilvitskii, “K-means++: The Advantages of Careful Seeding,” en,</footnote>
where $\mathbf{y}_{ij}^{ON}$ refer to the observations in the ON state, $K = (k(\mathbf{x}_{ij}, \mathbf{x}_i^c)) \in \mathbb{R}^{\sum_i n_i \times n_c}$ and $K^c = (k(\mathbf{x}_i^c, \mathbf{x}_i^c)) \in \mathbb{R}^{n_c \times n_c}$.

The solution is explicit as this problem is convex and the gradient can be computed to look for the critical point, which yields:

$$W^* = (K^T K + \lambda K^c)^{-1} (Y^{ON} - \eta) \qquad (4.2.13)$$

BKMR<sup>22</sup> is the result of adding a Bayesian framework to kernel regression. This is very similar to Gaussian processes, but here we present the model from the kernel regression point of view so we refer the reader to<sup>23</sup> for equivalences. We describe a single-dimension setting, as the extension to multidimensional output is straightforward.

The function $\phi$ is still approximated in a reproducing kernel Hilbert space (RKHS) where the basis functions are $\{k(\mathbf{x}, \cdot), \mathbf{x} \in \mathbb{R}\}$ for a chosen kernel $k: \mathbb{R}^m \times \mathbb{R}^m \to \mathbb{R}$ which is symmetric positive-definite. The representer theorem implies that the minimizing function in the RKHS is of the following form $\sum_{ij} \alpha_{ij} k(\mathbf{x}_{ij}, \cdot)$. We add a random intercept to the model which yields:

$$\phi(\mathbf{x}; \alpha, \beta) = \sum_{ij} \alpha_{ij} k(\mathbf{x}_{ij}, \mathbf{x}) + \beta \qquad (4.2.14)$$

with $\alpha \in \mathbb{R}^N$ with $N$ being the total number of observations and $\beta \in \mathbb{R}$. In the Bayesian setting, we add priors to those parameters: $\alpha_{ij} \sim \mathcal{N}(0, \sigma_\alpha)$, $\forall i \in [|1, n|]$, $\forall j \in [|1, n_i|]$, $\sigma_\alpha \sim \mathcal{N}(0, \bar{\sigma}_\alpha)$ and $\beta \sim \mathcal{N}(0, \sigma_\beta)$.

As for BLR we define the hyperparameter $\sigma_\beta = \text{Var}(y_{ij})$ and fix $\bar{\sigma}_\alpha = 10^3$. The parameter $\sigma_\alpha$ which corresponds to the regularization strength is automatically adjusted thanks to the Bayesian setting.

We chose the gaussian kernel, which formulates as:

$$k(\mathbf{x}, \mathbf{x}') = \prod_{l=1}^{m} \exp \left( \frac{(x_l - x_l')^2}{\sigma_l^2} \right) \qquad (4.2.15)$$

with $\sigma_l$ being the bandwidth for the l-th covariate. These bandwidths are hyperparameters. All hyperparameters were optimized using nested cross-validation.

However, due to the non-continuous aspect of some of the covariates used in our application, one can argue that the Gaussian kernel is not always adapted. Some ideas to go further would be to adapt kernels to the type of data used in the covariates, for instance a specific kernel for discrete data.

### Bayesian additive regression trees

BART was introduced in<sup>24</sup>. The model consists in a sum-of-trees:

<footnote n="22"><sup>22</sup>J. F. Bobb, L. Valeri, B. Claus Henn, et al., “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.</footnote>

<footnote n="23"><sup>23</sup>M. Kanagawa, P. Hennig, D. Sejdinovic, et al., Gaussian Processes and Kernel Methods: A Review on Connections and Equivalences, Jul. 2018.</footnote>

<footnote n="24"><sup>24</sup>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.</footnote>
$$ \phi(\mathbf{x}; \mathcal{T}, \mu) = \sum_{l=1}^{m} g(\mathbf{x}; \mathcal{T}_l, \mu_l) \qquad (4.2.16) $$

where m is the number of trees, $\mathcal{T}_l$ is the structure of the l-th tree (with interior nodes decision rules) and $\mu_l$ its set of leaf values, and g is the function returning the output of a given tree with specified parameters applied to an input $\mathbf{x}$. Specific priors are put on $\mathcal{T}_l$ and $\mu_l$ while the number of trees m and the parameter controlling depth $\alpha$ are hyperparameters. We refer the reader to the source paper for more details.

### 4.2.3 Estimation method

End-to-end estimation of the two models $\eta$ and $\phi$ would be too complicated since the values for the regression are inferred from the latent disease progression model. Therefore we used a two-step estimation.

First, the latent disease progression model needs to be calibrated on data coming from untreated individuals. In our case, this data corresponds to the observations of the patients prior to the beginning of the treatment. In most other applications where this data is unavailable, one can estimate the model parameters on a different dataset. For instance, in a new clinical trial, one could use a model estimated thanks to the control arm of a previous trial or thanks to an observational cohort of the disease progression. For the DCM, calibration is performed with the Monte-Carlo Markov chain stochastic approximation expectation maximization algorithm (MCMC-SAEM)$^{25}$. It results in estimating the population parameters $\mathbf{z}_{pop}$ on a training set. From there we can compute the individual parameters $\mathbf{z}_i$ of a new patient based on its inclusion visit or any amount of points before the beginning of the treatment while keeping $\mathbf{z}_{pop}$ fixed. Individual parameters are considered covariates that can be fed into the regression.

Since DCM is a Bayesian method, we have posterior distributions for parameters and predicted observations $\eta_k(t_{ij}, x; \mathbf{z}_i, \mathbf{z}_{pop})$. However, computing the regression based on distributions instead of point estimates is too expensive computationally as it would require many Monte-Carlo samples of those variables. We opted to use the maximum a posteriori (MAP) estimates for both $\mathbf{z}_i$ and $\eta_k(t_{ij}, x; \mathbf{z}_i, \mathbf{z}_{pop})$.

For the ODE model, the penalty method described previously is applied for the fit.

The Bayesian regression models were all estimated through Monte-Carlo sampling, with a No U-Turn sampler (NUTS)$^{26}$ for BLR and BKMR, and an ad-hoc sampler for BART$^{27}$.

The models were implemented as described next using the probabilistic library PyMC3 in Python for Bayesian models$^{28}$. We used an open-source implementation of BART in PyMC3$^{29}$. The penalty method was implemented in Python with the TensorFlow library.

<footnote>$^{25}$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>$^{26}$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>$^{27}$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.</footnote>

<footnote>$^{28}$J. Salvatier, T. V. Wiecki, and C. Fonnesbeck, “Probabilistic programming in Python using PyMC3,” en, *PeerJ Computer Science*, vol. 2, e55, Apr. 2016.</footnote>

<footnote>$^{29}$M. Quiroga, P. G. Garay, J. M. Alonso, et al., *Bayesian additive regression trees for probabilistic programming*, Oct. 2022.</footnote>
### 4.2.4 Experiments

#### Simulations

We already know that the disease progression models that we use are validated in previous works (³⁰ for the DCM and<fnref n="31" /> for the ODE model). Therefore we just need to validate the second part of the method proposed above. We performed some simulations to better understand how well the treatment effect was estimated under various scenarios.

**Scenario 1: True predictors** In this first scenario, we generated 10000 random observations $X$ of dimension 2 and assumed that the first dimension was the natural disease measure (OFF state). Then we computed the treatment effect as a deterministic function $f$ of $X$, mimicking the situation where the treatment effect would be increasing as the disease increases and as another covariate increases (which could be thought of as the treatment dose). We decided that $f$ would be linear so that the Bayesian linear regression would be able to fit the outcome: $f(x_1, x_2) = 10x_1 + x_2$. Finally we add a Gaussian noise to the output $Y = f(X) + \epsilon_Y$. In this very simple case we first recorded the $R$ squared metric for the different regressions on a test set while progressively increasing the noise level $\epsilon_Y$ in the training set. Results are shown in Table 4.1 below. This simply shows that regression models perform very well until a certain amount of noise, which is around the standard deviation of the treatment effect $Y$.

<table><caption>Table 4.1: $R$-squared for the regression models on simulated data depending on the noise level of the outcome. The noise level is multiplicative factor between the standard deviation of the noise and the standard deviation of the input data $X$.</caption><thead><tr><th>Noise level</th><th>0.0</th><th>0.1</th><th>1.0</th><th>10.0</th><th>100</th></tr></thead><tbody><tr><th>BLR</th><td>1.0</td><td>1.0</td><td>1.0</td><td>0.98</td><td>0.66</td></tr><tr><th>BKMR</th><td>0.99</td><td>0.98</td><td>0.98</td><td>0.94</td><td>0.41</td></tr><tr><th>BART</th><td>1.0</td><td>1.0</td><td>0.99</td><td>0.97</td><td>0.61</td></tr></tbody></table>

**Scenario 2: Noisy disease progression model** In this second scenario, we use the same data as previously, with a noise level on the outcome of 0.1 allowing the regressions to find the true CATE if the predictors are true. Except that the predictors coming from the disease progression model are noisy, that is we add a Gaussian centered noise on $x_1$: $\tilde{x}_1 = x_1 + \epsilon_X$. We can think of this as the uncertainty of the predictions of the disease progression model. Notice that in this case the predictors are still unbiased. We then have the true treatment effect $Y = f(X) + \epsilon_Y = y^{ON} - x_1$ and the inferred treatment effect which is based on the disease model predictions $\hat{Y} = y^{ON} - \tilde{x}_1$ (cf equation (4.2.3)). The regressions are trained to predict $\hat{Y}$ from $(\tilde{x}_1, x_2)$, but we test them on their ability to reconstruct the real treatment effect $Y$. Table 4.2 records the $R$-squared as we increase the noise level $\epsilon_X$ on the predictors. The results show that under a certain threshold, this noise is not an issue and the regression still perform very well.

<footnote n="30">³⁰I. Koval, A. Bône, M. Louis, et al., “AD Course Map charts Alzheimer’s disease progression,” eng, *Scientific Reports*, vol. 11, no. 1, p. 8020, Apr. 2021.</footnote>

<footnote n="31">³¹K. Lahouel, M. Wells, V. Rielly, et al., *Learning nonparametric ordinary differential equations from noisy data*, Feb. 2023.</footnote>
<table><caption>Table 4.2: R-squared for the regression models on simulated data depending on the noise level of the predictors.</caption><tr><td>Noise level</td><td>0.0</td><td>0.1</td><td>1.0</td><td>10.0</td><td>100</td></tr><tr><td>BLR</td><td>1.0</td><td>1.0</td><td>0.99</td><td>0.80</td><td>0.03</td></tr><tr><td>BKMR</td><td>0.98</td><td>0.98</td><td>0.98</td><td>0.80</td><td>0.05</td></tr><tr><td>BART</td><td>1.0</td><td>1.0</td><td>0.99</td><td>0.83</td><td>0.10</td></tr></table>

**Scenario 3: Biased disease progression model** In this last scenario, we add a noise to the disease progression model $x_1$ so that its predictions are biased. We tried several types of bias. The first case is a constant bias: $\tilde{x}_1 = x_1 + C + \epsilon_X$ with a certain constant $C$. Then we add a bias conditionally to the other covariate $x_2$ used in the regression: $\tilde{x}_1 = x_1 + 0.1x_2 + \epsilon_X$. The average marginal effect of each predictor is then estimated with its 95% credibility interval and we compare it to the ground truth in Figure 4.3.

The first and third rows of Figure 4.3 highlight the bias effect directly on the marginal effect of $x_1$. In both cases the regressions have learned the bias. This is clearly visible in the second row, where we see that the marginal effect on $x_2$ has learned a constant bias (which is exactly equal to the constant $C$). We also observe the high fluctuations of the estimated effect for $x_2$ on the BKMR and BART models, which is due to the smaller effect of $x_2$ relative to $x_1$ in $Y$, therefore the noise is felt more for this predictor.

What these simulations show is an overall high robustness to the noise, while being susceptible to the bias in the disease progression model predictions. Even in this case, the regressions learn with the bias, which means that we can still interpret the results if we have an idea of where the bias come from.
<figure><img src="image_40.png" /><figcaption>Figure 4.3: Marginal treatment effect compared to ground truth. First two rows correspond to the constant bias $\tilde{x}_1 = x_1+C+\epsilon_X$ and last two rows correspond to the linear bias $\tilde{x}_1 = x_1+0.1x_2+\epsilon_X$. Each column correspond to a regression model. The regression mean posterior value is plotted in plain blue line with the 95% credibility interval while the ground truth is plotted in red.</figcaption></figure>
## PPMI dataset

Our application concerns Parkinson's disease (PD), the second most prevalent neurodegenerative disease. It is characterized by the loss of neurons responsible for dopamine production in the brain, especially in the substantia nigra. Lack of dopamine leads mainly to motor symptoms, but this can be treated by a dopamine replacement strategy. Dopaminergic therapy has shown to be an efficient symptomatic treatment, even if no disease-modifying effect has been found<fnref n="32" />. All patients in the cohort follow this treatment.

Data used in the preparation of this article were obtained from the Parkinson's Progression Markers Initiative (PPMI) database<fnref n="33" /> (www.ppmi-info.org/access-data-specimens/download-data). For up-to-date information on the study, visit ppmi-info.org. <fnref n="34" />

The dataset comprises PD patients included within two years of their diagnosis. We selected all patients with at least two visits in the study for the longitudinal history, which resulted in 900 subjects for a total of 7918 visits. We included basic covariates such as sex, age, weight, BMI, treatment category, and dose and a clinical assessment: the MDS-UPDRS score. It is subdivided in 4 sections: non-motor aspects of experiences of daily living (I); motor aspects of experiences of daily living (II); motor examination (III) by the clinician; motor complications (IV) which measure the strength of secondary effects due to the dopaminergic treatment. As mentioned earlier, the third part is measured twice for a patient under treatment: an "OFF" version is assessed when the treatment effect has worn off; an "ON" version is measured shortly after the patient has taken its medication. This particularity of the cohort allows us to have direct access to the values of $Y_{ijk}|T_{ij} = 0$.

Dataset demographics can be found in Table 4.3.

Clinical imaging biomarkers are extracted from DAT-scan with single-photon emission computed tomography (SPECT) which traces the dopamine transporter molecule. The values are recorded for the areas where the signal is the most important in PD, which is the striatum. The caudate and putamen are subparts of the striatum. We refer to them not as left and right but ipsilateral or contralateral depending on the side of onset of the motor symptoms. 1704 visits were recorded with DAT-scan and but only 135 patients had imaging data prior to the start of treatment.

## Experiment on motor score

Our first experiment uses the MDS-UPDRS part III (motor score) as the outcome, which ranges from 0 to 132. We calibrated a DCM model on observations pre-treatment as our latent disease progression model. We then predicted the evolution of this motor score under no treatment with the DCM model. Since the cohort recorded the OFF-state values of the MDS-UPDRS part III, we can measure the prediction performance. We obtain a mean absolute error of 6.9, which has to be compared to the standard measurement error of this clinical score (the inherent noise), which according to<fnref n="35" /> is at 4.31. We estimated a BLR and a BART model on the set of observed covariates

<footnote n="32">³²C. Lungu, J. M. Cedarbaum, T. M. Dawson, et al., "Seeking progress in disease modification in Parkinson disease," en, *Parkinsonism & Related Disorders*, vol. 90, pp. 134-141, Sep. 2021.</footnote>

<footnote n="33">³³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.</footnote>

<footnote n="34">³⁴PPMI – a public-private partnership – is funded by the Michael J. Fox Foundation for Parkinson's Research and funding partners, including [list the full names of all of the PPMI funding partners found at www.ppmi-info.org/about-ppmi/who-we-are/study-sponsors].</footnote>

<footnote n="35">³⁵P. Martinez-Martin, C. Rodriguez-Blazquez, M. Alvarez-Sanchez, et al., "Expanded and independent validation of the Movement Disorder Society-Unified Parkinson's Disease Rating Scale (MDS-UPDRS)," en, *Journal of*</footnote>
<table><caption>Table 4.3: Demographics table for PPMI variables. Count is the number of visits where the variable is available. MDS-UPDRS: Movement Disorder Society - Unified Parkinson's Disease Rating Scale. TD: Tremor Dominant. PIGD: Postural Instability and Gait Disorder. Medication category: 0 = No medication, 1 = Levodopa, 2 = Dopamin agonist, 3 = Other, 4 = Levodopa and other, 5 = Levodopa and dopamin agonist, 6 = Dopamin agonist and other, 7 = Levodopa, dopamin agonist and other.</caption><thead><tr><th></th><th>MDS-UPDRS part I</th><th>MDS-UPDRS part II</th><th>MDS-UPDRS part III</th><th>MDS-UPDRS part IV</th><th>TD score</th></tr></thead><tbody><tr><td>name</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>count</td><td>7897</td><td>7897</td><td>7519</td><td>7922</td><td>7922</td></tr><tr><td>mean</td><td>7.998607</td><td>8.729264</td><td>25.512701</td><td>1.222545</td><td>5.447740</td></tr><tr><td>std</td><td>5.562359</td><td>6.387908</td><td>12.429672</td><td>2.494955</td><td>4.544258</td></tr><tr><td>min</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>median</td><td>7</td><td>7</td><td>24</td><td>0</td><td>5</td></tr><tr><td>max</td><td>38</td><td>48</td><td>100</td><td>17</td><td>30</td></tr><tr><td>name</td><td>PIGD score</td><td>Treatment dose</td><td>Cumulative treatment dose</td><td>Medication category</td><td>Sex</td></tr><tr><td>count</td><td>7922</td><td>5757</td><td>7922</td><td>6862</td><td>7922</td></tr><tr><td>mean</td><td>2.086342</td><td>546.896595</td><td>885.235715</td><td>-</td><td>55%M/45%F</td></tr><tr><td>std</td><td>2.239677</td><td>808.732513</td><td>1660.498211</td><td>-</td><td>-</td></tr><tr><td>min</td><td>0</td><td>0</td><td>0</td><td>0</td><td>-</td></tr><tr><td>median</td><td>1</td><td>421.2</td><td>325.0</td><td>2</td><td>-</td></tr><tr><td>max</td><td>16</td><td>24165</td><td>42400</td><td>7</td><td>-</td></tr><tr><td>name</td><td>Age since baseline</td><td>Weight</td><td>Body Mass Index</td><td></td><td></td></tr><tr><td>count</td><td>7922</td><td>3732</td><td>3728</td><td></td><td></td></tr><tr><td>mean</td><td>2.30794</td><td>78.506645</td><td>26.698898</td><td></td><td></td></tr><tr><td>std</td><td>2.07245</td><td>16.847416</td><td>4.677159</td><td></td><td></td></tr><tr><td>min</td><td>0</td><td>37.3</td><td>13.818027</td><td></td><td></td></tr><tr><td>median</td><td>1.625</td><td>77.8</td><td>26.178115</td><td></td><td></td></tr><tr><td>max</td><td>9.125</td><td>162</td><td>61.728395</td><td></td><td></td></tr></tbody></table>
to perform variable selection since BKMR does not fare well with too many covariates.

<figure><img src="image_41.png" /><figcaption>Figure 4.4: Variable importance in BART regression of ground truth treatment effect (ON measure minus OFF measure). MDS[X]: part X of the MDS-UPDRS. The baseline score corresponds to the last value of the score before treatment starts. LEDD: levodopa equivalent daily dose (levodopa is the main dopaminergic medication). TD SCORE: tremor dominant score. PIGD SCORE: postural imbalance and gait disorder score. TIME SINCE DOSE [ON/OFF]: time since the last medication at the moment of the examination [ON/OFF]. BMI: body mass index. PD MED: type of dopaminergic treatment. MDS1_11: gastrointestinal complications. AGE SINCE BL: time since baseline.</figcaption></figure>

Our variable selection relied on Figure 4.4. Using the whole range of covariates the model highlighted the MDS-UPDRS parts III and IV as the most important variables. Unfortunately, part IV only records the secondary effects of treatment, and a disease progression model trained on pre-treatment data will only predict zeros for the value of this covariate under the no-treatment hypothesis. However, this score is quite correlated to the treatment dose (LEDD) and the fact that LEDD is the third most important variable suggests that this variable be used instead for the regression. The $R^2$ with the full set of covariates reached 0.86 while reducing to those two only dropped to 0.82. The dose is interesting all the more so as we can *act* on this covariate, which means that we can *control* treatment effect. We thus reduced the covariates to the treatment dose and the MDS-UPDRS part III OFF-state representing the current disease state. This MDS-UPDRS part III is either observed directly, using the OFF measures or estimated, with the DCM.

We verified that the DCM model predictions were not biased conditionally to the selected covariates. There was no bias for the treatment dose, but a small bias for the predictions conditioned on the predicted value, as the model tends to under estimate at large values. However it does not affect the range of values where most of the patients are.

We estimated the three treatment effect models with two different settings. First, we regressed the observed treatment effect $\phi_{GT}$, which is the measured ON value minus the measured OFF value. We then reconstructed the ON score $\hat{y}_{ON} = y_{OFF} + \phi_{GT}(\mathbf{x}_{GT})$. Using the ground truth treatment effect and predictor provides an upper bound on the best metric scores we can achieve using inferred
<figure><img src="image_42.png" /><figcaption>Figure 4.5: Model illustration on data from a PPMI subject. DCM model is personalized with pre-treatment data (crosses) to predict post-treatment latent disease progression (dashline). Setting 1 ($\phi_{GT}$) models the difference between OFF and ON measures. Setting 2 ($\phi_{pred}$) models the difference between predicted OFF and ON measures.</figcaption></figure>

values for treatment effect and predictors, which was our second setting. In this second part, the regressions $\phi_{pred}$ used the predicted OFF value and the treatment dose as input and the inferred treatment effect (measured ON value minus predicted OFF value) as output. The reconstructed outcome becomes $\hat{y}_{ON} = \eta(t) + \phi_{pred}(\mathbf{x}_{pred})$. Figure 4.5 illustrates these two settings on actual data from PPMI.

We also provided a third setting for comparison without using the latent disease progression model predictions. We use the common approach used in treatment models, using baseline variables. This provides a lower bound error target on the reconstruction metrics. Here we estimated the regressions to predict the ON score directly. As predictors, we used the treatment dose, the baseline OFF score, i.e., the last measure without treatment, and the time since treatment started. We reported the prediction metrics on the reconstructed ON score in Table 4.4.

Compared to setting 1 using ground truth treatment effect and predictors, setting 2 is worse as it is expected since it uses inferred predictors and treatment effect. But setting 2 only lost around 2.3 points in mean absolute error while using the score inferred with DCM which had about 7 points of absolute error. It shows that the treatment effect models compensated for the DCM prediction errors. Comparing settings 2 and 3 we see that leveraging the information from the longitudinal predictive latent disease progression model improved the metrics by a significant margin, yet still far from the ground truth.

Overall BART performed best out of the three treatment effect models. The difference in maximum a posteriori predictions was low, but the log-likelihood was much better because BART has a much tighter credible interval for its predictions. However, the cost was paid in runtime.

Figure 4.6 shows the posterior distributions of the treatment effect models. The top and middle rows correspond to the treatment effect $\phi$ learned in settings 1 and 2 respectively. The treatment effect in inference setting 2 was underestimated by a shift of approximately 2 points, but it captured the same trends. The 2D interaction plot shows how the dose affects the treatment effect depending on the disease state. Specifically, the best treatment effect is obtained with an increasing dose as the disease progresses toward later stages.
<table><caption>Table 4.4: Reconstruction metrics for the ON-state MDS-UPDRS part III computed with four folds. Setting 1: ground truth (GT) treatment effect and ground truth covariates are used. This is the best possible performance. Setting 2: inferred treatment effect and predicted (pred) covariates are used. This is our method. Setting 3: baseline (BL) predictors are used with the time since baseline. This is the standard method. BLR: Bayesian Linear Regression. BKMR: Bayesian Kernel Machine Regression. BART: Bayesian Additive Regression Trees. $R^2$: R-squared coefficient of maximum a posteriori (MAP) prediction. MAE: Mean Absolute Error of MAP prediction. LL: posterior log-likelihood.</caption><thead><tr><th>Setting</th><th>Regression</th><th>R<sup>2</sup></th><th>MAE</th><th>LL</th><th>Runtime (s)</th></tr></thead><tbody><tr><td rowspan="3">1 (GT)</td><td>BLR</td><td>0.81±0.02</td><td>5.19±0.10</td><td>-6150±527</td><td>37.9±2.1</td></tr><tr><td>BKMR</td><td>0.81±0.03</td><td>5.11±0.09</td><td>-501±38</td><td>94.9±3.4</td></tr><tr><td>BART</td><td>0.82±0.03</td><td>5.11±0.10</td><td>-114±13</td><td>5240±294</td></tr><tr><td rowspan="3">2 (pred)</td><td>BLR</td><td>0.55±0.02</td><td>7.47±0.14</td><td>-7170±954</td><td>38.1±2.6</td></tr><tr><td>BKMR</td><td>0.56±0.02</td><td>7.42±0.13</td><td>-535±21</td><td>108±13</td></tr><tr><td>BART</td><td>0.57±0.02</td><td>7.37±0.19</td><td>-117±18</td><td>5190±110</td></tr><tr><td rowspan="3">3 (BL)</td><td>BLR</td><td>0.50±0.04</td><td>8.1±0.1</td><td>-5770±749</td><td>41.7±9.4</td></tr><tr><td>BKMR</td><td>0.48±0.03</td><td>8.2±0.2</td><td>-188±23</td><td>166±28</td></tr><tr><td>BART</td><td>0.52±0.03</td><td>7.98±0.16</td><td>-98±11</td><td>4550±143</td></tr></tbody></table>

Figure 4.6 only showed the mean effect of the dose on the treatment effect. Figure 4.7 includes the MDS-UPDRS part III (used as a proxy for the disease state). Notice how the effect seems to be proportional. This suggests an improvement of the model in this specific case of the MDS-UPDRS part III: measuring the treatment effect not as a difference between ON and OFF states, but as a difference relative to OFF score. The slope indicates that we observe a mean 35% improvement of the motor score under treatment.
<figure><img src="image_43.png" /><figcaption>Figure 4.6: Top row: posterior distribution of dose on real treatment effect (ON minus OFF) according to the three treatment effect models. The plain line is the mean, the colored area is the 95% credibility interval. Middle row: same plot but with inferred treatment effect (ON minus predicted OFF). Bottom row: interaction plot, displays the mean treatment effect learned by each regression as a function of the predicted disease state (predicted MDS-UPDRS part III OFF) and the dose.</figcaption></figure>
<figure><img src="image_44.png" /><figcaption>Figure 4.7: Top row: posterior distribution of disease state on real treatment effect (ON minus OFF) according to the three treatment effect models. The plain line is the mean, the colored area is the 95% credibility interval. Bottom row: same plot but with inferred treatment effect (predicted ON minus OFF).</figcaption></figure>
### Experiment on quality of life score

Our second experiment focused on the MDS-UPDRS part II, which ranges from 0 to 42. This score evaluates patients' quality of life (QoL), and QoL improvement is becoming one of the main criteria for treatment acceptance. However, there was no OFF measure for this score so the treatment effect could only be inferred. We thus replicated settings 2 and 3 from the last experiment. The best treatment effect model was BART in setting 2, with $R^2$:0.83 and MAE:2.56 on the reconstructed MDS-UPDRS part II in setting 2. This is much better than the best regression in setting 3 with $R^2$:0.71 and MAE:3.31. Table 4.5 reports full metrics on the reconstruction of the MDS-UPDRS part II.

The interest of the model, besides being able to predict the score of patients, is to understand how covariates influence the treatment effect on QoL. Our analysis of the full set of covariates revealed that the most important variables were the QoL without treatment as predicted by the DCM and the PIGD score (postural imbalance and gait disorder) as PIGD patients have a worsening of their QoL under treatment as opposed to TD (tremor dominant) patients. PIGD and TD are the two main sub-groups in Parkinson's disease, and TD patients tend to respond better to dopaminergic treatment than PIGD patients<fnref n="36" />.

<table><caption>Table 4.5: Reconstruction metrics for the MDS-UPDRS part II computed on a separate test set.</caption><thead><tr><th>Setting</th><th>Regression</th><th>R<sup>2</sup></th><th>MAE</th><th>LL</th><th>Runtime (s)</th></tr></thead><tbody><tr><td rowspan="3">2 (pred)</td><td>BLR</td><td>0.82</td><td>2.72</td><td>-1.64e+04</td><td>59.4</td></tr><tr><td>BKMR</td><td>0.80</td><td>2.65</td><td>-7.55e+02</td><td>989</td></tr><tr><td>BART</td><td><strong>0.83</strong></td><td><strong>2.56</strong></td><td><strong>-2.04e+02</strong></td><td>5e+03</td></tr><tr><td rowspan="3">3 (BL)</td><td>BLR</td><td>0.68</td><td>3.45</td><td>-1.69e+04</td><td>54.9</td></tr><tr><td>BKMR</td><td>0.63</td><td>3.56</td><td>-4.78e+02</td><td>913</td></tr><tr><td>BART</td><td>0.71</td><td>3.31</td><td>-2e+02</td><td>5e+03</td></tr></tbody></table>

### Experiment on imaging

We use the penalty method to fit an ODE to the imaging biomarkers. We use pre-treatment data only so that this model extrapolates the latent disease progression without the influence of treatment. The data and the fitted ODE curves are represented for the putamen in Figure 4.8. This figure shows that imaging data is quite noisy, so the model provides a smoothing of the trajectory.

We now have the predicted latent disease progression curve for each patient, which we will consider as our synthetic placebo arm. We can then compare it to the values observed under treatment. However, since imaging measures are noisy and sparse, we use a second ODE model to impute the trajectories of our treatment arm. This model uses the full set of points of patients, thus smoothing out the biomarkers and minimizing the noise. The difference between the treatment arm model and the placebo arm model is shown in Figure 4.9. The difference is consistent across individuals, which hints towards an effect of the dopaminergic treatment over the dopaminergic transporters in the putamen.

<footnote n="36">³⁶C. Marras, K. R. Chaudhuri, N. Titova, et al., “Therapy of Parkinson’s Disease Subtypes,” en, *Neurotherapeutics*, vol. 17, no. 4, pp. 1366–1377, Oct. 2020.</footnote>
<figure><img src="image_45.png" /><figcaption>Figure 4.8: Spaghetti plot for raw data (top) and model fit (bottom) on Putamen. Left column is ipsilateral (putamen on the side of motor signs onset), and right column is contralateral.</figcaption></figure>

In order to understand how the treatment effect is impacted by potential covariates we estimate the regression models. The CATE outcome is approximated as the difference between the treatment arm model and the placebo arm model as shown in figure 3. We started with all the regressors mentioned in the previous experiments plus the values of the imaging markers predicted by the placebo arm model.

Variable selection was first performed by computing variable importance with BART as shown in Figure 4.10. The selected variables which explained the most the treatment effect are the predicted values of the imaging biomarkers. We also included the treatment dose and the cumulative treatment dose as covariates of interest due to the possibility to control for those values.

The metrics for the regressions' fit are shown in Table 4.6. As we can see BART outperforms the other methods at the cost of a much higher runtime while the linear regression is not able to reconstruct the outcome as well, indicating that the interactions are mostly non-linear. To validate our approach, we also computed the reconstruction error between real values $y_{ijk}$ and model predictions $\eta_k(t_{ij}, x; z_i, \theta) + \phi(x_{ij}; \omega)$. For the putamen contralateral (resp. ipsilateral) the best method (BART) obtained a mean absolute error of 0.10 (resp. 0.17), which is good when compared to the noise of the real imaging markers.

To understand the interactions of these covariates on treatment effect we plot their average effect evaluated as the mean of the posterior of the models. Figure 4.11 shows the mean effect for some variables of interest. The first row tells us that the treatment effect is higher for low values of the putamen, which corresponds to later stages of the disease as the value for this marker decreases
<figure><img src="image_46.png" /><figcaption>Figure 4.9: Difference between the treatment arm model and the simulated placebo arm model. Each line is a separate individual. The plots show that the treatment seems to improve the dopaminergic transporter signal over time.</figcaption></figure>

<table><caption>Table 4.6: Prediction metrics for the regression models computed with four folds. $R^2$: R-squared coefficient of maximum a posteriori (MAP) prediction. MAE: Mean Absolute Error of MAP prediction. LL: posterior log-likelihood.</caption><thead><tr><th>Outcome</th><th>Regression</th><th>R<sup>2</sup></th><th>MAE</th><th>LL</th><th>Runtime (s)</th></tr></thead><tbody><tr><td rowspan="3">Putamen<br>Ipsilateral</td><td>BLR</td><td>0.61±0.02</td><td>0.053±0.003</td><td>-1.48±0.30×10<sup>3</sup></td><td>37±3</td></tr><tr><td>BKMR</td><td>0.72±0.04</td><td>0.047±0.003</td><td>-39±13</td><td>118±23</td></tr><tr><td>BART</td><td>0.77±0.02</td><td>0.044±0.003</td><td>-14.5±4.4</td><td>3.2±0.2×10<sup>3</sup></td></tr><tr><td rowspan="3">Putamen<br>Contralateral</td><td>BLR</td><td>0.69±0.03</td><td>0.061±0.006</td><td>-1.5±0.3×10<sup>3</sup></td><td>53±20</td></tr><tr><td>BKMR</td><td>0.77±0.04</td><td>0.050±0.005</td><td>-61±22</td><td>192±72</td></tr><tr><td>BART</td><td>0.82±0.02</td><td>0.046±0.005</td><td>-16±3</td><td>1.7±0.2×10<sup>3</sup></td></tr></tbody></table>

with time. On the second row, we see that the treatment effect seems to be correlated with the cumulative treatment dose, although the credibility interval is large due to the averaging over all the over covariates effect. However, it is much more interesting to look at the mean treatment effect as a function of both these two covariates as their effect is dependent on each other. The BART and BKRM in the last row show that at early disease stages (Putamen values around 1) the cumulative treatment dose has almost no influence as the treatment effect is close to 0. However, at later stages (when the putamen values are closer to 0.3) the treatment effect is correlated with the cumulative treatment dose. Overall the regression models suggest that the dopaminergic treatment is reducing the loss of dopaminergic transporters in the putamen compared to the natural disease progression. In practice, it might imply that the treatment, although it is thought to be purely symptomatic, might have a slight disease-modifying effect.
<figure><img src="image_47.png" /><figcaption>Figure 4.10: Variable importance computed with BART on full set of covariates. Variables predicted by the disease progression model are noted with "_pred". MDSX_TOTAL: Movement Disorder Society - Unified Parkinson's Disease Rating Scale, part X. TD: Tremor Dominant. PIGD: Postural Instability and Gait Disorder. PD_MED: medication category. LEDD: Levodopa equivalent daily dose (=treatment dose). AGE_SINCE_BL: age since baseline. BMI: Body Mass Index.</figcaption></figure>
<figure><img src="image_48.png" /><figcaption>Figure 4.11: Average treatment effect with 95% credibility interval. Top row is the average effect conditioned on predicted putamen value. Treatment effect is stronger as the putamen loses dopamine transporters. Middle row is the average effect conditioned on the cumulative treatment dose. Bottom row shows the interaction with both. At low putamen values the treatment effect correlates with the cumulative treatment dose.</figcaption></figure>
4.2.5 Discussion

We proposed a fully Bayesian method to infer treatment effects in observational cohorts without a control arm. Our approach relies on a disease progression model to infer the unobserved disease state of a patient under no treatment. The treatment effect is then modeled with a regression using various predictors, including disease markers predicted by the disease progression model. A subject disease state can then be predicted by adding the treatment effect to the latent disease progression.

In the simulation experiment we showcased the strength of the method while highlighting the main defects, namely the learning of the latent disease progression model's bias by the treatment effect model. This does not impact predictions of the patients under treatment, but this negates the main interest of this method which is to disentangle the conditional treatment effect *CATE* from the natural disease history.

We then showed that the model proved to be reliable in a real-life setting, with close to ground-truth data for ON and OFF states. Our approach also outperformed the baseline approach in predicting a patient's disease state in the second experiment.

In the last experiments, we extended the method to scores that were not observed without treatment effect, limiting our validation possibilities. We still obtained high predicting performances on the MDS-UPDRS part II (quality of life score) and on imaging data. More specifically, we applied our method to SPECT data for dopamine transporters in the striatum in Parkinson's disease. In this case, our results suggest that the dopaminergic treatment affects the disease progression. There is an indication that the cumulative treatment dose slows down the loss of dopamine transporters in the putamen.

The treatment effect model can be useful for understanding and optimizing treatment strategies for new patients. Another possible application concerns new treatments, as many are being tested for neurodegenerative diseases, with high expectations. The standard way to validate a new treatment is a clinical trial. However, a trial able to demonstrate a meaningful effect in a very slowly progressive disease is costly, and the statistical power is often very weak. To reduce the size of the trial one could leverage a latent disease progression model to create a synthetic control arm and estimate the treatment effect as we did in the observational study.

Our approach has inherent limits based on the latent disease progression model quality. The better the predictions for the latent disease progression, the better the treatment effect estimation. A proper theoretical study about the properties of this estimator remains to be done, especially if we wish this method to be applied to one-arm clinical trials with the disease progression model providing a synthetic control arm.

Finally we have seen that this approach is suited for simple treatment effects as symptomatic ones, but including disease-modifying effects requires to add time-dependent covariates with history information (cumulative treatment dose for instance). This is not the best way to model such an effect, so we will focus on a new model, more adequate for disease-modifying treatment effects.
## 4.3 Piecewise-geodesic model

We have seen a "simple" additive model, based on the hypothesis that a treatment is purely symptomatic. This allowed us to avoid any complex longitudinal modelling. Let's now use the ODE formulation in order to understand what is at stake with disease-modifying treatments:

$$ \frac{\partial \mathbf{y}}{\partial t}(t) = f(\mathbf{y}(t)) \qquad (4.3.1) $$

$f$ is the function we want to learn. It is possible to parametrize it as was done with the penalty method in the second experiment previously. It is notable that the introduction of covariates (such as the treatment $T$) in the ODE equation is not as easy as with mixed-effect models. In these mechanical models it is quite complex to introduce time-dependent covariates since the equations' integration won't be tractable. However we can choose a certain form for these time-dependent covariates in order to keep the properties of the ODE. Here we postulate that the treatment effect is a function of the position $\mathbf{y}$ and the treatment dosage $T(t)$ (which is supposed to be known at every timepoint). This echoes the work in the previous model where we saw that the most important variables for treatment effect models were the treatment dose and the latent disease state. The ODE equation becomes:

$$ \frac{\partial \mathbf{y}}{\partial t}(t) = f(\mathbf{y}(t)) + g(\mathbf{y}(t), T(t)) \qquad (4.3.2) $$

where $g$ is the treatment effect such that $g(\cdot, 0) = 0$. Note that if $T(t)$ is not fully known, for instance if we don't have prior knowledge about how treatment will be prescribed in future visits, or if $T$ is a control process which is latent but not directly observed, then the above equation cannot be solved. We would require further hypothesis about $T$, e.g a template curve or an interpolation method in between visits.

Another important remark: for symptomatic effects we can assume that $||g|| \gg ||f||$, and we have convergence to a stable point for the ODE $\frac{\partial \mathbf{y}}{\partial t}(t) = c+g(\mathbf{y}(t), T(t))$. Let me illustrate in the case of Parkinson's disease: L-dopa does not act instantly on the organism, it operates within one hour approximately, during which the dynamic of the symptoms changes until it stabilizes to a certain point. However one hour can be considered very small compared to the disease progression scale, which can span up to decades. Therefore the effect of the L-dopa as an "additive" symptomatic treatment can be seen as time-limit approximation of a very fast dynamic effect. However it is often not easy to learn on different time scales, therefore if the treatment has a symptomatic and a disease-modifying effect then it is easier to learn the two parts separately rather than the total effect on the dynamics directly.

### 4.3.1 Method

For this section it is preferable to understand the basics of Riemannian geometry. If needed, the reader can refer to the introduction to Riemannian geometry in the appendix 15.

The classical way to handle treatment in statistical models is to treat it as a time-dependent covariate. If we go back to a general formulation of a mixed-effect model, we write:

$$ \mathbf{y}_{ij} = \eta(t_{ij}, \mathbf{X}_{ij}; \mathbf{z}_{pop}, \mathbf{z}_i) + \epsilon_{ij} \quad [\text{Mixed-effect model with time-dependent covariates}] \qquad (4.3.3) $$

$$ \mathbf{z}_{ij} = f(T_{ij}, \mathbf{z}_{ij'}, j' < j) \qquad [\text{Covariate structure}] \qquad (4.3.4) $$
The main idea is that the time-dependent covariate $z_{ij}$ linked to the treatment is a function of the current dose but also of the history of the treatment. If we only used a function of the treatment dose, we would end up with a symptomatic model since when the treatment has a null dose the function $f$ with no history needs to be null too (no effect when no treatment).

All of the modelling task resides in choosing how to integrate previous treatment history into those covariates. Here we will specify our proposed approach with the DCM model in its generic framework. In this context we will need to specify $T$ over time as well. Indeed we want to avoid any integration as in a dynamic approach, so we will assume that $T$ is piece-wise constant. This is a fair hypothesis since at each visit the patient will be given a treatment prescription with a certain dose which should be constant until the next visit. The only drawback of this hypothesis is that it stops us from making future predictions without a new treatment policy function. Writing this hypothesis for one individual $i$ yields:

$$ T_i(t) = \sum_{j=1}^{n_i} a_{ij} \mathbb{1}_{t \in [b_{i,j}, b_{i,j+1}]} \qquad (4.3.5) $$

with $0 < b_{i,1} < ... < b_{i,j} < b_{i,j+1} < ... < b_{i,n_i+1} = \infty$ being the dates of the changes in treatment, called breakpoints. Note that $T = 0$ before $b_{i,1}$. Let's write the set of individual parameters in the absence of treatment (or before $b_{i,1}$) as $(\tau_i^0, \xi_i^0, \mathbf{w}_i^0)$. The main idea of our model is that on each interval $[b_{i,j}, b_{i,j+1}]$ the trajectory will follow a geodesic in the manifold (or more precisely be an Exp-parallelization of a geodesic). Hence the name: Piece-wise Geodesic model.

We will build a reference curve that is piece-wise geodesic, such that the trajectory of individual $i$ is obtained by Exp-parallelization 10 from it, as in the DCM. However, since each individual has now its own treatment function $T$ with different breakpoints, the reference curve will be different for each individual. We will denote $\gamma_i$ the reference curve for individual $i$. We will use the ODE characterization of the geodesics. Therefore we have $t_0, \mathbf{p}_0, \mathbf{v}_0$ as the parameters defining the geodesic $\gamma_0$. Before $b_{i,1}$, the Piece-wise Geodesic $\gamma_i$ is equal to $\gamma_0$. At $b_{i,1}$, the triplet $(b_{i,1}, \gamma_0(b_{i,1}), \dot{\gamma}_0(b_{i,1}))$ defines the same geodesic $\gamma_0$. The treatment applies a perturbation on the trajectory by changing the value of the derivative only, yielding a new triplet $(b_{i,1}, \gamma_0(b_{i,1}), \mathbf{v}_1)$ which defines a new geodesic passing by $\gamma_0(b_{i,1})$. The reference curve $\gamma_i$ follows this new geodesic until the next breakpoint $b_{i,2}$, where the same process is reiterated.

More formally: we fix $j \in [|1, n_i|]$. The piece of geodesic followed during the interval $(b_{i,j}, b_{i,j+1})$ is parametrized by $(b_{i,j}, \mathbf{p}_{i,j}, \mathbf{v}_{i,j})$ where $\gamma_i(b_{i,j}^+) = \mathbf{p}_{i,j}$ and $\dot{\gamma}_i(b_{i,j}^+) = \mathbf{v}_{i,j}$. At time $b_{i,j+1}, \gamma_i(b_{i,j+1}^-) = \mathbf{p}_{i,j+1}$ (we assume continuity of the trajectory) and $\dot{\gamma}_i(b_{i,j+1}^-) = \tilde{\mathbf{v}}_{i,j}$. The perturbation is encoded as a function $f$ such that:

$$ \dot{\gamma}_i(b_{i,j+1}^+) = \mathbf{v}_{i,j+1} = f(\mathbf{p}_{i,j+1}, \tilde{\mathbf{v}}_{i,j}, T_i(b_{i,j}^-), T_i(b_{i,j}^+)) \qquad (4.3.6) $$

We call $f$ the treatment function. We could also include covariates in this function, such as the weight (patients with a high weight often need a higher dose for a similar treatment effect) or treatment responsiveness, but we will keep it simple for now.

The individual trajectory then results from an Exp-parallelization 10 of the reference curve $\gamma_i$ with the space-shift $\mathbf{w}_i$. Exp-parallelisation means that the model will be "locally" following the same modification as $\gamma_i$ upon the curve $\eta^{\mathbf{w}_i}$. It is well defined thanks to the continuity of $\gamma_i$ and the fact that the Exp-parallelization is well defined on each interval. Figure 4.12 illustrates this procedure.
<figure><img src="image_49.png" /><figcaption>Figure 4.12: Left: Piewe-wise geodesic model scheme. $\gamma_i$ is the reference curve for individual i and $\eta_i^{w_i}$ is the Exp-parallelization of $\gamma_i$ with $w_i$. Right: samples of data generated following a Piece-wise Geodesic model on a logistic curves manifold. The black vertical lines correspond to the breakpoints, the individual trajectory after Exp-parallelization is shown in red and data points are in blue.</figcaption></figure>

One could argue that the model could be discontinuous at the breakpoints, and that we could learn the perturbation on the position as well. This is very sensible, as some treatments may have radical changes on the disease trajectory, for instance a surgery operation. However it becomes tricky to define the Exp-parallelization in a discontinuous setting: how do we transport $w_i$ from $\gamma_i(b_{i,j}^-)$ to $\gamma_i(b_{i,j}^+)$ ? We could do the transport along the geodesic between these two points, which would require to compute the Riemannian logarithm. This would lead to heavy computations and a lot of new parameters, which would be too much considering that the current choice of model already has issues with the estimation.

**Estimation** With this choice of model, treatment effect understood as the parametrization of *f* is a population parameter in the mixed-effect model. Estimation of the parameters of *f* is performed during the MCMC-SAEM algorithm, as described in the next chapter. The method has been implemented in Python, under the Leapy library.
### 4.3.2 Experiments

We only performed simulation experiments with this model for two reasons: the model is very data-hungry which is not a problem in simulations, and we have not yet a large neurodegenerative dataset with a disease-modifying treatment. We show that the model is able to recover the ground truth in a simplified setting.

We want to parametrize the treatment effect function $f$ as a function of only the position $p$ for simplification. In order to still keep a dependence on the velocity, we will write the perturbation as a function $g$ with a small norm on top of the identity: $f(\mathbf{p}, \mathbf{v}) = \mathbf{v} + g(\mathbf{p})$.

We then need to parametrize $g$. As often in this thesis, we resorted to a Reproducing Kernel Hilbert Space (RKHS) with kernel $k$, but one could use a simple linear function if needed. Using the representer theorem, the function $g$ writes:

$$g(\mathbf{p}) = \sum_{ij} \omega_{i,j} k(\mathbf{p}, \mathbf{p}_{i,j}) \qquad (4.3.7)$$

where the parameters are the weights $\omega_{i,j}$. They are as numerous as there are breakpoints in the dataset, which is very variable depending on the cohort. In a clinical trial for instance, we would only have one breakpoint per patient corresponding to the beginning of the treatment.

Data is generated by randomly sampling individual parameters, then for each individual we generate breakpoints and compute the trajectory. Finally we generate random timepoints for the visits, take the trajectory's value at the visits and add a small noise. We will adjust the settings with the number of breakpoints and the number of visits per individual to evaluate the limits of the learning capabilities of our model.

**Univariate setting** In a first setting, we generated each visit to be a breakpoint. Intuitively, the changes in $\mathbf{v}$ depending on $\mathbf{p}$ can be understood as the practitioner adjusting the treatment of the patient at each visit $b_{i,j}$ based on the current observed state $\mathbf{p}_{i,j}$. We now generate data using the Euclidean setting, where the geodesics are straight lines. This simplified setting amounts to the model being able to capture the change of slope at each breakpoint. The generating function for the change of slope is a heavyside. We chose the Gaussian kernel for the learning of $g$, which means that $g$ will be a linear combination of radial basis functions. We optimized the kernel bandwidth by random search.

Figure 4.13 shows the results of the simulation of 1024 individuals, 10 visits per individual, and the model obtained after convergence. We observe that the ground truth function is quite accurately modelled, the wavelets owning to the Gaussian kernel.

We repeated the experiment with other generating functions for $g$, such as linear or sinusoidal functions. In a setting with 1024 individuals and 10 visits per individual, the model is able to learn the target treatment function quite quickly, as shown in Figure 4.14. We also performed experiments with only one breakpoint per patient. In these experiments the ability of the model to learn the proper treatment function scaled with the number of individuals, i.e. the number of breakpoints, rather than the number of visits.
<figure><img src="image_50.png" /><figcaption>Figure 4.13: Piece-wise Geodesic model results on simulated data. The results are shown after 300 iterations of MCMC-SAEM, with an optimization of the weights $\omega_{i,j}$ every 10 iterations. The first row corresponds to three individuals, where the true trajectory is shown in red with the observations in blue, 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.</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>Figure 4.14: Convergence of the learned function $g$ during the iterations of the MCMC-SAEM. The results are shown after 10, 100 and 1000 iterations of MCMC-SAEM, with an optimization of the weights $\omega_{i,j}$ every 10 iterations.</figcaption></figure>