📊 statistics

A Bayesian location-scale joint model for time-to-event and multivariate longitudinal data with association based on within-individual variability

This paper proposes a novel Bayesian joint model that links time-to-event outcomes with multivariate longitudinal data by incorporating within-individual variability through a mixed-effect location-scale framework, thereby overcoming the limitations of summary statistics and regression dilution while demonstrating its efficacy through simulations and an application to cystic fibrosis mortality data.

Original authors: Marco Palma, Ruth H Keogh, Siobhán B Carr, Rhonda Szczesniak, David Taylor-Robinson, Angela M Wood, Graciela Muniz-Terrera, Jessica K Barrett

Published 2026-06-02

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Marco Palma, Ruth H Keogh, Siobhán B Carr, Rhonda Szczesniak, David Taylor-Robinson, Angela M Wood, Graciela Muniz-Terrera, Jessica K Barrett

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to predict how long a car will last before it breaks down. Traditionally, doctors and statisticians have looked at the average speed of the car over time. If the average speed is dropping, they assume the car is getting older and might break soon.

But this paper argues that looking at the average isn't enough. Sometimes, a car that drives at a steady 60 mph is fine, while another car that bounces wildly between 20 mph and 100 mph is in serious trouble—even if their average speed is the same. This "bouncing" or variability is what the authors call "within-individual variability."

Here is a breakdown of what the paper does, using simple analogies:

1. The Problem: The "Average" Trap

In the past, when studying diseases like Cystic Fibrosis (a condition that damages the lungs), researchers would take a patient's health data (like lung function or weight) and calculate a simple summary, such as the standard deviation (how much the numbers jump around).

The paper points out two big problems with this old way of doing things:

The "Snapshot" Problem: Simple summaries don't account for when the jumps happened. Did the patient fluctuate wildly last year and stabilize this year? A simple average misses that story.
The "Blurry Lens" Problem: When researchers used these simple summaries to predict death or disease progression, the results were often "blurred" or too weak. It's like trying to read a sign through a foggy window; you know something is there, but you can't see the details clearly. This is called "regression dilution."

2. The Solution: A "Double-Track" Model

The authors built a new mathematical tool (a "Joint Model") that looks at two tracks at the same time, rather than just one.

Track A (The Average): This tracks the patient's general health trend (e.g., is their lung function slowly going down?).
Track B (The Jitters): This tracks the variability or the "jitters" around that trend. It asks: "How much does this person's health swing up and down from their own average?"

Think of it like monitoring a heartbeat.

Track A is the heart rate (60 beats per minute).
Track B is the rhythm. Is it steady (steady 60), or is it erratic (60, then 40, then 90, then 50)?
The authors' model treats the "jitters" (Track B) as a real, important predictor of risk, not just noise.

3. How It Works: The "Shared Secret"

The model uses a clever trick called random effects. Imagine every patient has a hidden "personality" or "style" that influences both their average health and how much they fluctuate.

The model learns this hidden style for every single person.
It then uses this hidden style to predict the risk of a bad event (like death).
Crucially, it allows the "jitters" (variability) to talk directly to the prediction of death.

The authors wrote this model in a computer language called Stan, which is like a high-powered calculator that can solve these complex puzzles by simulating thousands of possible scenarios to find the most likely answer.

4. Testing the Theory: The Simulation

Before using real people, the authors tested their model with fake data (a simulation).

They created a world where "jitters" definitely caused problems.
They compared their new model against the old, standard models.
The Result: The old models missed the connection between "jitters" and the bad outcome. They were like a detective ignoring a smoking gun. The new model, however, correctly identified that the "jitters" were a major warning sign.

5. Real-World Test: Cystic Fibrosis Patients

The authors applied their model to real data from the UK Cystic Fibrosis registry, looking at 3,282 women. They tracked two main things:

Lung Function (FEV1): How much air they can blow out.
Body Mass Index (BMI): A measure of nutrition/weight.

What they found:

Lung Function: It wasn't just about having low lung function that was dangerous. It was also about having unstable lung function. Women whose lung numbers bounced around wildly were at a higher risk of death, even if their average wasn't terrible.
Weight (BMI): While having a low weight was dangerous, the bouncing around of weight numbers didn't show a clear link to death risk in this specific study.

The Bottom Line

This paper introduces a smarter way to look at health data. Instead of just asking, "How is the patient doing on average?", it asks, "How stable is the patient?"

By treating the "ups and downs" of a patient's health as a serious warning sign rather than just background noise, this model gives a clearer, more accurate picture of who is at risk. It proves that for some diseases, stability is just as important as the average.

Technical Summary: A Bayesian Location-Scale Joint Model for Time-to-Event and Multivariate Longitudinal Data

Problem Statement
In biomedical longitudinal studies, the within-individual variability (WIV) of health indicators is increasingly recognized as a critical predictor of disease progression and time-to-event outcomes, such as mortality. While simple summary statistics (e.g., standard deviation) are often used to quantify WIV, they suffer from significant limitations: they fail to account for time-varying changes, are sensitive to the number of visits, and can introduce biases such as "immortal time." Furthermore, when these summary statistics are used as covariates in survival models, the resulting estimates of hazard ratios are subject to regression dilution bias, leading to inaccurate associations. Existing joint modeling approaches have addressed WIV for single biomarkers or have used frequentist methods, but few have successfully integrated WIV for multivariate (correlated) longitudinal biomarkers into a Bayesian joint model framework where the variability itself drives the association with the event outcome.

Methodology
The authors propose a Bayesian Multivariate Joint Model with WIV-based association (JM-WIV). The model consists of two interconnected submodels:

Longitudinal Submodel (Mixed-Effect Location-Scale Model - MELSM):
- Instead of assuming constant variance, the model treats both the mean ( $\mu$ ) and the standard deviation ( $\sigma$ ) of $K$ longitudinal biomarkers as functions of time, covariates, and subject-specific random effects.
- The formulation is:
  $y_{ijk} = \eta^{\mu}_{ik}(t_{ij}) + \varepsilon_{ijk}, \quad g(\sigma_{ijk}) = \eta^{\sigma}_{ik}(t_{ij})$
  where $g(\cdot) = \log(\cdot)$ ensures non-negative variance.
- A shared random effect vector $\mathbf{b}_i$ links the mean and variability processes for each biomarker and allows for the modeling of correlations between different biomarkers (e.g., between the mean of biomarker A and the variability of biomarker B).
- The covariance structure is unstructured, allowing for complex correlations between random effects for means and standard deviations across multiple biomarkers.
Event Submodel (Proportional Hazard Model):
- The time-to-event is modeled using a Cox proportional hazard model with a flexible baseline hazard specified via B-splines.
- Crucially, the association term $\eta^{\alpha}_i$ incorporates information from the longitudinal submodel not just via current values of the means, but also via the current values of the standard deviations (variability).
- The authors define three potential association structures: functions of the current values (CV) of mean and SD, the linear predictors (LP) of mean and SD, or the random effects (RE) directly.
Implementation:
- The model is implemented in Stan using Hamiltonian Monte Carlo (HMC) for efficient parameter exploration.
- The likelihood is factorized into event and longitudinal processes conditional on random effects.
- Priors include normal distributions for fixed effects, half-t distributions for random effect standard deviations, and an LKJ prior for correlation matrices.

Key Contributions

Novel Framework: The paper introduces a Bayesian joint model capable of handling multivariate longitudinal data where the association with the time-to-event outcome is explicitly driven by the within-individual variability of the markers.
Addressing Regression Dilution: By modeling WIV directly within the joint framework rather than using two-stage summary statistics, the method aims to mitigate regression dilution bias.
Flexible Association: The model allows the hazard function to depend on the variability of biomarkers, a feature not present in standard joint models or previous location-scale extensions that did not include WIV in the survival submodel.
Software Availability: The authors provide an open-source implementation in Stan (via the rstanjmwiv package), extending existing joint modeling tools to handle bivariate (and potentially multivariate) location-scale settings.

Results

Simulation Study:
- A simulation study compared the proposed JM-WIV model against a standard joint model (implemented in JMbayes2) that assumes constant variance.
- When the true data-generating process included an association between WIV and the event, the standard joint model produced biased estimates for the survival coefficients and coverage probabilities well below the nominal 95% level.
- The JM-WIV model provided estimates with bias close to zero and coverage probabilities consistent with the nominal level, demonstrating its ability to correctly recover parameters when WIV is a relevant predictor.
Application to Cystic Fibrosis (CF) Data:
- The model was applied to data from the UK Cystic Fibrosis Registry (3,282 women, 29,094 visits) to evaluate the association between lung function (FEV1), malnutrition (BMI), and mortality.
- Longitudinal Findings: The model successfully captured the decline in FEV1 with age and the increase in BMI. It also quantified that FEV1 variability decreases with age, while BMI variability showed a specific association with F508 homozygosity.
- Survival Findings:
  - Lower mean FEV1 and lower mean BMI were associated with higher mortality risk (Hazard Ratios < 1).
  - Higher within-individual variability in FEV1 was significantly associated with increased mortality (HR = 1.374, 95% CI: 1.268–1.499).
  - No clear evidence was found for an association between BMI variability and mortality.
- The results confirm that fluctuations in lung function, independent of the average level, are a strong predictor of death in CF patients.

Significance and Claims
The paper claims that the JM-WIV model offers a methodological advancement by properly quantifying within-individual variability in a multivariate setting and linking it directly to survival outcomes without the biases inherent in summary-statistic approaches. The authors state that their simulation results demonstrate that neglecting the association of WIV with the event leads to bias and loss of coverage in survival coefficients. In the clinical application, the model confirms that higher lung function variability is a distinct risk factor for mortality, reinforcing the clinical utility of monitoring fluctuations rather than just mean trends.

The authors acknowledge limitations, noting that in location-scale models, WIV is indistinguishable from measurement error, and that the computational burden increases significantly with a large number of biomarkers or random effects. They suggest that future work could explore non-normal distributions, competing risks (e.g., transplant), and change-point structures to model the impact of new therapies on variability.