A Bivariate Copula–Driven Multi-State Model for Statistical Analysis in Medical Research

We develop and evaluate a copula-based multistate model for illness–death processes with dependent transition times. The framework couples Cox proportional hazards models for the marginal transition intensities with Archimedean copulas to capture dependence, and it is estimated via the Inference Functions for Margins (IFM) approach under right censoring. A Monte Carlo study shows that assuming independence between transitions can severely underestimate joint survival, yielding coverage as low as 40% under strong dependence, compared with 92% to 97% when copulas are used. We apply the method to a large Colombian cohort of COVID-19 patients (2021 to 2022) that includes sociodemographic, clinical, and vaccination data. The Gumbel copula best captures the strong positive dependence between hospitalization and death, producing more accurate joint survival estimates than independence-based models. Model diagnostics, including proportional hazards tests, Kaplan-Meier comparisons, hazard rate functions, and TTT plots, support the adequacy of the Cox margins. We also discuss limitations and avenues for extension, such as parametric or cure-fraction margins, nested or vine copulas, and full-likelihood estimation. Overall, the results underscore the methodological and applied value of integrating copulas into multistate models, offering a robust framework for analyzing dependent event times in epidemiology and biomedicine.