Bayesian latent time joint mixed effect models for multicohort longitudinal data

TL;DR

This paper introduces a Bayesian latent time joint mixed effects model to analyze long-term neurodegenerative disease progression from short-term multicohort data, enabling better understanding and intervention strategies.

Contribution

It presents a novel latent time modeling approach with MCMC estimation for long-term disease dynamics using short-term multicohort data.

Findings

Effective in characterizing disease progression

Accurate in estimating long-term disease trajectories

Applicable to Alzheimer's disease data

Abstract

Characterization of long-term disease dynamics, from disease-free to end-stage, is integral to understanding the course of neurodegenerative diseases such as Parkinson's and Alzheimer's; and ultimately, how best to intervene. Natural history studies typically recruit multiple cohorts at different stages of disease and follow them longitudinally for a relatively short period of time. We propose a latent time joint mixed effects model to characterize long-term disease dynamics using this short-term data. Markov chain Monte Carlo methods are proposed for estimation, model selection, and inference. We apply the model to detailed simulation studies and data from the Alzheimer's Disease Neuroimaging Initiative.

Figures

Tables

Table 1. Table 1 : Simulation study results. One hundred simulated data sets were generated under model M1 (univariate Gaussian random effects) and model M2 (multivariate Gaussian random effects). Both models were fit to all 200 data sets.

	True Model: M1						True Model: M2
	Fitted Model: M1			Fitted Model: M2			Fitted Model: M1			Fitted Model: M2
Parameter	Bias	MSPE	$C_{95}$	Bias	MSPE	$C_{95}$	Bias	MSPE	$C_{95}$	Bias	MSPE	$C_{95}$
${\beta }_1$	-0.0022	0.0086	0.96	-0.0005	0.0087	0.96	-0.0013	0.0064	0.96	-0.0016	0.0062	0.98
${\beta }_2$	-0.0043	0.0131	0.93	0.0003	0.0131	0.95	0.0006	0.0133	0.94	0.0020	0.0131	0.94
${\beta }_3$	-0.0059	0.0157	0.97	-0.0061	0.0160	0.97	-0.0041	0.0155	0.97	-0.0044	0.0147	0.97
${\beta }_4$	-0.0425	0.0481	0.97	-0.0409	0.0477	0.97	-0.0428	0.0550	0.95	-0.0420	0.0556	0.97
${\gamma }_1$	0.0274	0.0025	0.93	0.0321	0.0027	0.90	0.0305	0.0036	0.83	0.0449	0.0044	0.89
${\gamma }_2$	0.0194	0.0016	0.96	0.0222	0.0018	0.95	0.0246	0.0025	0.93	0.0360	0.0033	0.93
${\gamma }_3$	0.0265	0.0031	0.96	0.0324	0.0033	0.97	0.0388	0.0061	0.83	0.0551	0.0067	0.92
${\gamma }_4$	0.0398	0.0082	0.96	0.0504	0.0087	0.95	0.0529	0.0134	0.85	0.0621	0.0102	0.95
${\sigma }_1$	0.0006		0.91	0.0006		0.92	0.0006		0.92	0.0007		0.91
${\sigma }_2$	0.0020		0.95	0.0021		0.94	0.0009		0.97	0.0020		0.94
${\sigma }_3$	0.0023	0.0002	0.94	0.0028	0.0002	0.93	-0.0001	0.0002	0.93	0.0021	0.0002	0.94
${\sigma }_4$	0.0018	0.0002	0.96	0.0022	0.0002	0.96	0.0028	0.0002	0.96	0.0021	0.0002	0.96
${\sigma }_{\delta }$	-0.2383	0.4542	0.93	-0.3312	0.4846	0.91	-0.2502	0.7656	0.84	-0.4459	0.6715	0.89
${\sigma }_{01}$	0.0142	0.0037	0.93	0.0336	0.0048	0.87	0.0162	0.0021	0.91	0.0185	0.0022	0.89
${\sigma }_{02}$	0.0061	0.0052	0.97	0.0313	0.0064	0.95	0.0298	0.0057	0.96	0.0258	0.0055	0.95
${\sigma }_{03}$	0.0033	0.0060	0.94	0.0255	0.0067	0.94	0.0262	0.0058	0.91	0.0147	0.0051	0.96
${\sigma }_{11}$	0.0078	0.0066	0.94	0.0333	0.0079	0.91	0.0067	0.0066	0.89	0.0093	0.0061	0.93
${\sigma }_{12}$	0.0114	0.0095	0.97	0.0462	0.0119	0.94	0.0122	0.0206	0.94	-0.0089	0.0183	0.95
${\sigma }_{13}$	0.0003	0.0088	0.95	0.0297	0.0101	0.92	-0.0014	0.0111	0.97	-0.0283	0.0104	0.98
${\sigma }_{14}$	0.0105	0.0053	0.96	0.0356	0.0067	0.93	0.0114	0.0064	0.97	0.0076	0.0056	0.97
% best LOOIC	59%			41%			4%			96%
% best WAIC	65%			35%			19%			81%
Quartiles $(Q_1,Q_2,Q_3)$ of the criteria differences between M2 and M1
LOOIC ${\mathrm{diff}}_{\mathrm{M2}-\mathrm{M1}}$	$(-4.7313,2.2414,9.6980)$						$(-36.2904,-24.3240,-16.4950)$
WAIC ${\mathrm{diff}}_{\mathrm{M2}-\mathrm{M1}}$	$(-2.0153,2.5876,7.4679)$						$(-12.0306,-6.6466,-1.2454)$

Table 2. Table 2 : Results of posterior estimate of parameters for the proposed LTJMM fit to seven outcomes from the Alzheimer’s Disease Neuroimaging Initiative (ADNI).

Parameter	Posterior	95%	Parameter	Posterior	95%
	Mean	Credible Interval		Mean	Credible Interval
Hippocampal volume			CSF Amyloid Beta 1-42
Intercept	-4.102	$(-4.516,-3.688)$	Intercept	-1.454	$(-2.061,-0.872)$
Age	0.056	$(0.052,0.06)$	Age	0.016	$(0.011,0.022)$
APOE$\epsilon$4	0.317	$(0.258,0.394)$	APOE$\epsilon$4	0.802	$(0.704,0.888)$
Sex	-0.288	$(-0.364,-0.215)$	Sex	-0.028	$(-0.116,0.049)$
Education	-0.003	$(-0.015,0.008)$	Education	-0.012	$(-0.03,0.005)$
Latent time	0.035	$(0.03,0.04)$	Latent time	0.035	$(0.03,0.039)$
${\sigma }_1$	0.137	$(0.134,0.141)$	${\sigma }_5$	0.021	$(0.019,0.022)$
Alzheimer’s Disease Assessment Scale			Functional Activities Questionnaire
Intercept	-1.089	$(-1.57,-0.646)$	Intercept	-3.012	$(-3.934,-2.012)$
Age	0.025	$(0.02,0.03)$	Age	0.035	$(0.025,0.045)$
APOE$\epsilon$4	0.434	$(0.355,0.514)$	APOE$\epsilon$4	0.861	$(0.678,1.021)$
Sex	-0.246	$(-0.314,-0.168)$	Sex	-0.335	$(-0.471,-0.196)$
Education	-0.051	$(-0.064,-0.038)$	Education	-0.061	$(-0.089,-0.037)$
Latent time	0.045	$(0.039,0.05)$	Latent time	0.096	$(0.085,0.107)$
${\sigma }_2$	0.327	$(0.32,0.333)$	${\sigma }_6$	0.842	$(0.826,0.858)$
FDG PET			CSF Tau
Intercept	-1.820	$(-2.348,-1.3)$	Intercept	-0.766	$(-1.265,-0.202)$
Age	0.028	$(0.023,0.034)$	Age	0.012	$(0.007,0.017)$
APOE$\epsilon$4	0.453	$(0.363,0.54)$	APOE$\epsilon$4	0.621	$(0.505,0.72)$
Sex	-0.181	$(-0.259,-0.103)$	Sex	0.073	$(-0.036,0.19)$
Education	-0.026	$(-0.041,-0.01)$	Education	-0.029	$(-0.048,-0.007)$
Latent time	0.042	$(0.037,0.048)$	Latent time	0.032	$(0.028,0.036)$
${\sigma }_3$	0.294	$(0.283,0.304)$	${\sigma }_7$	0.016	$(0.014,0.017)$
Amyloid PET			Standard deviation of latent time
Intercept	-0.685	$(-1.328,-0.013)$	${\sigma }_{\delta }$	14.262	$(12.684,16.105)$
Age	0.004	$(-0.002,0.011)$
APOE$\epsilon$4	0.785	$(0.669,0.912)$
Sex	0.159	$(0.067,0.247)$
Education	-0.007	$(-0.029,0.014)$
Latent time	0.035	$(0.03,0.04)$
${\sigma }_4$	0.289	$(0.273,0.307)$

Table 3. Table 3 : Results of model selection for models fit to seven outcomes from the Alzheimer’s Disease Neuroimaging Initiative (ADNI).

Model	WAIC	LOOIC	DIC
${\text{LTJMM}}_1$	10246.31	14156.73	11070.65
${\text{LTJMM}}_2$	9684.75	13482.41	10528.48
MM	10023.88	13938.55	11010.95
JMM	9308.89	12385.48	10170.97
MM (diagnosis)	9792.05	13778.70	10761.32
JMM (diagnosis)	9420.06	13082.31	10286.82

Equations

$$h_{k}\left(\theta_{ijk}\right)=\eta_{ijk}=\textbf{x}^{\prime}_{t_{ijk}}\mbox{\boldmath$\beta$}_{k}+\gamma_{k}\left(t_{ijk}+\delta_{i}\right)+\alpha_{0ik}+\alpha_{1ik}t_{ijk}+\varepsilon_{ijk},$$

$$y_{ijk}=\textbf{x}^{\prime}_{ijk}\left(t_{ijk}\right)\mbox{\boldmath$\beta$}_{k}+\gamma_{k}\left(t_{ijk}+\delta_{i}\right)+\alpha_{0ik}+\alpha_{1ik}t_{ijk}+\varepsilon_{ijk},$$

$$y_{ijk}=\textbf{x}^{\prime}_{ijk}\left(t_{ijk}\right)\mbox{\boldmath$\beta$}_{k}+\gamma_{k}t_{ijk}+\tilde{\alpha}_{0ik}+\alpha_{1ik}t_{ijk}+\varepsilon_{ijk},$$

$${\tilde{\alpha }}_{0ik}$$

$$0$$

Full text

1 Introduction

Disentangling the effects of normal aging from the effects of pathophysiology on neurobiological trajectories is crucial for predicting who will develop Alzheimer’s or other neurodegenerative diseases. The determination of neurobiological trajectories is typically based on data obtained over a relatively short time frame relative to the time span of normal aging and neurodegeneration [1]. There are four key barriers to the accurate characterization of longitudinal trajectories of neurodegeneration and pathology in older adults: 1) late-life trajectories can span 40 or more years; 2) there is substantial variation in the timing and temporal dynamics of late-life trajectories; 3) study selection criteria and sampling have a substantially different impact across age groups; and 4) there are potentially significant levels of censoring due to death and disability in late-life that cannot be treated as independent of pathology trajectories in most cases. Accounting for these issues is crucial to studying the onset, course and fluctuations of neurodegenerative diseases [2, 3], and for guiding when and how to target interventions and preventative approaches [4]. A better understanding of how and when to intervene necessitates the determination of dynamic networks of trajectories of markers of cognition, function, and pathology. Long-term disease dynamics are of great interest and importance, and have been hypothesized without rigorous methods.

Methods for longitudinal data analysis are ubiquitous. For example, generalized linear or nonlinear mixed-effects models [5, 6] are commonly applied to model growth curves of height or weight over time since a meaningful “time zero” (e.g. birth). However, neurodegenerative diseases like Alzheimer’s and Parkinson’s that progress over long periods of time are typically studied by taking longitudinal samples of multiple cohorts at different stages of disease [7, 8]. The multiple cohorts have no common meaningful time zero. Time since dementia diagnosis might be a candidate time zero, however the transition to dementia is censored for many individuals in the pre-dementia cohorts. Also, the dementia diagnosis is subjective and may vary from one clinician to the next in a multicenter study. Further, dementia is not an absorbing state and reversion to pre-dementia diagnoses occur.

We propose a latent time joint mixed effects model (LTJMM) for characterizing biomarker trajectories in aging. The model extends joint mixed effects models [5] to include an individual-specific latent time shift. Each individual’s latent time shift is shared across all of their outcomes and represents the extent of their long-term disease progression. The model is similar to others proposed for analysis of the Alzheimer’s Disease Neuroimaging Initiative (ADNI)[7] (e.g. Jedynak et al.[9] and Donohue et al.[10]). However, the proposed model can also accommodate covariates for fixed-effects and the Bayesian implementation allows flexible but rigorous interrogation of the posterior distribution to make inferences about long-term disease dynamics and the potential propagation of treatment effects from early biomarkers to downstream cognitive and functional measures. Estimation and inference are accomplished using Hamiltonian Monte Carlo sampling as implemented in Stan [11] and the R package rstan (Stan Development Team[12], version 2.10.1). The Stan code model specifications and related tools for estimation and prediction are available as an additional R package from https://bitbucket.org/mdonohue/ltjmm.

2 Methods

2.1 Latent time mixed-effects models for joint longitudinal responses

Suppose that $p$ response variables are measured for $n$ individuals at different follow-up times. Follow-up may differ from subject to subject and outcome to outcome, so we denote the measured outcome $k$ for individual $i$ at time $j$ as $y_{ijk}$, where $i=1,\ldots,n$, $k=1,...,p$ and $j=1,\ldots,q_{ik}$. The $p$ measured outcomes could be mixtures of binary, count, ordinal and continuous. Assuming $y_{ijk}$ has a distribution in the exponential family with canonical parameter $\theta_{ijk}$, the latent time joint mixed effects model (LTJMM) is given by

[TABLE]

where $h_{k}\left(\cdot\right)$ is a monotonic differentiable link function specific to the $k$-th outcome, $\eta_{ijk}$ is the linear predictor for the $k$-th outcome from individual $i$ at time $j$. The vector $\textbf{x}^{\prime}_{t_{ijk}}$ represents possibly time-varying covariates, and the vector \mbox{\boldmath\beta}_{k} is their corresponding regression coefficients. The $\varepsilon_{ijk}\sim\mathcal{N}\left(0,\sigma_{k}^{2}\right)$ is a measurement error term, which accounts for outcome-specific variance.

The parameters $\alpha_{0ik}$ and $\alpha_{1ik}$ are the subject and outcome specific random intercept and slope, and $\delta_{i}$ is the subject-specific time shift shared among outcomes. “Short-term” observation time is represented by observed covariate $t_{ijk}$. The parameter $\gamma_{k}$ corresponds to the outcome-specific slope with respect to “shifted” or “long-term” time $t_{ijk}+\delta_{i}$. The time shift $\delta_{i}$ quantifies the progression of the $i$-th individual relative to the population, which is assumed to follow $\delta_{i}\sim\mathcal{N}\left(0,\sigma_{\delta}^{2}\right)$. And $\alpha_{1ik}$ provides the information on whether individual $i$ is evolving faster or slower than the average individual for outcome $k$.

When applying the model, we typically include a fixed effect for age, as a time-varying covariate, for each outcome. This admits two relevant rate-of-change parameters: one for biological age and one for “disease time”. This reflects the fact that individuals can reach the same stage of disease at different ages, and allows independent effects of healthy aging and disease progression. The random effects accommodate additional subject to subject variation in the level and rate of change of disease markers. The link function, $h_{k}$, can be used to accommodate any distribution from the exponential family.

A constraint must be placed on (1) to ensure model identifiability based on the observed data without relying on informative prior distributions. In particular, to ensure identifiability of $\delta$, the following constraint on random intercepts is sufficient: $\sum_{k=1}^{p}\alpha_{0ik}=0$ for all individuals $i=1,\ldots,n$. A proof of identifiability is available in the supplemental material. Model (1) can be modified to address more flexible nonlinear relationships over time by incorporating higher-order terms or splines [13], provided these extensions maintain monotonicity. Shape invariant models and self modeling regression (e.g. Kneip and Gasser[14]) similarly require positive first derivatives. Under model (1), we trivially assume that at least one of the outcome-specific slopes $\gamma_{k}$ is not equal to zero, $k=1,\ldots,p$. In our application, all the outcomes are oriented to be increasing, and thus we assume all $\gamma_{k}>0$, $k=1,\ldots,p$.

Finally, we explore two different distribution assumptions for random effects, namely:

1. Univariate Gaussian random effects: $\alpha_{0ik}\sim\mathcal{N}\left(0,\sigma_{0k}^{2}\right)$ and $\alpha_{1ik}\sim\mathcal{N}\left(0,\sigma_{1k}^{2}\right)$;

2. Multivariate Gaussian random effects: \mbox{\boldmath\alpha}_{i}\sim\mathcal{N}\left(\textbf{0},\mbox{\boldmath\Sigma}\right),

where $\Sigma$ is a covariance matrix of dimension $2p$. The latter assumption allows more direct exploration and inference regarding the correlation of response variables.

2.2 Prior specification

Since an improper prior may result in an improper posterior, we will use proper but weakly informative priors on all the model parameters [15]. The regression parameters $\beta$ and $\gamma$ are assigned independent weakly informative normal $\mathcal{N}\left(0,100\right)$ priors (truncated below by 0 for $\gamma$). Without specific information, the half-Cauchy prior is a good default choice for scale parameters [15]. The standard deviations, i.e., $\sigma_{\delta}$, $\sigma_{01},\cdots,\sigma_{0p}$, $\sigma_{11},\cdots,\sigma_{1p}$ and $\sigma$, are given weakly informative half-Cauchy priors with a small scale, i.e., $\text{half-Cauchy}\left(0,2.5\right)$. The prior variance 100 for the regression parameters is chosen to be sufficiently high to be vague enough, but sufficiently low to avoid slow mixing due to near impropriety of the posterior [16].

To ensure efficiency and arithmetic stability, we applied the Cholesky decomposition to the covariance matrix for the random effects, $\Sigma$, allowing us to model the standard deviations and correlations independently. We first decomposed the covariance matrix as \mbox{\boldmath\Sigma}=\mbox{\boldmath\Lambda}\mbox{\boldmath\Omega}\mbox{\boldmath\Lambda}, where $\Lambda$ is a diagonal matrix with diagonal elements $\sigma_{01},\cdots,\sigma_{0p}$, $\sigma_{11},\cdots,\sigma_{1p}$, and $\Omega$ is the correlation matrix with 1’s on the diagonal and off-diagonal elements $\rho$. A Choleksy decomposition gives \mbox{\boldmath\Omega}=\boldsymbol{L}\boldsymbol{L}^{\prime}, where $\boldsymbol{L}$ is a lower triangular matrix. Furthermore, \mbox{\boldmath\Sigma}=\mbox{\boldmath\Lambda}\boldsymbol{L}\boldsymbol{L}^{\prime}\mbox{\boldmath\Lambda}. During sampling, a draw is obtained from a multivariate Gaussian $\textbf{z}\sim\mathcal{N}\left(\textbf{0},\textbf{I}\right)$, and then the random effects are calculated as \mbox{\boldmath\Lambda}\boldsymbol{L}\textbf{z}. Here, z is regarded as a random reparameterization and independent of $\rho$. As suggested in Stan [17] manual, we imposed an LKJ prior on the correlation matrix [18].

2.3 Model comparison criteria

We use two model comparison criteria, namely, the widely applicable information criterion (WAIC) [19] and the leave-one-out cross-validation information criterion (LOOIC) [20, 21]. WAIC and LOOIC can be computed using the log-likelihood evaluated at the posterior simulations of the parameter values. Both have various advantages over the deviance information criterion (DIC). WAIC can be viewed as an improvement on DIC for Bayesian methods [22]. We refer to Gelman et al.[23] for the detailed rationale for the preference of WAIC to DIC. As WAIC is reported on the deviance scale [23], a difference in WAIC value of $2-6$ is considered as positive evidence, $6-10$ strong evidence, and $>10$ very strong evidence [24].

Approximate LOO can be computed using raw importance sampling (IS, Gelfand et al.[20]). Vehtari et al.[22] improved the LOO estimate using Pareto smoothed importance sampling (PSIS) method which provides a more accurate and reliable estimate than IS. The R package loo [25] provides tools for efficient computation of WAIC and LOOIC. The minimum WAIC and LOOIC indicate the best fit.

3 Results

3.1 Simulation study

We conducted a simulation study to assess the model performance and assess small sample bias. Data were generated from two candidate LTJMMs with identity link: (M1) Univariate Gaussian random effects; (M2) Multivariate Gaussian random effects. We simulated $n=100$ individuals with $p=4$ outcomes and $q=4$ time points each. The observation times for each individual were sampled from a Uniform distribution $t_{ijk}\sim\text{Uniform}\left(0,10\right)$. Additional simulation parameters were set to \mbox{\boldmath\beta}=\left(1,0.5,2,0.8\right)^{\prime}, \mbox{\boldmath\gamma}=\left(0.2,0.1,0.25,0.5\right)^{\prime}, $\sigma_{\delta}=4$, \mbox{\boldmath\sigma}=\left(0.1,0.2,0.3,0.25\right)^{\prime}, \mbox{\boldmath\sigma}_{\alpha_{0}}=\left(0.5,1,0.8\right)^{\prime}, and \mbox{\boldmath\sigma}_{\alpha_{1}}=\left(1,2,1.5,1\right)^{\prime}.

Let $\hat{\theta}_{m}$ denote the posterior estimate of a parameter $\theta$ for the $m$-th simulation and let $M$ denote the total number of simulations. We consider the following quantities to assess model performance: (i) total bias, $\text{Bias}=\sum_{m=1}^{M}\left(\hat{\theta}_{m}-\theta_{m}\right)/M$, (ii) total mean squared error of prediction, $\text{MSPE}=\sum_{m=1}^{M}\left(\hat{\theta}_{m}-\theta_{m}\right)^{2}/M$, and (iii) coverage rate of the 95% credible intervals, $C_{95}$. Results are reported and discussed in Section 3.2. Each model was simulated $M=100$ times and both models were fit to each simulated data set.

For each model fit, two parallel Markov chains were run with dispersed initial values to diagnosis convergence. For each single chain, we ran $2\>000$ iterations and discarded the first $1\>000$ iterations as a warm-up phase, yielding a total of $2\>000$ samples for posterior analysis. The remaining samples were used to calculate posterior summaries of the parameters of interest and model comparison measurements. The potential scale reduction statistic $\hat{R}$ [26], calculated by Stan was used to verify posterior convergence.

3.2 Simulation study results

The estimated potential scale reduction factors were below 1.1 for all parameters, indicating successful convergence. Table 1 summarizes the measures of total bias, MSPE, $C_{95}$, and the model comparison criteria for each scenario. When the data are generated from M2, the true model is most often preferred in terms of either LOOIC (96%) or WAIC (81%). The quartiles of the comparison criteria differences between M2 and M1 are significantly negative, indicating strong evidence that M2 outperforms M1. When the data are generated from M1, the models are similar, and the comparison criteria are performing as expected and guiding us to the true model on average (LOOIC: 59% and WAIC: 65%). Models with criteria differences smaller than two might be considered equivalent. This is expected since M2 accommodates the association between outcomes and includes M1 as a special case. For both scenarios, most parameters of the true model are estimated with better $C_{95}$, lower bias and MSPE, indicating good performance in terms of prediction error. Figure 1 plots the true versus estimated individual intercepts and slopes for the first outcome, and the true versus estimated time shifts of one simulated data set. Both plots demonstrate agreement between the true and estimated values.

3.3 Application to Alzheimer’s Disease Neuroimaging Initiative

The LTJMM was fit to data from the Alzheimer’s Disease Neuroimaging Initiative (ADNI), which has followed volunteers diagnosed as cognitively healthy or with varying degrees of cognitive impairment since 2005 [7]. The ADNI battery includes serial neuroimaging, cerebrospinal fluid (CSF), and other biomarkers; and clinical and neuropsychological assessments. Participants returned for repeated assessments at six months, one year, and every year thereafter. Seven outcome measures were included in the model: CSF tau and amyloid beta (A$\beta$) 1-42; PET imaging of amyloid deposition and glucose metabolism in the brain; volumetric magnetic resonance imaging (vMRI) of the hippocampus; the 13 item Alzheimer’s Disease Assessment Scale (ADAS13); and the Functional Activities Questionnaire (FAQ). Fixed effect covariates for each outcome included age, carriage of the APOE$\varepsilon$4 allele, sex, and education. The model was fit with and without assuming a multivariate Gaussian distribution for the random effects. Three parallel Markov chains were run for $8\>000$ iterations and the first $4\>000$ warm-up iterations were discarded. Every fourth value of the remaining part of each chain was stored to reduce correlation, yielding a total of $3\>000$ samples for posterior analysis. The posterior mean and the 95% credible intervals were calculated using the obtained samples for each parameter.

One goal of the analysis is to compare the long-term trends of the outcomes on a comparable scale and make conclusions about the temporal ordering of their emergence. With this goal in mind, the outcome measures were transformed first to quantiles, and the quantiles were then transformed by the inverse Gaussian quantile function. All the transformed outcomes were then modeled as Gaussian with identity link. ADNI subjects were diagnosed at their first visit with normal cognition (24%), subjective memory concern (6%), early mild cognitive impairment (17%), late mild cognitive impairment (32%), and mild to moderate dementia (19%). A weighted quantile transformation [27] was used so that quantiles approximate the quantiles from a sample with equal numbers of each diagnosis. Categorical diagnosis is based, in part, on the subjective interpretation of the clinical presentation (excluding CSF and imaging data) by the physician, and hence was not included as a covariate in LTJMMs.

We fit the LTJMMs to examine gains by incorporating correlated random effects. Furthermore, we compare our proposed LTJMMs with existing models: the commonly used linear mixed effects model (MM) with and without fixed effects for categorical baseline diagnosis and its interaction with time. Table 3 summarizes the model comparison results. We found that $\textrm{LTJMM}_{2}$ was preferred over $\textrm{LTJMM}_{1}$ and MM without diagnosis. However, MM with diagnosis outperformed our proposed $\textrm{LTJMM}_{2}$ in terms of WAIC, LOOIC and DIC.

Figure 2 shows the subject-level observations (top) and predictions (bottom) according to age. It is clear from the observations, that age explains only a small proportion of the variance in these outcomes. The bottom panel shows that the predictions provide a reasonable smooth of the observations and that latent time provides a reasonable ordering of individuals according to disease severity. The posterior mean (95% credible interval) for the latent time parameter was 14.2 (12.7 to 16.1) years (Table 2). Figure 3 shows a density plot for the posterior mean of the subject-specific latent time by diagnosis at first ADNI visit. The latent time estimates are temporally sorting individuals in a manner that is consistent with physician diagnosis. Latent time parameters provide a continuous alternative to diagnosis which is objectively derived from a comprehensive model of longitudinal measures of disease.

Figure 4 shows the posterior mean of correlation parameters for random intercepts and slopes. Not surprisingly, we see strong correlations between change in cognitive tests (ADAS13) and the function (FAQ). However we also, see strong correlations between measures of symptomatic change (cognition and function) and biomarkers (hippocampal atrophy and glucose metabolism [FDG PET]). The two amyloid measures, CSF A$\beta$ and amyloid PET, show a moderate positive correlation for random intercepts and weakly negative correlation of change.

Figure 5 displays the population-level predicted trajectories. The depicted curves are for female APOEe4 carriers with the ADNI mean education. Recall that age and latent time contribute to the model independently. For the sake of these predictions, age is calibrated so that the ADAS13 trajectory attains the ADNI mean ADAS13 at the ADNI mean age at the first visit. The bottom panel shows the same trajectories for progressive Alzheimer’s (red triangles) with contrasting trajectories for healthy aging (blue dots). To obtain estimates for healthy aging, the effect of latent time is forced to be zero to isolate the effect of age.

4 Discussion

We explored sampling the random effects from both univariate and multivariate Gaussian distributions, and found the multivariate is always preferred when it is the true model, and performs similarly well when the univariate is the true model. The model with multivariate random effects has the advantage of providing additional parameters for the correlation among outcomes (random intercepts) and among change in outcomes (random slopes). However, fitting the model with univariate random effects was often computationally faster, so it may have advantages in practice. We plan to explore extension of the model which include monotone smooth functions of latent time, and other random effects options. Furthermore, the variability of latent time shifts was assumed to be homogeneous (same for all individuals) in this work. We are currently working on an extension of the model to allow heterogeneous variances by explicitly modeling subject-specific latent time variance in terms of observed covariates.

Interestingly, the model selection criteria preferred LTJMM over MM when the MM did not include fixed effects for categorical baseline diagnosis. However, when the MM had the advantage of the additional covariates for baseline diagnosis, the MM was preferred over LTJMM. This analysis suggests LTJMM is a good alternative to MM when baseline diagnosis is not known or when novel insights about long-term progression across all the diagnostic categories are desired. These long-term insights are unattainable with MM. LTJMM also provides subject-specific latent disease-time estimates that could serve as a continuous objective alternative to categorical diagnosis.

The LTJMM provides a parsimonious extension of the existing family of joint generalized linear mixed effects models with a subject-specific latent time parameter that accommodates the temporal heterogeneity common to multicohort longitudinal data in late-life neurodegenerative diseases. The model provides key insights and inference regarding the evolution of disease markers over a period of time that is longer than the period of observation. The framework allows consideration of the independent effects of healthy aging and disease progression.

These novel features can be leveraged to improve subject-level prediction and better understand long-term disease dynamics. In particular we plan to leverage the model to help improve clinical trial design. For example, given an estimate of the short-term effect of a treatment on a biomarker, the model can be interrogated to estimate the downstream effects on cognition and function. The model can also be used to help identify populations expected to experience the maximum benefit from a given intervention. The LTJMM provides a flexible framework to begin to explore these and many other applications.

5 Supplementary material

The Stan code model specifications and related code and functions for simulation, estimation and prediction are available as an R package from https://bitbucket.com/mdonohue/ltjmm. Interactive convergence plots of the model fit to Alzheimer’s data is available from https://shiny.atrihub.org/public/adni_ltjmm/.

Acknowledgment

We are grateful to the ADNI study volunteers and their families. This work was supported by Biomarkers Across Neurodegenerative Disease (BAND-14-338179) grant from the Alzheimer’s Association, Michael J. Fox Foundation, and Weston Brain Institute; and National Institute on Aging grant R01-AG049750. Data collection and sharing for this project was funded by the ADNI (National Institutes of Health Grant U01 AG024904) and DOD ADNI (Department of Defense award number W81XWH-12-2-0012). ADNI is funded by the National Institute on Aging, the National Institute of Biomedical Imaging and Bioengineering, and through generous contributions from the following: AbbVie, Alzheimer’s Association; Alzheimer’s Drug Discovery Foundation; Araclon Biotech; BioClinica, Inc.; Biogen; Bristol-Myers Squibb Company; CereSpir, Inc.; Eisai Inc.; Elan Pharmaceuticals, Inc.; Eli Lilly and Company; EuroImmun; F. Hoffmann-La Roche Ltd and its affiliated company Genentech, Inc.; Fujirebio; GE Healthcare; IXICO Ltd.; Janssen Alzheimer Immunotherapy Research & Development, LLC.; Johnson & Johnson Pharmaceutical Research & Development LLC.; Lumosity; Lundbeck; Merck & Co., Inc.; Meso Scale Diagnostics, LLC.; NeuroRx Research; Neurotrack Technologies; Novartis Pharmaceuticals Corporation; Pfizer Inc.; Piramal Imaging; Servier; Takeda Pharmaceutical Company; and Transition Therapeutics. The Canadian Institutes of Health Research is providing funds to support ADNI clinical sites in Canada. Private sector contributions are facilitated by the Foundation for the National Institutes of Health (www.fnih.org). The grantee organization is the Northern California Institute for Research and Education, and the study is coordinated by the Alzheimer’s Therapeutic Research Institute at the University of Southern California. ADNI data are disseminated by the Laboratory for Neuro Imaging at the University of Southern California.

Appendix A Proof of the Identifiability

Suppose that $p$ response variables are measured for $n$ individuals at different follow-up times. The measured outcomes for individual $i$ at time $j$ are denoted as $\textbf{y}_{ij}=\left(y_{ij1},\cdots,y_{ijp}\right)^{\prime}$. The latent time joint mixed effects model is given by

[TABLE]

where $i=1,\ldots,n$, $k=1,\ldots,p$ and $j=1,\ldots,q_{ik}$.

The model (2) can be re-expressed as the form of linear mixed effects model

[TABLE]

where $\tilde{\alpha}_{0ik}=\gamma_{k}\delta_{i}+\alpha_{0ik}$. Identifiability, up to $\tilde{\alpha}_{0ik}$ follows from the identifiability of standard linear mixed effects models (see Wang (2013)222Wang, W. (2013). Identifiability of linear mixed effects models. Electronic Journal of Statistics, 7, 244-263. for a discussion) provided $(\textbf{x}^{\prime},\textbf{t})$ has full column rank, where t is the vector formed by $t_{ijk}$; and t is linearly independent from the vector of 1’s.

By placing the sum-to-zero constraints on $\alpha_{0ik}$, i.e., $\sum_{k=1}^{p}\alpha_{0ik}=0$, it remains to show the system below:

[TABLE]

has a unique solution, where $\gamma_{k}$ and $\tilde{\alpha}_{0ik}$ are known, since they are identifiable in (3); $\delta_{i}$ and $\alpha_{0ik}$ are unknown. There are $n\times p$ unknown subject and outcome specific parameters $\alpha_{0ik}$, and $n$ subject-specific latent time parameters $\delta_{i}$, for a total of $np+n$ unknown parameters. There are $np+n$ constraints (or equations). The system (4 & 5) with the same number of equations and unknown parameters has a unique solution. Therefore, $\delta_{i}$ and $\alpha_{0ik}$ are uniquely determined. Thus the model is identifiable.