Joint modelling of time-dependent biomarker variability and time-to-event outcomes, a two-step approach

Technical Summary: Joint Modelling of Time-Dependent Biomarker Variability and Time-to-Event Outcomes

Problem Statement
Standard joint models for longitudinal biomarkers and time-to-event outcomes typically focus on the expected value (mean trajectory) of the biomarker, assuming constant residual variability across individuals and time. However, emerging evidence in oncology and cardiology suggests that the variability (instability) of longitudinal biomarkers carries independent prognostic information for clinical outcomes. While fully joint extensions exist that model subject-specific residual variance, they are computationally demanding, often require dedicated software packages, and lack the flexibility to handle multiple biomarkers, competing risks, or dynamic prediction. Furthermore, existing methods often rely on restrictive association structures or simple summary statistics (e.g., standard deviation over the whole period) that fail to capture time-dependent fluctuations, potentially introducing biases such as immortal time bias and regression dilution.

Methodology
The authors propose a flexible, computationally efficient two-step approach to incorporate time-dependent biomarker variability into standard joint modeling frameworks without requiring custom software extensions.

1. Step 1: Estimation of Subject-Specific Residuals
   A standard linear mixed-effects model is fitted to the longitudinal biomarker data to estimate the subject-specific mean trajectory ($m_i(t)$). This model assumes constant residual variance initially. The individual residuals ($\hat{\epsilon}_i(t)$) are calculated as the difference between observed values and the fitted subject-specific trajectories. These residuals represent the within-subject fluctuations around the expected trajectory.

2. Step 2: Modeling Variability as a Longitudinal Outcome
   The residuals (or their transformations, such as squared residuals $\hat{\epsilon}_i(t)^2$ or absolute residuals $|\hat{\epsilon}_i(t)|$) are treated as a new longitudinal outcome representing biomarker variability.

   • A second mixed-effects model is fitted to these residuals to capture both population-level trends and subject-specific patterns in variability.
   • These variability measures are then incorporated as time-varying covariates in a standard joint model alongside the original biomarker's mean trajectory.
   • The time-to-event submodel (typically a Cox proportional hazards model) links the event risk to both the current value of the biomarker mean trajectory and the current value of the estimated variability measure.

The approach is implemented using existing Bayesian joint modeling software (specifically the JMbayes2 package in R), utilizing Markov chain Monte Carlo (MCMC) sampling for parameter estimation.

Key Contributions

• Implementation Accessibility: The method allows researchers to model biomarker variability using standard joint modeling software, avoiding the need for complex, custom-coded extensions or specialized packages that often lack features like multi-biomarker handling or dynamic prediction.
• Time-Dependent Variability: Unlike previous two-stage approaches that use static summary statistics (e.g., overall standard deviation), this method treats variability as a time-varying quantity, preserving temporal information and mitigating biases associated with sparse observations.
• Scalability: The framework supports the simultaneous inclusion of multiple biomarkers, each with their own mean and variability components, a feature often limited in existing fully joint variance models.

Results

• Simulation Study: The authors evaluated the method across scenarios with linear and non-linear biomarker trajectories and varying strengths of association between variability and event risk.
  • The approach demonstrated reasonable performance with minimal bias for weak to moderate associations.
  • For stronger associations, the method exhibited some bias and reduced coverage rates, likely due to the inherent challenges of estimating time-varying variance from residuals and the fact that sample standard deviation is a biased estimator of population standard deviation.
  • Convergence was generally good, though smaller sample sizes (n=100) presented computational challenges.
• Motivating Example (Glioblastoma Trial): The method was applied to white blood cell (WBC) count data from the EORTC 1709/MIRAGE phase III trial (n=680).
  • Results indicated that both the mean level and the within-subject variability of WBC counts carry prognostic information for overall survival (OS).
  • Specifically, higher mean WBC counts were associated with a reduced hazard of death, while greater within-subject variability in WBC counts was associated with an increased hazard of death.
  • The computational efficiency was highlighted: the two-step approach took approximately 8.6 minutes, compared to over 3.6 hours for a fully joint frequentist approach (FlexVarJM) on the same data.

Significance and Claims
The paper claims that this two-step approach offers a practical and computationally efficient strategy for incorporating subject-specific biomarker variability into joint models. It enables broader application in clinical research by leveraging existing software infrastructure. The authors acknowledge that, as a two-stage method, it does not fully account for the uncertainty introduced in the first-stage estimation when fitting the second-stage survival model, which may lead to underestimated variability of association parameters compared to fully joint models. However, they posit that the trade-off is acceptable for settings where flexibility, scalability (handling multiple biomarkers), and computational feasibility are priorities. The method provides a valuable alternative for researchers seeking to explore the prognostic value of biomarker instability without the barriers of complex, custom modeling.