Joint Models with Multiple Markers and Multiple Time-to-event Outcomes Using Variational Approximations

Joint models are well suited to modelling linked data from laboratories and health registers. However, there are few examples of joint models that allow for (a) multiple markers, (b) multiple survival outcomes (including terminal events, competing events, and recurrent events), (c) delayed entry and (d) scalability. We propose a full likelihood approach for joint models based on a Gaussian variational approximation to satisfy criteria (a)-(d). We provide an open-source implementation for this approach, allowing for flexible sets of models for the longitudinal markers and survival outcomes. Through simulations, we find that the lower bound for the variational approximation is close to the full likelihood. We also find that our approach and implementation are fast and scalable. We provide an application with a joint model for longitudinal measurements of dense and fatty breast tissue and time to first breast cancer diagnosis. The use of variational approximations provides a promising approach for extending current joint models.

💡 Research Summary

The paper addresses a critical gap in joint modeling of longitudinal biomarkers and time‑to‑event outcomes: the simultaneous handling of multiple markers, multiple survival endpoints (including terminal, competing, and recurrent events), delayed entry, and computational scalability. Existing software such as JMbayes, INLA, or MCMC‑based tools either lack flexibility for high‑dimensional random effects or become prohibitively slow on large datasets. To overcome these limitations, the authors propose a full‑likelihood approach that relies on a Gaussian variational approximation (GVA) to replace the intractable marginal likelihood with a tractable evidence lower bound (ELBO).

Model Specification
The longitudinal component consists of L continuous markers modeled jointly by a multivariate linear mixed‑effects model. For individual i, each marker l has a random effect vector (U_{il}\in\mathbb{R}^{r_l}) and a fixed‑effects design matrix (x_{il}(s)). The conditional mean is (\mu_{il}(s,U_{il}) = x_{il}(s)^\top\beta_l + m_{il}(s)^\top U_{il}) with residual errors (\epsilon_{ij\ell}\sim N(0,\sigma_\ell^2)). Stacking all markers yields a multivariate normal distribution with covariance matrix (\Sigma) for the residuals and (\Psi) for the random effects.

The survival component accommodates E different event types. For event e the hazard is specified as
\