A joint model for multiple dynamic processes and clinical endpoints: application to Alzheimer's disease
Cécile Proust-Lima Viviane Philipps Jean-François Dartigues










A joint model for multiple dynamic processes and clinical endpoints: application to Alzheimer's disease
    Cécile Proust-Lima^∗, Viviane Philipps, Jean-François Dartigues


INSERM, UMR1219, Univ. Bordeaux, ISPED, Bordeaux, France



cecile.proust-lima@inserm.fr
    December 30, 2023
==============================================================================================================================================================









Abstract: As other neurodegenerative diseases, Alzheimer's disease, the most frequent dementia in the elderly, is characterized by multiple progressive impairments in the brain structure and in clinical functions such as cognitive functioning and functional disability. Until recently, these components were mostly studied independently since no joint model for multivariate longitudinal data and time to event was available in the statistical community. Yet, these components are fundamentally inter-related in the degradation process towards dementia and should be analyzed together.
We thus propose a joint model to simultaneously describe the dynamics of multiple
correlated components. Each component, defined as a latent process, is measured by one or several continuous markers (not necessarily Gaussian). Rather than considering
the associated time to diagnosis as in standard joint models, we assume
diagnosis corresponds to the passing above a
covariate-specific threshold (to be estimated) of a pathological process which is modelled as a combination of the component-specific latent processes.
This definition captures the clinical complexity of diagnoses such as dementia diagnosis but also benefits from simplifications for the computation of Maximum Likelihood Estimates. We show that the model and estimation procedure can also handle competing clinical endpoints. The estimation procedure, implemented in a R package, is validated by simulations and the method is illustrated on a large French population-based cohort of cerebral aging in which we focused on the dynamics of three clinical manifestations and the associated risk of dementia and death before dementia.

Keywords: aging; competing risks; dynamic model; joint model; latent process
model; multivariate longitudinal data














A joint model for multiple dynamic processes and clinical endpoints: application to Alzheimer's disease
    Cécile Proust-Lima^∗, Viviane Philipps, Jean-François Dartigues


INSERM, UMR1219, Univ. Bordeaux, ISPED, Bordeaux, France



cecile.proust-lima@inserm.fr
    December 30, 2023
==============================================================================================================================================================










§ INTRODUCTION
 

Dementia is a syndrome which mostly affects individuals 60 years and older, the most frequent type of dementia (approximately 60%) being Alzheimer's disease. Dementia is characterized by a very long degradation process before clinical diagnosis which lasts a few decades (). As other neuro-degenerative diseases, the degradation process in dementia has multiple anatomo-clinical components () including 

an accumulation of biomarkers in the brain (β-amyloid, τ protein), an alteration of the brain structure (atrophy of some regions such as the hippocampus), impaired clinical manifestations with the decline of several cognitive functions (e.g., executive functions, processing speed, episodic memory), an increase of functional limitations in the daily life (e.g. shopping, toileting, transferring from bed to chair), and possibly increased depressive symptoms among other behavioral alterations. 
Although inter-related, these components were mostly studied independently with a large focus on cognitive decline for which repeated measures have been available for long in cerebral aging cohorts (e.g., ). 

The limitation to a single component (also called domain thereafter) was partly explained by the statistical complexity of analyzing multiple longitudinal components in link with a clinical endpoint. It requires the joint modelling of multivariate longitudinal markers and a survival process. A series of joint models for multiple longitudinal and survival data were proposed (see review in  and examples in ) but their application was limited to very simple cases. 
Indeed numerical complexity is considerable in joint models involving multiple longitudinal components due notably to an integral over the random effects in the likelihood which can not be solved analytically. The integral computation becomes cumbersome with more than a couple of random effects () and not necessarily accurately solved (). This entails that most applications focused on only two longitudinal markers (e.g. ) and/or random intercepts only () which is not sensible in complex diseases such as dementia. Contributions also mostly relied on a Bayesian estimation to circumvent the numerical issues () or a two-stage estimation ().

An additional problem in dementia as in many other diseases (especially in psychiatry and neurology) is that each domain is not directly measured; it is approached by multivariate markers. 

For example, cognitive level is usually measured by a battery of psychometric tests to apprehend all the coexisting cognitive functions, and brain structure comprises various regional volumes including the hippocampus volume but not limited to (). Joint models for multivariate longitudinal markers measuring the same underlying process were proposed but were systematically limited to a univariate latent process () such as cognition in dementia ().
Finally, for psychometric and behavioral data at least, the usual assumption of normality does not hold and normalizing transformations have to be incorporated in the model to be able to rely on standard linear mixed regressions ().

In this context, our work aimed at developing a novel joint model for multivariate longitudinal markers measuring several latent processes and one or several clinical events. Instead of considering a standard proportional hazard model for the time to event, as mainly done in joint modelling framework, we considered 

a degradation process model for each clinical endpoint: the clinical endpoint is defined as a binary variable repeatedly collected at visits which becomes positive when its underlying continuous degradation process has passed above an unknown threshold. Using this definition, clinically relevant for dementia diagnosis and many other clinical endpoints, the joint model has an exact likelihood. As such, it can be applied to cases with more than just a couple of longitudinal markers and/or random effects. The replacement of the classical survival model by another model to provide exact likelihood had already been proposed in joint models with a univariate longitudinal Gaussian marker (). In their work, the authors had opted for a sequential probit model defined in discrete time with which our approach has some similarities as explained in discussion. 


Section <ref> describes the joint model for multiple latent domains and one clinical endpoint. Section <ref> details the likelihood computation and Section <ref> extends the model to other types of endpoints and competing endpoints. Section <ref> validates by simulations the estimation procedure which is implemented in a R package. Section <ref> details an application to three clinical manifestations in dementia, cognitive functioning, functional dependency and depressive symptomatology, in link with dementia diagnosis and in the presence of competing risk of dementia-free death. Finally, we conclude in Section <ref>. 




§ THE JOINT MODEL FOR MULTIPLE LATENT DOMAINS AND A CLINICAL ENDPOINT



The joint model for multiple latent domains and a clinical endpoint is described in Figure <ref> and formalized in Subsections <ref>, <ref> and <ref>.









 §.§ Latent domains measured by multivariate longitudinal markers


In a population of N individuals, let consider D latent domains (e.g., cognition, brain structure), each one defined for individual i (i=1,...,N) as a latent process in continuous time (Λ^d_i(t)) with t ∈ℝ and d=1,...,D. Each latent domain d is measured through a battery of K_d markers repeatedly collected over time (e.g., several cognitive tests, volumes of different brain regions). Let define Y^d_kij the measure of marker k (k=1,...K_d) and individual i for domain d collected at time t_dkij with j=1,...n_dki. 

To handle markers that are not necessarily Gaussian, we rely on previous works () and assume that each marker, normalized by a parameterized link function, is a noisy measure of the underlying latent domain:


    H^d_k(Y^d_kij ; η^d_k) = Λ^d_i(t_dkij) + ϵ^d_kij


where ϵ^d_kij are centered independent Gaussian measurement errors with variance σ^d_k^2 and H^d_k(.; η^d_k) is the link function which transforms the marker into a Gaussian framework. This link function depends on parameters η^d_k that are estimated on the data along with the other parameters. Any family of monotonically increasing parameterized transformations can be chosen for H^d_k, including the family of linear transformations which reduces to the standard Gaussian case. When departures from normality are suspected, we recommend the flexible and parsimonious family of linear combinations of quadratic I-splines (IS_l) so that H^d_k(x; η^d_k) = η^d_k0 + ∑_l=1^n_dk+1η^d_kl^2 IS_l(x) with n_dk the number of knots kept small for parsimony (I-splines are integrated M-splines which ensure the monotonicity of the link functions, see ). Previous simulation studies demonstrated that using such nonlinear parameterized link functions yielded correct inference in mixed models when marker distribution deviated from normality (). Note that although not detailed here for sake of brevity, marker-specific effects of covariates and/or random effects can be easily added to this equation of observation (). 



 §.§ Correlated trajectories of latent domains


Each latent domain trajectory is described via a linear mixed model ():


    Λ^d_i(t) = X_id(t)^⊤β^d + Z_id(t)^⊤ b^d_i + w^d_i(t)    , ∀ t ∈ℝ


where X_id(t) is a vector of covariates associated with the vector of fixed effects β^d at the population level and Z_id(t) is a vector of covariates associated with the q_d-vector of random effects b^d_i at the individual level with b^d_i∼𝒩(0,B_dd). A zero-mean Gaussian process w^d_i(t) can be added to make the trajectory more flexible at the individual level; for example a Brownian motion with covariance structure cov(w^d_i(t),w^d_i(u)) = σ_w^d^2min(t,u) is often relevant in aging studies (). We chose to keep the zero-mean Gaussian processes independent across domains. The correlation between the latent domains is only captured by correlations between the domain-specific random effects so that each individual is characterized by an overall vector of random effects: 



    b_i =  ( [ b_i^1;     ⋮; b_i^d;     ⋮; b_i^D;       ] ) 
    ∼𝒩 (
     ( [ 0; ⋮; 0; ⋮; 0;   ] ) ,
    B= ( [   B_11    ...   B_1d    ...   B_1D;      ⋮      ⋱      ⋮      ⋱      ⋮; B_1d^⊤    ...   B_dd    ...   B_dD;      ⋮      ⋱      ⋮      ⋱      ⋮; B_1D^⊤    ... B_dD^⊤    ...   B_DD;        ] )  )



where B_de is the covariance matrix between random effects of latent domains d and e ((d,e) ∈{1,D}^2). The unstructured B matrix was parameterized through the q=∑_d=1^D q_d variances of random effects (parameter σ. to estimate (σ. ∈ℝ) for a variance of σ.^2) and the q×(q-1)2 correlations between pairs of random effects (parameter ρ. to estimate (ρ. ∈ℝ) for a correlation of exp(ρ.)-11+exp(ρ.). 

As in any latent variable model, the dimension (location and scale) of the latent processes need to be defined to reach identifiability. This is achieved by standardizing each latent process d with a zero mean intercept in β^d and a unit variance of the random intercept in b^d_i. 



 §.§ Degradation process toward a clinical endpoint


We assume there exists a degradation process in continuous time for each individual i denoted (Δ^diag_i(t)) with t ∈ℝ. This process is measured by repeated observations of the clinical status of interest, dementia diagnosis in our case, Y^diag_ij at time t^diag_ij with j the occasion (j=1,...,n^diag_i) so that an individual has the positive clinical status if the degradation process with an additional noise has reached a threshold (to estimate) at the visit time:

    Y^diag_ij = 1    ⇔   Δ^diag_i(t^diag_ij) + ϵ^diag_ij ≥ ζ^diag_0  +  X^diag_ij^⊤ζ^diag


where ϵ^diag_ij is a zero-mean Gaussian independent random variable (ϵ^diag_ij∼𝒩(0,1)), ζ^diag_0 is the threshold parameter in the reference group and X^diag_ij is a vector of covariates associated with parameters ζ^diag that can modulate the threshold defining the positive clinical status. Note that observations are no more considered after a change to a positive clinical status. The variance of ϵ^diag_ij is constrained to 1 to ensure identifiability as shown in supplementary material, Section 1.1.

The degradation process is defined as a linear combination of the latent domains: 

    Δ^diag_i(t) = γ_1iΛ^1_i (t) + ... + γ_diΛ^d_i (t) + ... + γ_DiΛ^D_i (t)


where γ_di is the intensity of latent domain d contribution for individual i. It can be contrasted according to covariates with γ_di =  X^diag_di^⊤γ_d in which X^diag_di is a vector of covariates (including the intercept) associated with the vector of parameters γ_d. For instance, in dementia, the contributions of latent domains may be modulated by factors linked to cognitive reserve.



§ MAXIMUM LIKELIHOOD ESTIMATION



The entire vector of parameters denoted θ comprises parameters {η_k^d,d=1,...,D,k=1,...,K_d }, {β^d,σ_w^d, σ^d_k, γ_d, d=1,...,D}, ζ^diag_0, ζ^diag and all the parameters σ and ρ constituting the B matrix. It is estimated in the maximum likelihood framework. 



 §.§ Likelihood


Let denote Y^diag_i= (Y^diag_ij)_j=1,...,n^diag_i the repeated observed diagnoses, Y_i=(Y^d_i)_d=1,...,D the repeated and multivariate observed markers (with Y^d_i the vector of all the repeated and multivariate observations of latent domain d) and Λ_i =  ( (Λ^d_i(t^diag_ij))_j = 1,...,n^diag_i )_d=1,...,D the set of latent processes at the observed diagnosis visits. By denoting f(X) the generic probability density function of a random variable X, the likelihood is L(θ) = ∏_i=1^N f(Y^diag_i,Y_i; θ) with 

    f(Y^diag_i,Y_i)    = ∫_ℝ^D n^diag_i f(Y^diag_i,Y_i|Λ_i) f(Λ_i)dΛ_i = ∫_ℝ^D n^diag_i f(Y^diag_i|Λ_i)f(Y_i|Λ_i)f(Λ_i)dΛ_i
       = f(Y_i)  ∫_ℝ^D n^diag_i  f(Y^diag_i|Λ_i)  f(Λ_i|Y_i)dΛ_i

where the vector of parameters θ is omitted for simplicity and 


  * the marginal density of the repeated markers f(Y_i) can be decomposed into the multivariate normal density of the transformed observations of the markers through the link functions ϕ_H(Y)(H(Y_i);μ_H(Y),V_H(Y)) times the Jacobian of the link functions J(H(Y_i))=∏_d=1^D ∏_k=1^K_d∏_j=1^n_dki J(H^d_k(Y^d_kij ; η^d_k)). The multivariate normal density of the transformed observations H(Y) has mean μ_H(Y)=X_i β and variance V_H(Y)= Z_iBZ_i^⊤ + R_i + Σ_i with β= (β^d)_d=1,...D, X_i and Z_i the D-block diagonal matrices with blocks X^d_i and Z^d_i that comprise row vectors X_id(t)^⊤ and Z_id(t)^⊤ for all the outcome-specific occasions t,  R_i the covariance matrix of the Gaussian processes (w^d_i(t)) at the same visits and Σ_i the diagonal variance matrix of measurement errors;


  * Following equations (<ref>) and (<ref>), the conditional distribution of the repeated diagnoses is defined as follows (detailed calculations are provided in supplementary material, Section 1.2):



  * f(Y^diag_i|Λ_i)=Φ^(n^diag_i) (ζ_i - Γ_i Λ_i ; 0, I_n^diag_i ) for subjects not diagnosed with dementia where Φ^(n)(x;m,V) is a n-dimensional Gaussian cumulative distribution function with mean m and variance V computed in x, ζ_i=  ( ζ^diag_0  +  X^diag_ij^⊤ζ^diag )_j = 1,...,n^diag_i, Γ_i =  ([ γ_1iI_n^diag_i            ... γ_diI_n^diag_i            ... γ_DiI_n^diag_i ] ) and I_n is the n× n identity matrix;

  * f(Y^diag_i|Λ_i)=Φ^(n^diag_i-1) (ζ_i^- - Γ_i^- Λ_i ; 0, I_n^diag_i-1 ) - Φ^(n^diag_i) (ζ_i - Γ_i Λ_i ; 0, I_n^diag_i-1 ) for those diagnosed with dementia; ζ_i^- and Γ_i^- are respectively ζ_i and Γ_i without the last row for t^diag_in^diag_i;



  * the conditional density f(Λ_i|Y_i) is a multivariate normal density function ϕ_Λ_i(Λ_i;μ_Λ_i,V_Λ_i) with mean μ_Λ_i= X̃_i β + V_Λ H(Y)V_H(Y)^-1 (H(Y_i) -X_i β) and variance V_Λ_i = Z̃_iBZ̃_i^⊤ + R̃_i -
V_Λ H(Y)V_H(Y)^-1V_Λ H(Y)^⊤ where X̃_i and Z̃_i are the D-block diagonal matrices with blocks X̃_i^d and Z̃_i^d that comprise row vectors X_id(t)^⊤ and Z_id(t)^⊤ at diagnosis visits t^diag_ij (j=1,..., n^diag_i), R̃_i the covariance matrix of the Gaussian processes (w^d_i(t)) at the same visits and R̃̃̃_i the covariance matrix of the Gaussian processes (w^d_i(t)) between diagnosis visits and marker visits so that V_Λ H(Y)= Z̃_i B Z_i^⊤ + R̃̃̃_i. 




The joint density of the observations can thus be rewritten:


    f(Y^diag_i   ,Y_i) =ϕ_H(Y)(H(Y_i);μ_H(Y),V_H(Y))J(H(Y_i))   ×
       ∫_ℝ^Dn_i^diagΦ^(n^diag_i) (ζ_i - Γ_i Λ_i ; 0, I_n^diag_i )
    ϕ_Λ_i(Λ_i;μ_Λ_i,V_Λ_i) dΛ_i   if  Y^diag_ij = 0  ∀ j 
    ∫_ℝ^Dn_i^diag ( Φ^(n^diag_i-1) (ζ_i^- - Γ_i^- Λ_i ; 0, I_n^diag_i-1 ) . -    
    . Φ^(n^diag_i) (ζ_i - Γ_i Λ_i ; 0, I_n^diag_i )  )   ϕ_Λ_i(Λ_i;μ_Λ_i,V_Λ_i)dΛ_i   if  Y^diag_in^diag = 1


Instead of approximating the multivariate integral over Λ_i by numerical integration techniques (as mostly done in joint modelling area with Gaussian quadratures), we use the properties of the skew-normal variables which provide a closed form for the integrals in equation (<ref>) () (see Appendix for the general formula). The joint density becomes:


    f(Y^diag_i,Y_i)    = ϕ_H(Y)(H(Y_i);μ_H(Y),V_H(Y))J(H(Y_i))   ×
       Φ^(n^diag_i) (ζ_i - Γ_i μ_Λ_i ; 0, I_n^diag_i  + Γ_i V_Λ_iΓ_i ^⊤ )     if  Y^diag_ij = 0  ∀ j 
    Φ^(n^diag_i-1) (ζ_i^- - Γ_i^- μ_Λ_i ; 0, I_n^diag_i-1  + Γ_i^- V_Λ_iΓ_i^-^⊤ )  -    
    Φ^(n^diag_i) (ζ_i - Γ_i μ_Λ_i ; 0, I_n^diag_i  + Γ_i V_Λ_iΓ_i ^⊤ )     if  Y^diag_in^diag = 1





 §.§ Likelihood accounting for delayed entry


When necessary, delayed entry can be accounted for by dividing the likelihood by the probability of being at risk of the clinical endpoint at study entry that is Y^diag_i1=0 at time t^diag_i1:

    L^delayed(θ) = L(θ)∏_i=1^N P(Y^diag_i1=0)

where 

    P(Y^diag_i1=0)    = ∫_ℝ^DΦ^(1) (ζ_i^0 - Γ_i^0 Λ_i^0 ; 0, 1  )
    ϕ(Λ_i^0; X_i^0 β,Z_i^0 B Z_i^0^⊤+R_i^0)dΛ_i^0
       = Φ^(1) ( ζ_i^0 - Γ_i^0 X_i^0 β; 0; 1 +  Γ_i^0 (Z_i^0 B Z_i^0^⊤+R_i^0) Γ_i^0^⊤ )



with Λ_i^0 = (Λ_i^d(t^diag_i1))_d=1,...,D the latent processes at entry, ζ_i^0 = ζ^diag_0 + X^diag_i1^⊤ζ^diag the threshold at entry, Γ_i^0=(γ_1i,...,γ_Di) the vector of domain-specific contributions, X_i^0 and Z_i^0 the D-block diagonal matrices with blocks X^d_i(t^diag_i1) and Z^d_i(t^diag_i1), and R_i^0 the variance of the Gaussian processes at entry w^d_i(t^diag_i1) for d=1,...,D.




 §.§ Likelihood Optimization and implementation


The log-likelihood is maximised using a modified Marquardt algorithm, a Newton-like algorithm, with convergence criteria based on parameter and log-likelihood stability and derivatives size  (). The latter is defined at iteration l as ∇(L(θ^(l)))^⊤ℋ^(l)-1∇(L(θ^(l)))/n_θ≤ω with ∇(L(θ^(l))) and ℋ^(l) the gradient and the Hessian matrix of the log-likelihood at iteration l, n_θ the length of θ, and ω the convergence threshold (fixed here at 0.001). This criterion is very stringent and ensures convergence toward a maximum. The program was implemented in C++ with an interface in R and parallel computations to fasten the estimation procedure. The R package called  can be downloaded on github: (<https://github.com/VivianePhilipps/multLPM>). The program fits any model that has the same structure as defined in Section <ref> (and possibly a second competing event as defined in Section <ref>). The multivariate normal cumulative distribution functions were computed by using Genz routines (). 




§ EXTENSION TO MULTIPLE CLINICAL ENDPOINTS AND CONTINUOUS-TIME EVENTS




 §.§ Multiple endpoints at predefined visits


The definition of the degradation process and measurement model for a clinical event observed at specific visits such as diagnoses (Subsection <ref>) extends naturally to the case of two causes of diagnosis (with cause-specific degradation processes). The log-likelihood has the exact same structure except that the dimension of the cumulative distribution functions is augmented to the number of visits for each cause of clinical event, 2× n^diag_i if the two diagnoses are made at the same times. 



 §.§ Event in continuous time


The model definition relies on discrete repeated visits for the clinical event and as such does not directly extend to clinical events that are intrinsically defined in continuous time such as death. With events in continuous time, a preliminary discretization of time is necessary: we partition time into S intervals ℐ_s=[l_s,u_s] with midpoints m_s for s=1,...,S (l_s=u_s-1 for s>1). A subject is followed from interval s_0i, the interval containing the exact entry time in the study, and is followed until interval s_i: if the subject has the event, s_i is such that the time of event is in ℐ_s_i; if the subject is censored, s_i is such that the censoring time is in ℐ_s_i+1. Without loss of generality, we focus on death and define the repeated death status Y^death_ij in each interval j (j=s_0i,...,s_i) with Y^death_ij=0 for all j except j=s_i for those who die (Y^death_is_i=1). 









Using the exact same definition as for clinical diagnoses, we consider a latent degradation process defined as a linear combination of the latent domains:

    Δ^death_i(t) = δ_1iΛ^1_i (t) + ... + δ_diΛ^d_i (t) + ... + δ_DiΛ^D_i (t)


We then define the probability of dying in interval s as the probability that the noisy underlying degradation process is above a specific threshold at the midpoint of interval s:

    P(Y^death_is = 1) = P(Δ^death_i(m_s) + ϵ^death_is ≥ ζ^death_0  +  X^death_is^⊤ζ^death)

As previously, δ_di characterizes the intensity of latent domain d contribution for individual i and it can be contrasted according to covariates with δ_di = X^death_di^⊤δ_d in which X^death_di is a vector of covariates (including the intercept) associated with the vector of parameters δ_d; ϵ^death_is is a zero-mean Gaussian independent random variable (ϵ^death_is∼𝒩(0,1)); ζ^death_0 is the threshold in the reference group and X^death_is is a vector of covariates associated with parameters ζ^death that can modulate the threshold defining the death status. In particular, X^death_is may include time. 

When this discretized event replaces clinical diagnosis, the likelihood has the same structure except that the dimension of the cumulative distribution function becomes s_i-s_0i+1. In case of delayed entry, the likelihood is divided by the probability to have survived until l_s_0i.





 §.§ Endpoint at predefined visits in competition with an event in continuous time


Without loss of generality, we take the example of a clinical diagnosis for which death before diagnosis constitutes a competing event. The methodology handles this by combining the two degradation processes of sections <ref> and <ref>, and by slightly changing the observations for death. To focus on death before diagnosis, we only consider death in the x years following a negative diagnosis (x=3 years is realistic in the case of dementia). Otherwise, time to death is censored at the last visit t_in_i^diag^diag. 

With these changes, the joint log-likelihood has again the exact same structure as in section <ref> except that the dimension of the cumulative distribution functions is augmented to the number of visits for the clinical event and the number of intervals for the discretized event n^diag_i + s_i-s_0i+1. In case of delayed entry, the likelihood is divided by the probability to have survived until s_0i-1 and still be at risk of clinical event at entry.




§ NUMERICAL EVALUATION OF THE METHODOLOGY




We evaluated the methodolgy with two series of simulations. We first validated the estimation procedure under a correctly specified model by investigating different combinations of key simulation parameters (i.e., number of subjects, proportion of events, nature of the event and uncertainty in the markers measurement). We then evaluated the behavior of the method under different types misspecification (i.e., distribution of errors in the degradation model, correlation between the latent domains, second clinical endpoint or additional domain not modelled). For all scenarios, design and parameter values were inspired by the application data, and 200 replications were done. All the scenarios are summarized in supplementary Table S1.



 §.§ Simulation study I: Validation of the estimation procedure




We considered two longitudinal domains, namely cognition and functional dependency, measured by two psychometric tests and one scale, respectively. We generated linear trajectories of the two domains according to age and adjusted for binary educational level (0.5-probability Bernoulli) and their interaction for domain 1. Correlated random intercept and slope with age took into account the correlation within and between domains at the individual level. Entry in the cohort was generated from a normal distribution (mean 75, standard deviation 3). Observed visits were generated with a uniform distribution in [-1,1] years around the theoretical visits every 2.5 years up to 20 years, and a 15% dropout was considered at each visit. We primarily considered that the two domains contributed to the clinical event model and that a clinical diagnosis was made at each visit. Longitudinal data were censored after the event. We considered Gaussian longitudinal markers with different ranges inspired by psychometric tests IST and DSST for cognition and IADL sum-score for function (see application Section <ref>). More details on the generating procedure are given in supplementary section 2.1.

In the main scenario (I.1.a), generating parameters roughly corresponded to those obtained on the application data and samples included 500 subjects. Over 200 replicates, this lead on average to 4.3 visits per subject and 150 diagnoses (30%). As reported in Table <ref>, parameters were well estimated without bias and no departure from the expected 95% coverage rate of the 95% confidence interval. In additional simulations, we also considered: (scenario I.1.b) smaller samples of 200 subjects which lead in mean to 59 (29.5%) diagnoses; (scenarios I.2.a and I.2.b) smaller proportion of diagnoses by changing the values of each domain contribution. This lead on average to 56 (11%) and 22 (11%) diagnoses in samples of 500 and 200 subjects, respectively; (scenario I.3) larger measurement errors for the markers; (scenario I.4) continuous event discretized into 10 intervals instead of a clinical diagnosis made at each visit. In all these scenarios, inference remained correct with negligible bias and satisfying coverate rates (see supplementary Tables S2 to S6). With 200 subjects, variances were systematically higher though. Overall, these simulations validated the estimation procedure. 






 §.§ Second simulation study: Behavior under misspecification


Based on the same overall design of simulations, we investigated 4 types of misspecification. Results are summarized in supplementary material, section 2. When generating logistic errors in the degradation process toward diagnosis (in equation (<ref>)) instead of Gaussian errors (supplementary Table S7) or when neglecting the intra-individual correlation between latent domains (supplementary Table S8), inference quality was not much affected. When a second event was not modelled (supplementary Table S9), results depended on the level of association between domains and neglected competing event. Inference was not affected with very small contributions of the domains but as expected, parameters of the degradation process model became slightly biased with increased contributions (as it lead to informative censoring). When a third correlated latent domain was not modelled, the inference at the latent domain (and marker) level remained correct. The estimates of the degradation process differed depending on the intensity of association with the neglected domain but this was expected as the interpretation of the contributions (adjusted or not for other domains) differ (supplementary Table S10).





§ APPLICATION TO CLINICAL MANIFESTATIONS IN DEMENTIA



Changes in various clinical measures such as cognitive tests, dependency scales or depressive symptomatology have been separately observed in prodromal dementia () suggesting a possible concomitant role in the dementia process with a modulation of the intensity by educational level which illustrates a potential compensatory mechanism. Our objective was to precisely investigate the role of cognition, dependency and depression in the degradation process toward dementia by jointly analyzing their trajectories and their determinants in link with dementia diagnosis and accounting for the competing death. 



 §.§ The PAQUID data

We relied on the data from the population-based PAQUID cohort which included 3777 individuals aged 65 years and older and living at home in southwestern France in 1989-1990. Individuals were then followed for up to 25 years with repeated neuropsychological evaluations and clinical diagnoses of dementia every 2 or 3 years () and death continuously recorded. We focused on a subsample of 646 individuals who were tested for ApoE4, the main genetic factor associated with aging. We excluded 22 individuals diagnosed with dementia at baseline and 31 individuals who did not have at least one observation of each clinical measure. The final sample consisted of 593 individuals among which 180 developed a dementia and 283 died in the three years following a negative dementia diagnosis. The sample comprised 332 (56%) women, 438 (74%) individuals with higher education level (EL+; individuals who graduated from primary school) and 130 (22%) carriers of at least one copy of the APOE4 allele. The mean age at entry was 73.6 (SD=6.1) years.

The three clinical manifestations were: 


  * Cognitive impairment; it was assessed by four psychometric tests (inverted so that higher levels indicated higher impairment). The Mini-Mental State Examination (MMSE) provides an index of global cognitive performance, the Benton Visual Retention Test (BVRT) assesses visual memory, the Isaacs Set Test (IST) measures verbal semantic memory and speed, and the Digit Symbol Substitution Test (DSST) provides a global measure of executive functioning and processing speed.

  * Functional dependency; it was assessed by the French version of the Instrumental Activities of Daily Living (IADL). We summed the grades of dependency for four activities, telephone use, transportation, medication and domestic finances.

  * Depressive symptomatology; it was assessed by the sum-score of the Center for Epidemiologic Studies Depression Scale (CES-D).


Individuals had between 1 and 12 repeated measures of each marker with a median of 6 (Interquartile range IQR=4-8) for MMSE, IADL and dementia diagnosis, 5 (IQR=4-8) for CES-D, 5 (IQR=3-8) for BVRT, 4 (IQR=3-7) for IST and 4 (IQR=2-6) for DSST. 



 §.§ Model specification


The structure of the model is summarized in Figure <ref> and specifics are given below. 




Quadratic trajectories according to age were assumed for each domain with an adjustment for age at entry, EL+, ApoE4 and gender (and their interactions with age and age squared), and three individual random effects (on intercept, slope and slope squared) correlated within and between domains. An additional time dependent variable indicating the baseline evaluation was included due to evidence of a primo passation negative effect. The selection of covariates and interactions was determined in separate analyses by domain with a 20% significance level. 

Marker-specific link functions used in equation (<ref>) to normalize the markers were quadratic I-splines functions with 3 internal knots placed at the quartiles of the marker distribution over follow-ups. The relevance of the splines transformations was checked in domain-specific analyses by comparing the Akaike criterion (AIC) with the AIC of the model assuming a linear link. 

Degradation process toward dementia was modelled according to the three domains with an adjustment for age at entry, gender, ApoE4 carriers and EL+ and a potential change in the threshold of dementia after the 10-year follow-up visit (due to a new drug on the market that implied earlier diagnoses). 

Degradation process toward dementia-free death was also modelled according to the three domains. We discretized death into 8 intervals with boundaries at 65, 70, 74, 77, 80, 83, 86, 90 and 104 years chosen according to the distribution of death times. The threshold defining the probability of dying was modelled according to age using natural cubic B-splines with knots at 70, 83, 90 and 100 years.

As for any complex model, we recommend to estimate the model progressively by first fitting domain-specific mixed models separately then jointly and finally with dementia and death models. This provides at each step reasonable initial values and reduces computational time. 



 §.§ Results


Fixed effects of the multivariate mixed model for cognition, functional dependency and depressive symptoms are given in Table <ref> and predicted trajectories by covariate profile are displayed in Figure S1 of Supplementary Materials. In summary, each domain had a quadratic trajectory with age characterized by an acceleration in older ages. Individuals included at an older age were systematically more impaired than those included at an younger age (whatever the current age) highlighting a cohort effect. The first passing indicator confirmed that evaluations at baseline underestimated cognitive level and depressive symptoms while overestimating functional dependency. More educated individuals had better cognitive and functional levels. ApoE4 carriers had a faster increase of cognitive and functional impairments. Finally, men had a lower level of depressive symptoms at 65 years but the difference with women reduced when age increased. Men also had a slower increase of functional dependency. 









Estimates of the submodel for the degradation process toward dementia are given in Table <ref>. 

Adjusted for covariates, cognitive and functional domains significantly contributed to the degradation process toward dementia with increased impairments associated to higher degradation of the dementia process and a higher weight of function among ApoE4 carriers. Interestingly, depressive symptoms did not contribute at all to the degradation process toward dementia among individuals with a high educational level (θ=-0.447+0.452=0.005, p=0.957), but it highly contributed among those with a low educational level (θ=-0.447, p=002); in this group, higher depressive symptoms induced a lower level of the degradation process for the same level of cognition and function. This result may be explained by the cognitive reserve linked with education. Among individuals who did not reach a sufficient educational level, having higher depressive symptoms alters the cognitive evaluation while individuals who reached a sufficient educational level are able to compensate and properly pass the cognitive evaluation. Finally, the threshold at which dementia is diagnosed differed according to covariates, with dementia diagnosed at a lower level of degradation for those with a higher educational level, a younger age or a diagnosis made after the 10year visit. No significant difference was found according to gender or ApoE4 status. 

To illustrate these differences in the structure of the degradation process, Figure <ref> displays the mean predicted trajectories of the degradation process according to education, age at entry and ApoE4 status. Note that here, the degradation process was recentered by absorbing the covariate-specific threshold so that dementia is diagnosed when the process is above 0. ApoE4 status constitutes the main modulating factor of the trajectory toward dementia with ApoE4 carriers diagnosed with a dementia 5 to 6 years before ApoE4 non carriers. In contrast, despite a different structure of the degradation process according to education, very limited differences were found which highlights that the differential contribution of depressive symptoms mainly served as a compensating factor in the degradation process definition. 


We investigated the added value of simultaneously considering multiple domains in association with dementia and death rather than considering each domain separately. As longitudinal information differed, we focused on information criteria for dementia and death conditional to longitudinal information with conditional Akaike and Bayesian Information Criteria (AIC, BIC) () and Universal Approximate Cross Validation criterion (UACV) () (See calculations in supplementary Section 3.1). All criteria concluded to the large superiority of the multivariate model over univariate models with a difference of at least 0.1 in UACV, 130 points in AIC and 109 points in BIC (Table  <ref>). From an epidemiologic perspective, conclusions also differed. In univariate models, each domain contributed positively to the degradation process toward dementia even depressive symptomatology although, adjusted for cognition and function, depression was no more positively associated with the degradation process but contributed only among low educated people and negatively as a compensatory factor.







 §.§ Supplementary results and analyses



The estimated splines transformations linking each marker to its underlying domain (Figure S2 in Supplementary Materials) exhibited clear nonlinear relationships except for IST and BVRT cognitive tests. This illustrates the varying sensitivity to change of scales observed in previous works () and the relevance of taking such nonlinear relationships into account. The correlation matrix between the 9 individual random effects also exhibited very high correlations between the level at 65 years, age and age^2 within and between cognitive, functional and depression domains with correlations up to 0.90 between slopes of cognitive and functional domains (Figure S3 in Supplementary Materials). In a secondary analysis, we considered the domains as independent. Although epidemiological results regarding the degradation processes toward death and dementia did not change, the fit of dementia and death data was substantially impacted (conditional AIC=2021.0, conditional BIC=2099.9, UACV=1.710). 

The goodness-of-fit of the model was check by verifying that the individual predictions were closed enough to the observations in the multivariate mixed submodel (Figure S4 in Supplementary Materials). We also verified that the model for dementia-free death considering thresholds approximated by splines was flexible enough by comparing it with a model considering interval-specific probability of death (Figure S5 in Supplementary Materials). 







§ DISCUSSION


We proposed a novel joint model for multiple longitudinal dimensions and clinical endpoints. Initially motived by the study of one repeated binary clinical endpoint measured at predefined visits, the model also applies to continuous endpoints such as death (provided they can be discretized) and to competing risk setting as shown in the application. 

The complexity in joint models, as induced for instance by multiple longitudinal markers, usually orientates the model development toward Bayesian approaches which may better tackle numerical problems. An originality of this work is that the estimation in the frequentist framework was made possible thanks to properties of skew-normal distributions which avoided the cumbersome multiple numerical integration over the random effects usually encountered in joint shared random effect models (). As such, this approach can be applied to contexts where the number of random effects becomes substantial (such as 9 in the application) and/or correlated Gaussian processes are added to relax the model. Being able to jointly model a substantial number of markers becomes indeed a real challenge with the growing availability of dynamic markers in longitudinal health studies.


The properties of skew-normal distributions had already been exploited to simplify inference in joint models in the case of a single repeated biomarker and a single continuous time to event (). The authors had opted for a sequential probit model for the discretized time where we opted for a degradation process model (Section <ref>). Although different in their definition, the two approaches are numerically equivalent in the absence of delayed entry as shown in Section 1.3 of Supplementary Materials. 


In Alzheimer's disease, this model gives for the first time an opportunity to describe the multidimensional pathological process toward dementia. We considered cognitive, functional and depressive symptomatology measures to understand how they contributed to the degradation process toward dementia. Taken independently, each component was associated with the degradation toward dementia. However, when considered jointly, we found an interesting compensatory mechanism of depressive symptomatology among lower educated participants. Depressive symptomatology seems to counterbalance their cognitive evaluation which does not correctly translate their actual cognitive level in the case of depressive symptomatology, probably due to a poorer cognitive reserve. This mechanism impacted only the definition of the dementia degradation process, not its level. We also confirmed and quantified the higher propensity of ApoE4 carriers to develop dementia. Next step will be to also consider neurodegeneration information to capture the anatomic component of dementia ().


The approach however has several limits. Although elegant, the technique is limited to manageable dimensions of repeated clinical endpoints at the individual level. The likelihood calculation relies on algorithms to compute the multivariate normal cumulative distribution function. Their good precision may become questionable beyond very large dimensions which probably limits the methodology to no more than two competing events. Such limitation remains at the individual level, not at the population level which may still include a lot more intervals and/or diagnosis visits. In addition, we believe this is a reasonable concession in many applications where multivariate longitudinal dimensions with a large amount of random effects are to be modelled and could not by using the numerical approaches previously proposed in the literature. Another limitation is that the methodology does not apply to repeated binary, ordinal or count markers. However, by using parameterized link functions, we can handle continuous non Gaussian markers which permits to correctly analyze many asymmetric scales (). 
For the latent process models, we chose to assume that one marker was the manifestation of only one latent process as usual in latent variable methodology and that the possible correlation between the markers was handled at the latent processes level. This was relevant in Alzheimer's disease where domains under study are defined from distinct families of markers. However, an interesting alternative may be to consider independent principal latent processes built from all the markers in the spirit of principal component analyses. Finally, to adapt to the context of Alzheimer's disease, we defined latent processes measured by multiple markers which may not be necessary in other contexts where markers directly constitute the components under study. The methodology and the program made available under R apply in such a case as the standard multivariate linear mixed model constitutes a specific case of the approach. 


To conclude, by handling multivariate repeated markers and clinical endpoint, this joint model opens to many applications in which identifying markers of clinical progression are of importance. The definition of the clinical endpoint differs from most joint model proposals but it is actually clinically relevant in many diseases characterized, as in neurodegenerative diseases, by a body of progressive impairments. 







 §.§ Acknowledgements

This work was supported by project SMALA (ANR-15-CE37-0002) of the French National Research Agency (ANR) and project 2CFaAL of France Alzheimer Association. Computer time was provided by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) at the Université de Bordeaux and the Université de Pau et des Pays de l'Adour.




 §.§ Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request with exception to the application raw data from PAQUID cohort. The software is openly available in github at https://github.com/VivianePhilipps/multLPM.



biom






We used the following equality obtained by <cit.> (equation 59 of the paper) for skew-normal variables to compute an exact formula of the joint model likelihood:  


    Φ^(m) ( λ_0 ; 0, Δ + ΛΛ^T  ) = ∫_ℝ^kΦ^(m) ( λ_0 + Λ z ; 0, Δ ) ϕ^(k) (z) dz

where λ_0 is a m-vector, Δ is a m× m variance covariance matrix and Λ is a m× k matrix. As in the main text, Φ^(n)(x;m,V) denotes the n-dimensional Gaussian cumulative distribution function with mean m and variance V computed in x, and ϕ^(k) is the k-dimensional standard Gaussian density function. 


