Bayesian Joint Modeling of Multivariate Longitudinal and Time-to-Event Data Using Multivariate Skew-Normal Random Errors

Joint models have become the standard framework for simultaneously analyzing longitudinal biomarkers and time-to-event outcomes because they account for their underlying association. However, most existing joint models are restricted to a single longitudinal outcome and assume multivariate normal random errors. These assumptions are often unrealistic in chronic disease studies where multiple correlated biomarkers are measured repeatedly and exhibit asymmetric distributions. Ignoring these features may lead to biased estimation and reduced efficiency. A Bayesian multivariate joint model is proposed that simultaneously accommodates multiple correlated longitudinal outcomes and survival data by incorporating a skewed multivariate normal distribution for random errors. The methodology is motivated by chronic disease fasting blood sugar (FBS), systolic blood pressure (SBP), and time to chronic complications among individuals with type 2 diabetes (T2D) and hypertension, and its performance is assessed via a simulation study. Bayesian inference is conducted using Markov chain Monte Carlo (MCMC) algorithms implemented in Just Another Gibbs Sampler (JAGS). The limited simulation comparison suggests differences in estimation performance between the models but does not establish uniform superiority. In the application dataset, the proposed model captured dependence and asymmetry in the longitudinal biomarkers. The estimated association between the SBP random-slope component and the hazard of chronic complications was positive, although the 95% credible interval included zero, indicating substantial uncertainty regarding the direction and magnitude of the association.