A two-stage joint modeling approach for multiple longitudinal markers and time-to-event data

Joint modeling of multiple longitudinal markers and time-to-event outcomes is common in clinical studies. However, as the number of markers increases, estimation becomes computationally challenging or infeasible due to long runtimes and convergence difficulties. We propose a novel two-stage Bayesian approach for estimating joint models involving multiple longitudinal measurements and time-to-event outcomes. The proposed method is related to the standard two-stage approach, which separately estimates longitudinal submodels and then incorporates their outputs as time-dependent covariates in a survival model. Unlike the standard method, our first stage estimates separate one-marker joint models for the event and each longitudinal marker, rather than relying on mixed-effects models. From these models, predictions of expected current values and/or slopes of individual marker trajectories are obtained, thereby avoiding bias due to informative dropout. In the second stage, a proportional hazards model is fitted that includes the predicted current values and/or slopes of all markers as time-dependent covariates. To account for uncertainty in the first-stage predictions, a multiple imputation strategy is employed when estimating the survival model. This approach enables the construction of prediction models based on a large number of longitudinal markers that would otherwise be computationally intractable using conventional multi-marker joint models. The performance of the proposed method is evaluated through simulation studies and an application to the public PBC2 dataset. Additionally, it is applied to predict dementia risk using a real-world dataset with seventeen longitudinal markers. To facilitate practical use, we developed an R package, TSJM, which is freely available on GitHub: https://github.com/tbaghfalaki/TSJM.