Joint models for the longitudinal analysis of measurement scales in the presence of informative dropout

10/06/2021 ∙ by Tiphaine Saulnier, et al. ∙ Université de Bordeaux 0

In health cohort studies, repeated measures of markers are often used to describe the natural history of a disease. Joint models allow to study their evolution by taking into account the possible informative dropout usually due to clinical events. However, joint modeling developments mostly focused on continuous Gaussian markers while, in an increasing number of studies, the actual marker of interest is non-directly measurable; it consitutes a latent quantity evaluated by a set of observed indicators from questionnaires or measurement scales. Classical examples include anxiety, fatigue, cognition. In this work, we explain how joint models can be extended to the framework of a latent quantity measured over time by markers of different nature (e.g. continuous, binary, ordinal). The longitudinal submodel describes the evolution over time of the quantity of interest defined as a latent process in a structural mixed model, and links the latent process to each marker repeated observation through appropriate measurement models. Simultaneously, the risk of multi-cause event is modelled via a proportional cause-specific hazard model that includes a function of the mixed model elements as linear predictor to take into account the association between the latent process and the risk of event. Estimation, carried out in the maximum likelihood framework and implemented in the R-package JLPM, has been validated by simulations. The methodology is illustrated in the French cohort on Multiple-System Atrophy (MSA), a rare and fatal neurodegenerative disease, with the study of dysphagia progression over time truncated by the occurrence of death.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

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


In health cohort studies, repeated measures of markers are often used to describe the natural history of a disease. Joint models allow to study their evolution by taking into account the possible informative dropout usually due to clinical events. However, joint modeling developments mostly focused on continuous Gaussian markers while, in an increasing number of studies, the actual marker of interest is non-directly measurable; it consitutes a latent quantity evaluated by a set of observed indicators from questionnaires or measurement scales. Classical examples include anxiety, fatigue, cognition. In this work, we explain how joint models can be extended to the framework of a latent quantity measured over time by markers of different nature (e.g. continuous, binary, ordinal). The longitudinal submodel describes the evolution over time of the quantity of interest defined as a latent process in a structural mixed model, and links the latent process to each marker repeated observation through appropriate measurement models. Simultaneously, the risk of multi-cause event is modelled via a proportional cause-specific hazard model that includes a function of the mixed model elements as linear predictor to take into account the association between the latent process and the risk of event. Estimation, carried out in the maximum likelihood framework and implemented in the R-package JLPM, has been validated by simulations. The methodology is illustrated in the French cohort on Multiple-System Atrophy (MSA), a rare and fatal neurodegenerative disease, with the study of dysphagia progression over time truncated by the occurrence of death.


joint model, time-to-event data, longitudinal data, latent process, ordinal data


  • [label=•]

  • Joint models are a useful tool to analyze longitudinal data with informative dropout

  • joint Latent Process Models (LPM) allow to describe a latent quantity measured by noisy indicators over time

  • Maximum Likelihood Estimation of the joint LPM is available in the R-package JLPM

  • The case study quantifies the progression of dysphagia during MSA disease and its association with death risk



Graded Response Model


Hazard Ratio


Item Response Theory


Joint Latent Process Model






Magnetic Resonance Imaging


Multiple-System Atrophy (disease)


MSA with predominant Cerebellar impairment


MSA with predominant Parkinsonism


Unified Multiple-System Atrophy Rating Scale

1 Introduction

In health cohort studies, markers measured at multiple times are often used to describe the natural history of a disease, monitor patients or predict the risk of clinical progression. Classical examples include T-cell CD4 counts and viral load for the progression of HIV [1] or PSA for prostate cancer evolution [2]. Due to the intrinsic intra-subject correlation between the repeated measures of markers, their evolution can not be modelled using classical regression models, and mixed models which include individual random effects to account for this serial correlation over time are now adopted worldwide [3].

During the follow-up, clinical events (e.g. diagnosis, recurrence, death) may disrupt the progression of the markers and induce a dropout. When this truncation of the follow-up is predictible by the observed marker data, the missing data mechanism is called missing-at-random (MAR), and inference provided by the mixed model is still valid [4]. However, in many cases, the dropout is likely to depend on the underlying (unobserved) disease mechanism rather than only on the strictly observed data. The missing data mechanism becomes missing-not-at-random (MNAR), and mixed models may not provide correct inference anymore [5, 6]. During the last twenty years, the statistical community has massively embraced the issue of dropout in longitudinal analysis which lead to the development of joint models for the simultaneous analysis of repeated markers and clinical events [7, 8, 1, 9]

. Joint models combine a mixed model describing the progression of the markers and a survival model for the time of occurence of the clinical events, while associating the two models through shared random variables, usually the random effects from the mixed model

[7]. Beyond the study of the marker progression in the presence of MNAR dropout, this method also more generally assesses the association between a marker trajectory and an associated event of interest [10, 11, 2].

Despite many developments in joint models in the recent years [12, 13], most works are dedicated to classical continuous biomarkers stemmed from blood samples, MRI, etc. Yet, in an ever increasing number of health studies, the actual marker of interest is not directly measured. It is a latent construct which is assessed using a set of noisy indicators usually stemmed from questionnaires or measurement scales. Examples include health related quality of life in Cancer research and beyond [14], cognitive functioning and functional dependency in neurodegenerative diseases [15, 16], or many other constructs such as anxiety or depressive symptomatology [17, 18]. The analysis of latent constructs stemmed from measurement scales and questionnaires requires a specific attention. Measurement scales usually translate into multiple categorical and/or continuous items that measure different aspects of the underlying construct of interest. Furthermore, when aggregating the item information into a sumscore, the resulting univariate marker may not have classical Gaussian properties for continuous markers; they are usually bounded with floor/ceiling effects and possible unequal-interval scaling [16, 19].

In this contribution, we show how the joint modeling methodology can be extended to handle repeated data from measurement scales in the presence of an informative clinical event, and provide a software solution with the R package JLPM. We illustrate the methodology though simulations and in Multiple-System Atrophy (MSA), a rare and deadly neurodegenerative disease, which progression is almost exclusively described by measurement scales [20].

2 Methodology

Let consider a cohort of individuals followed-up over time. We define the time of occurence of an event of interest. This event may be due to different causes, defining a competing risk setting. As some individuals may be censored before this occurence, we observe where is the censoring time, and the event indicator such that when the event of cause occurs first and before censoring (), and otherwise. We also collect repeated measures of markers measuring the same construct of interest (for instance items of a measurement scale or only 1 isolated item, or 1 sumscore). The marker values are noted for subject () and marker () at time with the repeated occasion (). The event of interest truncates the collection of the repeated markers so that necessarily . In this methodology, the markers do not need to be measured at the same visit times across markers and across individuals.

The joint model aims to analyze the trajectory of the construct of interest over time and the risk of event by accounting for their correlation. As shown in Figure 1, the joint model includes two submodels, one for the longitudinal process (on the left) and one for the time-to-event process (on the right). They are detailed below.

[width = 15cm, trim= 0cm 7cm 15cm 0cm, clip, page = 1]schema_methodo.pdf

Figure 1: Schematic diagram of the joint model structure for repeated markers measuring the same construct and a time of event (possibly of multiple causes)

2.1 Longitudinal model

When analyzing data from measurement scales or questionnaires, and more generally psychometric data, one makes the distinction between the construct of interest which is a latent process, denoted and defined in continuous time (), and its noisy observations, denoted and measured at discrete time visits . As emphasized in Figure 1, the longitudinal model thus consists of a structural model describing the latent process over time, and measurement models defining its link with each marker [16].

2.1.1 Structural model

The trajectory over time of

can be described using a linear mixed model

[21] as follows :


with and

two vectors of covariates associated to

the vector of fixed effect parameters, and the

-vector of normally distributed individual random effects (

), respectively. We note that as the quantity described here is not directly observed or measured, we do not consider any measurement error in the equation. The term defines the mean trajectory of the latent process at the population level while captures the individual deviation to the mean trajectory.

2.1.2 Measurement model

The measurement models define the nature of the link between the latent process and each marker, considered as a noisy measure of . Depending on the study to be conducted and the data available, the markers may vary in number and nature (e.g., continuous, ordinal). We review here different measurement models according to the nature of the markers. In any case, one measurement model is to be defined for each marker .

Continuous Gaussian marker

Classically, continuous markers are considered as having a multivariate normal distribution with Gaussian errors. In this case, the measurement model for marker is:


where and are two marker--specific parameters transforming the continuous marker into the latent process scale, and is the independent Gaussian measurement error ().

In the case of one single Gaussian marker (), and reduces the equation to the standard definition of a linear mixed model:


Continuous non-Gaussian marker

When analyzing data from measurement scales, the sum-score of all or a part of the items is frequently considered as a relevant marker of the underlying process of interest. This score is often treated as continuous because of its large number of possible levels. However its distribution is very rarely Gaussian due to ceiling and floor effects induced by its lower and upper bounds, but also due to the fact that, as based on the sum of ordinal items, a one-point change does not necessarily have the same meaning depending on the level of the score. This phenomenom, called curvilinearity, is not compatible with the assumptions of the linear model and linear mixed model [22, 19]. In this context, a continuous but non-gaussian marker can be described using a curvilinear measurement model where a parameterized transformation normalizes marker :


with and the marker--specific parametric link function. This function is monotonically increasing and is often approximated by splines [16]. We can notice that if is a linear function, this model turns into the classical linear mixed model defined in the previous subsection.

Discrete marker

Ordinal items are frequently considered for measuring complex constructs. They can be described using proportional odds logistic models or cumulative probit models

[3, 23]. Such models are based on the assumption of increasing monotonicity indicating that a higher level of the item reflects a higher degradation of the latent dimension that it measures. This is formalized as follows:


where is an additional noise either Logistic in the proportional odds logistic model or Gaussian in the cumulative probit [3]. We consider here the latter with .


are to be estimated. Also called thresholds or locations, they correspond to the level of the latent construct at which the probability of observing

lower/higher is 0.5. These measurement models are usual in Item Response Theory (IRT) where they are called Graded Response Models (GRM) [24].

2.1.3 Caution note

In the case of multiple markers, this framework relies on important assumptions:

  • unidimensionnality of the latent construct. The model assumes that all markers should measure the same phenomenon, the latent dimension of interest. In the case of a multidimensional measurement scale, items should be grouped by unidimensional underlying dimensions, and different dimensions analyzed separately.

  • inter-marker and intra-marker conditional independence given the latent construct. The latent construct should constitute the only source of shared information across pairs of markers and within repeated measures of the same marker. In the case of residual dependency across markers, redundant markers may be deleted or their correlation modeled in the measurement models. In the case of residual dependency within marker, additional random-effect specific to the marker and/or marker-specific association with covariate may be added [16].

These assumptions should be verified in preliminary analyses (see for instance [25]) to prevent any bias in the results. Note that intra-marker conditional independence given the latent construct may be assessed in sensitivity analyses.

2.1.4 Identifiability Constraints

With exception for the specific case of equation (3), the longitudinal submodel as defined in this section is not identifiable. Due to the introduction of a latent process, identifiability constraints must be added. In this work, we chose to center and reduce the latent process at time 0 in the reference category. This translates into an intercept fixed to

(location constraint) and a variance of the first random-effect (usually the random intercept) fixed to

(dispersion constraint) [3]. This implies that the unit of the latent process corresponds to the residual inter-individual standard deviation of the latent quantity when time .

2.2 Survival model

The survival model describes the risk of occurrence of each cause of event, usually using a cause-specific proportional hazard model [3]. In shared random effect joint models, the instantaneous risk function depends on a (possibly time-dependent) function of the elements of the structural model (1) noted . This element, which notably depends on the random effects , is added as a linear predictor in the survival model to quantify the association/intensity of the link between the latent dimension and the considered event .

The cause-specific hazard model is described in continuous time () for each cause and each subject () as follows:


with the parametric baseline risk function associated to the vector of parameters for cause . This function is usually chosen among Weibull hazards, piecewise constant hazards or approximated using spline functions. The vector of exogenous time-independent covariates is associated with the vector of parameters , and the function is associated with the vector of parameters which quantifies the association between the two processes. Many different functions can be specified [26, 7]. We focus here on two structures only:

  • the vector of random effects (dimension ): the survival model is adjusted on the individual deviations to the mean latent process trajectory with

  • the current latent process level (dimension ): the instantaneous risk is adjusted on the level of the underlying process (defined in Equation 1) at the same time with


2.3 Estimation

Let denotes the vector of all sub-model parameters including :

  • for the structural model: vector of fixed effects, parameters of the random-effects variance-covariance matrix ;

  • for the measurement model: for all , the vector of parameters for marker , it can be regression parameters if marker is Gaussian, link function parameters if marker is continuous but not necessarily Gaussian, or thresholds if marker is ordinal;

  • for the survival model: for all , the vector of parameters in the baseline hazard function for cause , vector of parameters for all covariates included in the survival model, and , the association parameters.

The joint model parameters can be estimated by maximizing the log-likelihood where is the individual contribution to the likelihood.

2.3.1 Individual contribution to the likelihood

Thanks to the assumption of independence between the latent dimension , measured by the markers , and the time-to-event conditionnally to the random effects , the individual contribution to the likelihood for subject can be developed as follows:


with the generic notation for a density function.

The density of the repeated marker depends on its nature. The density of a continuous marker can be expressed using the change of variable theorem:


where designates the value of marker transformed by its link function. Functions and are the density function of a normal distribution and the jacobian associated to the link function, respectively.

For a discrete marker, it can be expressed as follows:


where is the standard normal distribution function.

The density of the time-to-event is:


with the hazard function specific to cause given in Equation 6 and the survival function combining all causes:


Finally, the density of the random effects is the density of a centered multivariate normal distribution:


In the presence of delayed entry, that is when a subject is included in the study some time after time zero, the time-to-event is left-truncated. This is accounted for by considering instead the individual contribution to the likelihood for truncated data in which the naive individual contribution is divided by the probability to survive until the entry time :


2.3.2 Software

This model is implemented on the R package JLPM ( The optimization of the log-likelihood is numerically carried out by a Marquardt-Levenberg algorithm with stringent convergence criteria based on the first and second derivatives of the log-likelihood (see [27] for details). Example scripts are given in the package.

The integral on the random effects in Equation 9 does not have analytical solution so that the computation of the likelihood involves numerical integrations. This can be carried out with a Monte-Carlo algorithm which consists in simulating many draws of the random effects (generally around 100000) and to compute and then average the results of the function to integrate. In this package, this is carried out with a quasi Monte-Carlo algorithm. This integration method uses a deterministic sequence rather than a random one. Then, smaller integration errors are obtained due to the low discrepancy of the chosen deterministic sequence, so that the number of required points is usually reduced to 1000.

Moreover, when the association structure is time-dependent (not the case for the dependence on the random effects), the survival function in Equation 13 also involves an integral with no analytical solution; it is approximated by a Gauss-Kronrod gaussian quadrature with 15 points [28].

3 Simulation study

A small simulation study was performed to illustrate the methodology and validate the estimation procedure implemented in the R-package JLPM.

3.1 Simulation design

We considered the setting of the time to a unique cause of event, and the repeated measures to 4 ordinal markers measuring the same underlying construct. The longitudinal part constitutes a dynamic Item Response Theory model for graded responses (see [23] in the same special issue). The longitudinal and survival processes were associated through the current latent process level (i.e., dependence structure defined in Equation (8)).

We simulated 500 samples of 300 subjects each. The structural model for the underlying process consisted in a linear function of time at the population () and individual level (). The time-to-event was defined by a Weibull baseline risk function with . The unique linear predictor was the current level of the underlying process with an association parameter fixed to .

Visit times were generated every year (or time unit) from year 0 and up to the minimum between year 4 (administrative censoring) and the time-to-event. This lead to 36.9% of censoring on average and 3.9 repeated measures of each marker on average. The 4 markers had 4 ordinal levels each. The marker data were generated according to equation (5) with thresholds , , and and measurement error variances fixed to .

3.2 Results

We report in Table 1 the mean estimate, relative bias, variance estimate (empirical or asymptotic) and the coverage rate of the 95% confidence interval for each parameter. The estimation procedure provides very good results on this example: for all parameters, the relative bias does not exceed 2%, the mean asymptotic variance is close to the empirical variance, and the coverate rate of the 95% confidence interval is very close to the nominal value.

parameters true mean relative empirical mean asymptotic 95%CI
value estimate bias standard standard coverage
(in %) deviation deviation rate (in %)
survival model
      baseline risk function
          Weibull scale 0.447 0.447 -0.1 0.013 0.014 95.0
          Weibull shape 2.236 2.243 0.3 0.112 0.113 94.8
      association parameter, 0.1 0.099 -0.6 0.061 0.065 96.0
structural model
      adjustment covariates,
          time 1 1.010 1.0 0.076 0.078 96.2
      random effect covariance parameters,
          choleski 1 0 0.005 - 0.053 0.052 95.0
          choleski 2 0.447 0.445 -0.5 0.052 0.048 91.4
outcome-specific measurement model
          threshold 1 for item 1 0.5 0.501 0.3 0.087 0.087 95.2
          threshold 2 for item 1 0.707 0.709 0.2 0.039 0.038 96.2
          threshold 3 for item 1 0.707 0.712 0.6 0.040 0.038 92.6
          threshold 1 for item 2 0.25 0.249 -0.6 0.084 0.085 96.2
          threshold 2 for item 2 0.707 0.709 0.3 0.036 0.039 96.2
          threshold 3 for item 2 0.224 0.222 -0.6 0.033 0.033 94.8
          threshold 1 for item 3 0.1 0.102 1.9 0.082 0.085 95.8
          threshold 2 for item 3 0.316 0.313 -0.9 0.038 0.036 93.6
          threshold 3 for item 3 0.447 0.449 0.4 0.036 0.036 94.4
          threshold 1 for item 4 0.2 0.200 0.1 0.083 0.085 97.4
          threshold 2 for item 4 0.447 0.449 0.3 0.035 0.036 95.2
          threshold 3 for item 4 0.632 0.634 0.3 0.038 0.037 94.4
      standard deviation from measurement error,
          sd for item 1 1 1.008 0.8 0.081 0.084 96.0
          sd for item 2 1 1.008 0.8 0.082 0.085 96.2
          sd for item 3 1 1.005 0.5 0.089 0.087 93.0
          sd for item 4 1 1.008 0.8 0.085 0.084 95.2
Table 1: Summary of the estimation on 500 replicates for a joint model with 4 repeated ordinal markers, and samples of 300 individuals.

4 Illustration in Multi-System Atrophy

We illustrate our methodology to the Multiple-System Atrophy (MSA), a rare neurodegenerative disease characterized by various combinations of parkinsonism, cerebellar ataxia and dysautonomic symptoms. The disease progresses very fast and is fatal with a median survival between 8 and 10 years after the first symptoms onset [20]. The occurrence of death suddenly interrupts the follow-up of MSA patients who are usually the most affected. This constitutes an informative dropout that needs to be accounted for in the statistical analyses to avoid biasing the model estimates.

As many other neurodegenerative diseases, MSA clinical progression is studied almost exclusively using scales that measure the motor degradation, the dysautonomic dysfunction or the pathology-related quality-of-life. In particular, the Unified Multiple-System Atrophy Rating Scale (UMSARS) assesses disease severity in MSA patients in 4 sub-scales. The first one, reported by the patient, assesses functional impairments. In this illustration, we focus on the 2nd item of the first UMSARS subscale which measures dysphagia, the difficulty to swallow food and saliva. Clinical assumptions are that, among the motor and functional impairments of MSA, dysphagia could represent an important aspect of the prognosis and have a substantial impact on the risk of death. This ability progressively deteriorates during the course of the disease and most severe cases requires nasogastric tube or gastrostomy feeding. Dysphagia is measured as a 5-level likert scale (from 0 Normal state to 4 Severe incapacity) leading to repeated ordinal data.

Finally, the MSA diagnosis can be made several years after the first symptoms since the first MSA symptoms may not be specific (e.g., some are similar to those a Parkinson’s disease). This delayed diagnosis yields to a delayed entry in the cohort.

By accounting for these three challenges (ordinal repeated data, informative dropout and delayed entry), our methodology is particularly useful to study MSA progression.

4.1 Multiple-System Atrophy French cohort

Since 2007, the university hospitals of Bordeaux and Toulouse, the two French reference centres for MSA, have constituted the MSA French cohort which includes all the patients diagnosed with a MSA, and follows them every year with a complete clinical examination including the UMSARS completion. For this application, we considered all the MSA cases included in the cohort who had at least one completed dysphagia item during the follow-up before the administrative censoring on December 31st 2019. A total of 634 MSA patients were included with a total of 1819 dysphagia assessments (see description in Table 2). Patients were included in the cohort 4.58 (SD=2.61) years on average after their first symptoms. They were 65.1 (SD=8.2) years old in mean at entry. Two subtypes of MSA are distinguished: MSA-P for predominant parkinsonism and MSA-C for predominant cerebellar impairment. MSA-P (66.1%) were more frequent than MSA-C in the cohort (33.9%) as reported in Caucasian cohort [29]. As for many neurodegenerative diseases, the definite diagnosis of MSA is based on post-mortem neuropathologies. In clinical practice, the MSA diagnosis is therefore based on criteria established in 2008 by Gilman et al. [30] and giving two degrees of diagnostic certainty: probable and possible. At first visit, 151 (23.8%) and 483 (76.2%) met consensus for possible MSA and probable MSA, respectively.

Patients had in mean 2.9 (SD=2.1) clinical visits, and remained in the cohort up to 6.9 (SD=3.3) years after first symptoms onset on average. More than half of the patients (51.89%) died and 17.82% were not seen in the 18 months preceding administrative censoring. As shown in Figure 2 (bottom right), the survival probability dropped very quickly from approximately 2.5 years after the first symptom onset reaching a survival probability of 0.38 [0.33 - 0.43] after 10 years, 0.094 [0.061 - 0.144] after 15 years and 0 after 20 years. This illustrates the importance of the truncation of the follow-up by the occurence of death in MSA. Previous analyses suggest that this truncation is related to the course of the disease [20].

Regarding dysphagia, at baseline/entry, 180 (28.39%) patients had a normal state (no particular difficulty to swallow), 242 (38.17%) a mild impairment (choking less than once a week), 139 (21.92%) a moderate impairment (occasional food aspiration with choking more than once a week), 62 (9.78%) a marked impairment (frequent food aspiration) and 6 (0.95%) a severe state requiring assistance to be feeded. Dysphagia item degradation substantially progressed over disease time as illustrated with the individual observed trajectories (Figure 2 -bottom left).

Characteristic N (%) mean sd
At baseline
          Male 311 (49.05%)
          Female 323 (50.95%)
          Bordeaux 320 (50.47%)
          Toulouse 314 (49.53%)
          MSA-C, with predominant cerebellar impairment 215 (33.91%)
          MSA-P, with predominant parkinsonism 419 (66.09%)
      Diagnosis certainty
          Possible 151 (23.82%)
          Probable 483 (76.18%)
      Age at first symptom onset 60.48 8.32
      Age at entry 65.06 8.19
      Years since first symptom onset 4.58 2.61
      Dysphagia item 1.16 0.98
          0. Normal 180 (28.39%)
          1. Mild impairment 242 (38.17%)
          2. Moderate impairment 139 (21.92%)
          3. Marked impairment 62 (9.78%)
          4. Nasogastric tube or gastrostomy feeding 6 (0.95%)
During follow-up
      Visits with dysphagia item completed 1819 (98%)
      Dysphagia item observations per patient 2.87 2.09
      Years of follow-up 6.94 3.34
      Patients with nasogastric tube or gastrostomy feeding 31 (4.89%)
      Drop-out 113 (17.82%)
      Death 329 (51.89%)
Table 2: Description of the MSA sample at baseline and over time (N=634)

4.2 Specification of the model

This study aimed at better understanding progression of swallowing difficulties over disease time and at quantifying its association with death. To do so, we built a shared random effect joint model as illustrated in Figure 2, and considered time since first symptoms as the study timescale to describe the natural history of the item progression.

[width = 14cm, trim= 0cm 0cm 8.5cm 0cm, clip, page = 2]schema_methodo.pdf

Figure 2: Schematic diagram of the shared random effect joint model applied to MSA dysphagia study, along with observed individual observed trajectory of dysphagia, and Kaplan-Meier estimate of survival probability (N=634)

The observed 5-level item probability was linked to the underlying continuous dysphagia process at the exact same time by a cumulative probit model (cf. Equation 5). The trajectory of the underlying dysphagia process was simultaneously described over time by a latent linear mixed model (cf. Equation 1). We considered a quadratic trajectory to take into account a possible nonlinear evolution, with individual random effects on the intercept, slope and quadratic slope to take into account the within patient correlation. The association between dysphagia and death was simultaneously modelled using a survival model where the instantaneous risk of death was a function of the current underlying level of dysphagia (cf. Equation 6 with option 8). The baseline hazard function was defined by cubic M-splines with knots placed at , , , and years. The structural longitudinal model and the survival model were adjusted for sex (male or female), diagnosis (MSA-C or MSA-P), certainty of the diagnosis (possible or probable) and age at the first symptoms onset.

4.3 Results

The estimates of the joint model are provided in supplementary table.

The predicted trajectories of dysphagia item are displayed in Figure 3 according to profiles of patient differing by the four main covariates of interest (sex, age at first symptom onset, MSA subtype and certainty degree of diagnosis). The reference profile corresponded to a male patient, 60 years old at first symptoms onset, diagnosed MSA-C with possible certainty (black line).

Overall, the trajectories of dysphagia did not significantly differ according to sex, age and MSA subtype. Male and female patients seemed to endure the same dysphagia progression over at least the first 5 years of the disease. Female patients seemed to face a more important degradation afterwards although the dysphagia progression was not statistically different by sex (interaction on time and quadratic time, p=0.407). Dysphagia evolved slighly faster for older patients with substantial differences during the follow-up. At 5 years of disease, the dysphagia level was 0.51 (95%CI=1.29,1.75), 1.27 (95%CI=1.06,1.50) and 1.04 (95%CI=0.83,1.26) for patients aged 70, 60 and 50 at the beginning of the disease, respectively. This highlights an approximate gap of 0.5 points between patients with 20 years difference, stable between 5 and 10 years of disease.

No difference in dysphagia progression was observed between MSA-C and MSA-P patients. At the beginning of the disease, MSA-P patients seemed to have a little higher level of dysphagia than MSA-C patients, but after 5 years of disease, this trend reversed and MSA-C patients surpassed MSA-P patients in terms of degradation degree. This could be explained by the fact that MSA-P patients are initially more affected on the motor level as opposed to MSA-C patients who suffer from non-motor and dysautonomia symptoms [20].

Trajectories of dysphagia substantially differed according to the diagnosis certainty. Patients diagnosed with possible MSA had a significantly higher level of dysphagia at the beginning of the disease compared to patients diagnosed with more confidence (probable) with a mean level of 0.59 (95%CI=0.37,0.90) and 0.36 (95%CI=0.20,0.60) at onset for possible and probable, respectively. However, patients with probable MSA progressed much more rapidly reaching a mean level of 2.78 (95%CI= 2.49,3.02) after only 10 years while possible MSA patients still had in mean 2.17 (95%CI=1.85,2.49).


Figure 3: Mean trajectories of dysphagia item predicted by the shared random effect joint model according to the four main covariates (sex, age at first symptom onset, MSA subtype and certainty degree of diagnosis). The reference profile is a male patient, 60 years old at first symptoms onset, diagnosed MSA-C with possible certainty. Shades represent the 95% confidence intervals obtained by Monte Carlo approximation with 1000 draws. The reported p-values are those of Wald tests for the association at inclusion and the association with the functions of time.

The trajectory of dysphagia was also substantially associated with the risk of death through its current underlying level (HR = 3.96, CI95% =2.75,5.70) after adjustment on sex, age at first symptom onset, MSA subtype and certainty degree of diagnosis (Table 3). This means that the level of dysphagia at a certain time is strongly associated with the instantaneous risk of death at the exact same time: the higher the level of dysphagia, the higher the risk of death. Taking into account the association between progression of dysphagia and risk of death substantially changed the association estimates for sex and diagnosis certainty as shown in Table 3 which reports the estimates of the proportional hazard model accounting for dysphagia or not accounting for it (i.e., fixing its association to ). The risk of death was further reduced for women compared to men after taking into account the underlying level of dysphagia (HR = 0.82, CI95% = 0.65,1.04) and the higher risk of death of probable cases was slightly attenuated when accounting for the level of dysphagia impairment (HR = 1.40, CI95% = 1.03,1.90). This was expected as we saw that probable patients had a significantly higher level of dysphagia than possible patients (Figure 3). In contrast, the association with age at first symptoms was not impacted by the adjustment on dysphagia suggesting that the effect of age on the risk of death in this pathology may be independent of the level of clinical progression. Older patients at first symptoms had a higher risk of death (HR = 1.31, 95%CI = 1.13,1.52 for a 10-year difference). Finally, with or without adjustment for dysphagia, there was no evidence of difference of death risk according to MSA subtype (HR = 1.07, 95%CI = 0.83,1.38).

Hazard ratios [95% CI]
Covariate, modality (reference modality) Without dysphagia With dysphagia
    Women (ref Men) 0.923 [ 0.742 ; 1.147 ] 0.821 [ 0.648 ; 1.04 ]
Age at 1st symptom onset,
    10 years gap (ref 60 years old) 1.377 [ 1.209 ; 1.568 ] 1.309 [ 1.131 ; 1.515 ]
    MSA-P (ref MSA-C) 1.031 [ 0.817 ; 1.3 ] 1.073 [ 0.834 ; 1.379 ]
Diagnosis certainty,
    Probable MSA (ref Possible MSA) 1.869 [ 1.401 ; 2.494 ] 1.398 [ 1.027 ; 1.902 ]
Current level of dysphagia,
    unit = 1 SD at first symptoms - 3.962 [ 2.752 ; 5.702 ]
Table 3: Hazard ratios for the risk of death in proportional hazard models when accounting or not for the progression of dysphagia over time (N=634)

5 Concluding remarks

The joint modelling of longitudinal data and time-to-event data has become a standard methodology to address the problem of informative dropout and more generally investigate the association between longitudinal markers and clinical events [31]. In this work, we have shown how this methodology could be extended to longitudinal data of different natures by separating the longitudinal model for the quantity of interest defined as a latent process from the measurement model that links the quantity of interest to its repeated observations. This model thus naturally handles univariate or multivariate markers of the same underlying quantity but also continuous (Gaussian or not) or discrete markers. Specific examples include the case of a unique continuous Gaussian marker which corresponds to the classical joint model framework [7] or the case of various ordinal items measuring the same underlying construct. In the latter case, the longitudinal submodel is a dynamic Item Response Theory (IRT) submodel (see [23] in this special issue, [32]).

All joint models assume a given structure of correlation between the longitudinal data and the time-to-event data. In the literature, two main structure types can be found, the ”shared random effect joint models” which consider that the random effects from the longitudinal model capture all the correlation between the longitudinal process and the time-to-event process, and the ”joint latent class models” which assume that a latent group structure captures all the correlation between the two processes [9]. In this work, we focused only on the extension of the shared random effect joint model which enables a quantification of the association with hazard ratios, and have been favored for treating informative dropout.

Within the framework of shared random effect joint models, we considered a dependence structure on the current level of the latent process or on the individual random-effects of the structural longitudinal model. However, other dependence structures could also be worth exploring such as the current slope of the latent process, an accumulation of the latent process over time [33] or nonlinear functions of them [26]. The choice of the dependence structure depends on the application framework and is part of the model specification and building.

Although not illustrated in this work, the joint model handles competing risks of events, a situation which is often encountered in cohorts when different types of clinical progression and/or dropout may occur (e.g., types of recurrence in cancer research, multiple causes of death).

To conclude, measurement scales and questionnaires become more and more central in health studies to measure the physical, mental or psychological state of patients. With this joint model implemented in the R-package JLPM, we provide a new and relevant solution to study their trajectories over time and their association with events of interest.

Funding This work was funded by the French National Research Agency 455 (Project DyMES - ANR-18-C36-0004-01).


  • [1] A. A. Tsiatis, M. Davidian, Joint modeling of longitudinal and time-to-event data : an overview, Statistica Sinica 14 (3) (2004) 809–834, publisher: Institute of Statistical Science, Academia Sinica.
  • [2] L. Ferrer, V. Rondeau, J. Dignam, T. Pickles, H. Jacqmin-Gadda, C. Proust-Lima, Joint modelling of longitudinal and multi-state processes: application to clinical progressions in prostate cancer, Statistics in Medicine 35 (22) (2016) 3933–3948. doi:10.1002/sim.6972.
  • [3] D. Commenges, H. Jacqmin-Gadda, Dynamical Biostatistical Models, Chapman and Hall/CRC, New York, 2015. doi:10.1201/b19109.
  • [4] C. Thomadakis, L. Meligkotsidou, N. Pantazis, G. Touloumi, Longitudinal and time-to-drop-out joint models can lead to seriously biased estimates when the drop-out mechanism is at random, Biometrics 75 (1) (2019) 58–68. doi:10.1111/biom.12986.
  • [5] R. J. A. Little, Modeling the Drop-Out Mechanism in Repeated-Measures Studies, Journal of the American Statistical Association 90 (431) (1995) 1112–1121, publisher: [American Statistical Association, Taylor & Francis, Ltd.]. doi:10.2307/2291350.
  • [6] A. Rouanet, C. Helmer, J.-F. Dartigues, H. Jacqmin-Gadda, Interpretation of mixed models and marginal models with cohort attrition due to death and drop-out, Statistical Methods in Medical Research 28 (2) (2019) 343–356. doi:10.1177/0962280217723675.
  • [7] D. Rizopoulos, Joint Models for Longitudinal and Time-to-Event Data. With Applications in R, 2012.
  • [8] R. Henderson, P. Diggle, A. Dobson, Joint modelling of longitudinal measurements and event time data, Biostatistics (Oxford, England) 1 (4) (2000) 465–480. doi:10.1093/biostatistics/1.4.465.
  • [9] C. Proust-Lima, M. Séne, J. M. G. Taylor, H. Jacqmin-Gadda, Joint latent class models for longitudinal and time-to-event data: a review, Statistical Methods in Medical Research 23 (1) (2014) 74–90. doi:10.1177/0962280212445839.
  • [10] B. He, S. Luo, Joint modeling of multivariate longitudinal measurements and survival data with applications to Parkinson’s disease, Statistical Methods in Medical Research 25 (4) (2016) 1346–1358. doi:10.1177/0962280213480877.
  • [11] J. Wang, S. Luo, Joint modeling of multiple repeated measures and survival data using multidimensional latent trait linear mixed model, Statistical Methods in Medical Research 28 (10-11) (2019) 3392–3403. doi:10.1177/0962280218802300.
  • [12] M. Alsefri, M. Sudell, M. García-Fiñana, R. Kolamunnage-Dona, Bayesian joint modelling of longitudinal and time to event data: a methodological review, BMC Medical Research Methodology 20 (1) (2020) 94. doi:10.1186/s12874-020-00976-2.
  • [13] G. L. Hickey, P. Philipson, A. Jorgensen, R. Kolamunnage-Dona, Joint Models of Longitudinal and Time-to-Event Data with More Than One Event Time Outcome: A Review, The International Journal of Biostatistics 14 (1), number: 1 (2018). doi:10.1515/ijb-2017-0047.
  • [14] P. Fayers, A. Bottomley, EORTC Quality of Life Group, Quality of Life Unit, Quality of life research within the EORTC-the EORTC QLQ-C30. European Organisation for Research and Treatment of Cancer, European Journal of Cancer (Oxford, England: 1990) 38 Suppl 4 (2002) S125–133. doi:10.1016/s0959-8049(01)00448-8.
  • [15] A. Edjolo, C. Proust-Lima, F. Delva, J.-F. Dartigues, K. Pérès, Natural History of Dependency in the Elderly: A 24-Year Population-Based Study Using a Longitudinal Item Response Theory Model, American Journal of Epidemiology 183 (4) (2016) 277–285, number: 4. doi:10.1093/aje/kwv223.
  • [16] C. Proust-Lima, H. Amieva, H. Jacqmin-Gadda, Analysis of multivariate mixed longitudinal data: a flexible latent process approach, The British Journal of Mathematical and Statistical Psychology 66 (3) (2013) 470–487, number: 3. doi:10.1111/bmsp.12000.
  • [17] C. James, M. Powell, A. Seixas, A. Bateman, S. Pengpid, K. Peltzer, Exploring the psychometric properties of the CES-D-10 and its practicality in detecting depressive symptomatology in 27 low- and middle-income countries, International Journal of Psychology: Journal International De Psychologie 55 (3) (2020) 435–445. doi:10.1002/ijop.12613.
  • [18] P. M. Kruyen, W. H. M. Emons, K. Sijtsma, Shortening the S-STAI: consequences for research and clinical practice, Journal of Psychosomatic Research 75 (2) (2013) 167–172. doi:10.1016/j.jpsychores.2013.03.013.
  • [19]

    C. Proust-Lima, V. Philipps, J.-F. Dartigues, D. A. Bennett, M. M. Glymour, H. Jacqmin-Gadda, C. Samieri, Are latent variable models preferable to composite score approaches when assessing risk factors of change? Evaluation of type-I error and statistical power in longitudinal cognitive studies, Statistical Methods in Medical Research 28 (7) (2019) 1942–1957.

  • [20] A. Foubert-Samier, A. Pavy-Le Traon, F. Guillet, M. Le-Goff, C. Helmer, F. Tison, O. Rascol, C. Proust-Lima, W. G. Meissner, Disease progression and prognostic factors in multiple system atrophy: A prospective cohort study, Neurobiology of Disease 139 (2020) 104813. doi:10.1016/j.nbd.2020.104813.
  • [21] N. M. Laird, J. H. Ware, Random-effects models for longitudinal data, Biometrics 38 (4) (1982) 963–974.
  • [22] C. Proust-Lima, J.-F. Dartigues, H. Jacqmin-Gadda, Misuse of the linear mixed model when evaluating risk factors of cognitive decline, American Journal of Epidemiology 174 (9) (2011) 1077–1088. doi:10.1093/aje/kwr243.
  • [23] C. Proust-Lima, V. Philipps, B. Perrot, M. Blanchin, V. Sébille, Continuous-time modeling of self-reported outcome data: a dynamic Item Response Theory model, Methods (2021).
  • [24] F. B. Baker, S.-H. Kim, Item Response Theory: Parameter Estimation Techniques, Second Edition, CRC Press, 2004, google-Books-ID: gUdZDwAAQBAJ.
  • [25] B. B. Reeve, R. D. Hays, J. B. Bjorner, K. F. Cook, P. K. Crane, J. A. Teresi, D. Thissen, D. A. Revicki, D. J. Weiss, R. K. Hambleton, H. Liu, R. Gershon, S. P. Reise, J.-s. Lai, D. Cella, PROMIS Cooperative Group, Psychometric evaluation and calibration of health-related quality of life item banks: plans for the Patient-Reported Outcomes Measurement Information System (PROMIS), Medical Care 45 (5 Suppl 1) (2007) S22–31. doi:10.1097/01.mlr.0000250483.85507.04.
  • [26] M. Sène, C. Bellera, C. Proust-Lima, Shared random-effect models for the joint analysis of longitudinal and time-to-event data: Application to the prediction of prostate cancer recurrence, Journal de la Sociètè Française de Statistique 155 (Oct. 2013).
  • [27] V. Philipps, B. Hejblum, M. Prague, D. Commenges, C. Proust-Lima, Robust and efficient optimization using a Marquardt-Levenberg algorithm with R package marqLevAlg, RJournal (2021).
  • [28] P. Gonnet, A review of error estimation in adaptive quadrature, ACM Computing Surveys 44 (4) (2012) 22:1–22:36. doi:10.1145/2333112.2333117.
  • [29] G. K. Wenning, F. Geser, F. Krismer, K. Seppi, S. Duerr, S. Boesch, M. Köllensperger, G. Goebel, K. P. Pfeiffer, P. Barone, M. T. Pellecchia, N. P. Quinn, V. Koukouni, C. J. Fowler, A. Schrag, C. J. Mathias, N. Giladi, T. Gurevich, E. Dupont, K. Ostergaard, C. F. Nilsson, H. Widner, W. Oertel, K. M. Eggert, A. Albanese, F. del Sorbo, E. Tolosa, A. Cardozo, G. Deuschl, H. Hellriegel, T. Klockgether, R. Dodel, C. Sampaio, M. Coelho, R. Djaldetti, E. Melamed, T. Gasser, C. Kamm, G. Meco, C. Colosimo, O. Rascol, W. G. Meissner, F. Tison, W. Poewe, European Multiple System Atrophy Study Group, The natural history of multiple system atrophy: a prospective European cohort study, The Lancet. Neurology 12 (3) (2013) 264–274, number: 3. doi:10.1016/S1474-4422(12)70327-7.
  • [30] S. Gilman, G. K. Wenning, P. A. Low, D. J. Brooks, C. J. Mathias, J. Q. Trojanowski, N. W. Wood, C. Colosimo, A. Dürr, C. J. Fowler, H. Kaufmann, T. Klockgether, A. Lees, W. Poewe, N. Quinn, T. Revesz, D. Robertson, P. Sandroni, K. Seppi, M. Vidailhet, Second consensus statement on the diagnosis of multiple system atrophy, Neurology 71 (9) (2008) 670–676. doi:10.1212/01.wnl.0000324625.00404.15.
  • [31] D. Rizopoulos, E. Lesaffre, Introduction to the special issue on joint modelling techniques, Statistical Methods in Medical Research 23 (1) (2014) 3–10, number: 1. doi:10.1177/0962280212445800.
  • [32] A. Barbieri, Méthodes longitudinales pour l’analyse de la qualité de vie relative à la santé en cancérologie, phdthesis, Université Montpellier (Jun. 2016).
  • [33] E.-R. Andrinopoulou, D. Rizopoulos, Bayesian shrinkage approach for a joint model of longitudinal and survival outcomes assuming different association structures, Statistics in Medicine 35 (26) (2016) 4813–4823. doi:10.1002/sim.7027.