What is joint modelling? Longitudinal and survival data

Joint modelling fits two models at once: one for a measurement taken repeatedly over time, such as blood pressure, prostate-specific antigen (PSA) or a quality-of-life score, and one for the time to an event, such as death, relapse or a transplant. The two are linked, so that the model describes how the risk of the event depends on the underlying course of the measurement, and allows for the fact that the measurements stop when the event happens. Analysed separately, the measurement error and the dropout can bias both analyses.

The most common form combines a linear mixed effects model with a proportional hazards survival model, linked through the random effects that describe each person’s underlying trajectory (Crowther et al., Stat Med 2016). Ibrahim, Chu and Chen (J Clin Oncol 2010) give an introductory overview aimed at clinical trials.

Schematic, simulated for illustration. One patient in a joint model. Top: the biomarker is measured with error (dots) around the patient’s true trajectory (line), and the measurements end with the event. Bottom: the hazard follows the trajectory, here with a positive association, so the risk rises as the biomarker rises.

The two submodels

For person $i$ , the longitudinal submodel says that the value observed at time $t$ is a true value plus measurement error:

y_{i} (t) = m_{i} (t) + ϵ_{i} (t),

m_{i} (t) = x_{i} (t)^{⊤} β + z_{i} (t)^{⊤} b_{i},

b_{i} ~ N (0, Σ) .

The trajectory $m_{i} (t)$ is the underlying value of the biomarker. It combines fixed effects $β$ that everyone shares with random effects $b_{i}$ that belong to person $i$ , such as their own starting level and rate of change, with $x_{i} (t)$ and $z_{i} (t)$ functions of time and covariates and $Σ$ the covariance matrix of the random effects. The measurement error $ϵ_{i} (t)$ has mean zero and constant variance $σ^{2}$ . The survival submodel then makes the hazard of the event depend on that trajectory:

h_{i} (t) = h_{0} (t) e^{γ^{⊤} v_{i} + α m_{i} (t)} .

Here $h_{0} (t)$ is the baseline hazard, the hazard when the true value and the covariates are all zero, $v_{i}$ are baseline covariates such as treatment, with log hazard ratios $γ$ , and $α$ is the association parameter: the log hazard ratio for a one-unit higher true current value of the biomarker, among people with the same baseline covariates. If the biomarker is systolic blood pressure in mmHg and $α = 0.03$ , a true value 1 mmHg higher goes with a hazard ratio of 1.03, and one 10 mmHg higher with $\exp (0.3) = 1.35$ . It is the hazard ratio for the true value, not for the noisy measurement.

The model assumes that, given the random effects and the covariates, the measurements are independent of one another and of the event time, and, in its basic form, that the error and the random effects are normally distributed. That suits a continuous biomarker whose errors are roughly symmetric. A skewed one such as PSA is usually modelled on the log scale, and, with natural logarithms, $α$ is then per unit of log PSA, that is per 2.7-fold change in PSA; a binary or count biomarker needs a generalised longitudinal submodel. It also assumes that the times at which the biomarker is measured carry no information about the outcome beyond the model; where patients are seen more often when they are sicker, as in health-care records, that can fail (Gasparini et al., Stat Neerl 2020, who study the problem for longitudinal outcomes). With the trajectory in the hazard, the coefficient for treatment is its effect at the same true biomarker value, not its total effect on survival.

Why not fit the two models separately?

Each shortcut goes wrong in its own way.

Measurement error. A biomarker entered into a survival model as measured is the true value plus noise, so its effect is typically biased towards no effect, and the last observed value may be out of date. Sweeting and Thompson compared a joint model with a naive proportional hazards model with a time-dependent covariate, and with a two-stage approach that puts fitted values from a mixed model into a survival model. They showed that such approximations can severely underestimate the association between the current underlying value and the hazard (Biom J 2011).
A biomarker that stops at death. For a terminal event such as death, the biomarker can only be measured while the person is alive, so it is an internal time-varying covariate (Kalbfleisch and Prentice, 2002), and its path carries information about when the event will happen. Its future path is not known in advance and stops at the event, so a Cox model with the last observed value as a time-dependent covariate gives a hazard ratio for that value but cannot, on its own, give a survival curve or a prediction for a new patient.
Dropout can be informative. Analysed on their own, the measurements stop for people who die, and death can be related to the course of the biomarker. In the Vasopressin and Septic Shock Trial, organ-failure scores compared between groups could be biased if death differed between them, and 28-day mortality was numerically higher with noradrenaline (39% against 35%) (Harhay et al., Crit Care Explor 2020).
Change scores in two steps. Estimating each person’s change first and using it as a fixed covariate in a survival model ignores the dropout and the error in that estimate, and a change taken over follow-up conditions on surviving to its end (immortal time bias). A joint model can describe the long-term trend while allowing for dropout and measurement error (Andersson et al., Am J Epidemiol 2017).

How the two models are linked

The link is the association structure, and the clinical question should choose it.

Current value. The hazard depends on the true value of the biomarker now, $m_{i} (t)$ , as above. In the scenarios simulated, estimates for this structure were little affected by an over-simple model for the trajectory (Crowther et al., Stat Med 2016), but the fit should still be checked.
Current value and rate of change. Add the current slope of the trajectory, for the question of whether a rising or falling biomarker matters beyond its level. The rate of change was highly sensitive to assuming a simpler trend than the true one, which can bias the estimates substantially (Crowther et al., Stat Med 2016).
Shared random effects. Link the hazard to a latent process shared with the measurements; the most familiar special case is the random effects themselves, such as a person’s own starting level and slope (Henderson, Diggle and Dobson, Biostatistics 2000). Linking to the intercept alone adjusts a baseline biomarker for measurement error, using repeated measurements that registry data often hold but analyses ignore (Crowther, Lambert and Abrams, BMC Med Res Methodol 2013).
Other forms, such as the cumulative effect, where the hazard depends on the area under the trajectory so far, are described in Rizopoulos, 2012 and offered, for example, by stan_jm.

What a joint model is used for

The effect of a biomarker on survival, with the measurement error and the intermittent measurement handled in the model.
The course of the biomarker itself, allowing for informative dropout such as death, under the model’s assumptions.
Dynamic prediction. For a new patient, the model gives the probability of staying event-free for a further period, given that they are still event-free at the latest measurement, and updates it with each new measurement (Rizopoulos, Biometrics 2011; Proust-Lima and Taylor, Biostatistics 2009, for recurrence of prostate cancer from repeated PSA). Landmarking is a simpler alternative that makes fewer assumptions about the biomarker. Simulation studies have found it less efficient than a correctly specified joint model, but the joint model can do worse when the model for the biomarker is wrong (reviewed by Putter and van Houwelingen, 2022).
Quality-adjusted life-years in health technology assessment, estimated in continuous time from a joint model of quality of life and survival (Crowther et al., Med Decis Making 2026). It is implemented in our Stata package jmqaly.

Fitting a joint model

The likelihood has no closed form for most joint models. It integrates over the random effects, and, when the hazard depends on the current value, over time to get the cumulative hazard. Gauss–Hermite quadrature handles the random effects and Gauss–Kronrod quadrature the integral over time, and fully adaptive quadrature reduced bias and computation time compared with non-adaptive quadrature in a simulation study (Crowther et al., Stat Med 2016). A flexible parametric survival submodel, with restricted cubic splines for the log baseline cumulative hazard, avoids committing to a Weibull or other fixed shape (Crowther, Abrams and Lambert, Stat Med 2012); when the hazard depends on the current value, the cumulative hazard still needs numerical integration. The Bayesian alternative fits the same kind of model by MCMC and carries the uncertainty in the estimates into predictions (Sweeting and Thompson, 2011), at the cost of choosing priors and checking convergence. Several random effects or several biomarkers enlarge the integration and slow the fit, and more than one type of event needs a hazard for each (see what are competing risks?).

In R, joint models can be fitted with JM (Rizopoulos, J Stat Softw 2010), JMbayes2, which takes several longitudinal outcomes and competing risks in a Bayesian framework, joineR and the stan_jm function of rstanarm. In Stata, stjm (Crowther, Abrams and Lambert, Stata J 2013) and our merlin command (Crowther, Stata J 2020) fit them; our introductory tutorial shows how, with code.

Common mistakes

A trajectory that is too simple. A straight line for a biomarker that curves biases the association structures that depend on its slope. Use polynomials or splines for the trajectory and check the fit against the data (Crowther et al., Stat Med 2016).
Choosing the association structure by habit. The current value and the rate of change answer different questions. Choose the one the clinical question needs before looking at the results, and report the others as sensitivity analyses.
Leaving out delayed entry. If people enter follow-up after the start of the time scale, for example patients already diagnosed when the study begins, the likelihood must condition on their having survived to entry, which needs a second set of numerical integration (Crowther et al., Stat Med 2016).
Modelling a skewed biomarker on its raw scale. The normal error in the basic model is a poor description of a biomarker such as PSA. Transform it, usually with the logarithm, and say what unit $α$ is per.
Not checking the integration. Increase the number of quadrature points until the estimates stop changing, and report that you did.
Reading $α$ as a causal effect. It is an association between the true current value and the hazard. Whether changing the biomarker would change the risk is a causal question that the model does not answer by itself, and the treatment coefficient is conditional on a trajectory that is measured after randomisation, so it is not protected by randomisation.
Using a joint model when a simpler analysis answers the question. With a single baseline measurement, a survival model will do (one measurement per person cannot separate error from true variation). A mixed model is enough if dropout depends only on the measurements already observed. For predictions at a few fixed times, landmarking may be enough.

Where joint models are used

Our team has published methods and applications of joint models since 2012, and a selection of the papers is on our joint models page. Some examples:

Breast cancer. Mammographic density and cause-specific mortality in 1,740 Swedish women, with 6,317 density measurements in the first five years of follow-up, including baseline (Andersson et al., Am J Epidemiol 2017).
Septic shock. Organ-failure scores in 763 patients in a randomised trial of vasopressin against noradrenaline (norepinephrine), with scores recorded until discharge, death or day 28 (Harhay et al., Crit Care Explor 2020).
Early pregnancy. Daily urinary human chorionic gonadotrophin (hCG) and miscarriage in 127 women, used to ask whether early pregnancy outcome could be monitored as measurements arrive (Ashra et al., Sci Rep 2020).
HIV. Viral load subject to a lower detection limit, modelled with a two-part hurdle submodel, and the duration of the first treatment regimen (Brilleman et al., Stat Med 2016).
Haemodialysis. Changes in body mass index (BMI) and the two competing events of transplant and death without transplant in 16,414 patients, with latent classes of BMI trajectory (Brilleman et al., Epidemiology 2019).
Prognostic biomarkers. Baseline systolic blood pressure and stroke in obese patients with type 2 diabetes, adjusting the blood pressure for measurement error (Crowther, Lambert and Abrams, 2013).

A repeated measurement in your own data?

Deciding whether a study needs a joint model, choosing how the biomarker and the event are linked, and explaining the result to clinicians and reviewers are the hard parts. We can help: see our joint modelling service, or talk to us about your data.

Want to learn more?

MethodsJoint longitudinal–survival models Longitudinal data analysis Survival analysis

What is joint modelling?