A bivariate cure copula model with zero-inflated gamma frailty: dependence in both cure fractions and survival times

Abstract: In biomedical studies, paired survival data arise naturally when two event times are observed within the same subject. Existing statistical models seldom accommodate both cure fractions and complex dependence structures. In this paper, we propose a novel bivariate cure frailty-copula model for paired survival data with a cure fraction. By incorporating a zero-inflated gamma frailty, the proposed framework simultaneously accommodates a cure fraction and continuous unobserved heterogeneity among uncured subjects. Dependence between cure statuses is modeled naturally via an odds-ratio parameter, while dependence between survival times conditional on frailty is captured through a copula. We show that the proposed model includes existing bivariate cure models as special cases. Population-level rank correlation coefficients are derived for the proposed model, namely tie-corrected versions of Kendall's tau and Spearman's rho. For suitable choices of marginal distributions and copula, the joint survival function admits a closed-form expression, enabling maximum likelihood estimation and likelihood ratio testing. Simulation studies and a real data application demonstrate the practical utility of the proposed approach. An R package, curecopula, implementing the proposed methods is publicly available on GitHub.

• The paper introduces a unified bivariate cure copula model integrating zero-inflated gamma frailty to jointly model cure indicators and survival times.
• It provides closed-form likelihood solutions and robust inference, demonstrating unbiased estimation and strong hypothesis testing in simulations.
• The model’s application to diabetic retinopathy data validates its practical utility, offering clearer insights into dependencies between cure status and survival outcomes.

Bivariate Cure Copula Model with Zero-Inflated Gamma Frailty: Dependence in Both Cure Fractions and Survival Times

Introduction

The paper presents a comprehensive statistical framework for modeling paired survival data with cure fractions, introducing a bivariate cure frailty-copula model with a zero-inflated gamma frailty. Traditional bivariate survival models typically focus on capturing dependence either through discrete frailty or copula-based techniques, and often treat cure rates as secondary effects or as marginal distributions without fully characterizing their joint behavior. This work systematically integrates cure fractions, continuous unobserved heterogeneity (through frailty), and flexible copula-driven dependence, allowing for explicit and independent parameterization of each component. Crucially, the model facilitates dependence between both the cure indicators and the survival times among uncured subjects, achieved through an odds-ratio parameter for cure status and a copula for survival time dependence.

Model Formulation and Properties

The adopted model considers two possibly correlated survival outcomes $(T_1, T_2)$ per subject, indexed either by individual organs or paired individuals. Each marginal survival process incorporates a latent binary cure indicator $(X_1, X_2)$, parameterized by population cure fractions $p_1$ and $p_2$. Continuous heterogeneity among uncured subjects is handled via a shared gamma frailty variable $W$. The model employs a zero-inflated gamma frailty, such that cured margins receive zero hazard, while uncured margins are subject to positive frailty and continuous hazard.

The joint survival function is decomposed according to cure configurations, with closed-form solutions attainable for specific copulas and Weibull marginals. The model generalizes prior approaches, reducing to notable special cases:

• The bivariate cure-copula model as in Chatterjee & Shih (2001) when frailty vanishes.
• The correlated gamma frailty model of Wienke et al. (2003) under restricted copulas.
• The frailty-copula model for multivariate survival without cure fractions when cure parameters approach zero.

Parametric flexibility is preserved by decoupling copula and frailty specifications, allowing independence or dependence (Gumbel, FGM, BB1 copulas) to be imposed in modeling residual correlation among uncured subjects. Dependence between the cure indicators is rigorously introduced via an odds ratio $R$, offering separate identification of cure and susceptible-time parameters.

Population Rank Correlation Coefficients

A technical advance is the derivation of tie-adjusted rank measures for bivariate survival with cure, namely Kendall’s $\tau_b$ and Spearman’s $\rho_b$. These coefficients account for mass at infinity (cured margins) and provide interpretable dependence measures at both population and subpopulation levels. Closed-form expressions are given for the independence (Clayton) and Gumbel (BB1) copulas, establishing their attainable range and asymptotic regimes. The model admits perfect concordance and negative association as limits, and the dependence structure is fully characterized by $R$, frailty variance $\gamma$, and copula parameters.

Likelihood-Based Estimation and Testing

The paper develops likelihood-based inference, including closed-form log-likelihoods for censored bivariate data and explicit derivations for maximum likelihood estimation under different copula choices. Wald-type confidence intervals are constructed for all positive parameters using suitable transformations, and the odds ratio $(X_1, X_2)$0 is handled carefully, allowing for region-specific estimation and testing.

A likelihood ratio test is proposed for independence in cure status ($(X_1, X_2)$1), with asymptotic $(X_1, X_2)$2 calibration, and its finite-sample behavior is evaluated via simulation.

Simulation Studies

Extensive simulation studies assess finite-sample behavior for both parameter estimation and hypothesis testing. For MLE performance, the estimator demonstrates unbiasedness and strong coverage across a range of sample sizes and model specifications (e.g., positive vs. negative cure dependence, high or low frailty). Increasing sample size improves bias and MSE as expected by asymptotic theory.

For the LRT, empirical power curves exhibit desirable properties: Type I error is tightly controlled near nominal levels, and power increases smoothly as the odds ratio departs from unity, with asymmetric sensitivity depending on direction of association. Censoring rates and covariate effects are chosen to mirror realistic biomedical scenarios and are incorporated seamlessly into the regression structure for cure fractions and survival times.

Figure 1: Empirical power of the likelihood ratio test for $(X_1, X_2)$3 across parameter settings, displaying controlled Type I error and efficient detection as $(X_1, X_2)$4 deviates from 1.

Real Data Application

The proposed approach is applied to paired eye survival data from the Diabetic Retinopathy Study, modeling treated and untreated eyes per patient. Multiple copula and odds ratio regimes are compared via AIC and BIC, showing that the independence copula with covariates and $(X_1, X_2)$5 achieves optimal fit. The frailty parameter yields moderate within-susceptible dependence, while the cure odds ratio is not statistically distinguishable from unity, indicating that the cure process for the two margins operates independently in this cohort.

Kaplan–Meier and model-based marginal survival functions are overlaid for diagnostic purposes. The proposed estimates for long-term cure fractions ($(X_1, X_2)$6 for treated, $(X_1, X_2)$7 for untreated) correspond well with nonparametric bounds, and the fitted model tracks observed survival dynamics accurately.

Figure 2: Kaplan–Meier curves (solid), fitted model-based survival (dashed), and estimated cure fractions (dot-dashed) for both treated and untreated eyes.

Practical and Theoretical Implications

The model’s generality enables direct handling of cure dependence and susceptible-time association in diverse biomedical contexts, including bilateral and familial disease, clustered failure mechanisms, and complex survival datasets with covariates. The zero-inflated gamma frailty yields better fit and flexibility than discrete frailty or mixed Poisson approaches, and separating cure and susceptible-time parameters enhances interpretability and clinical inference.

From a theoretical standpoint, the derivation of tie-adjusted rank measures, explicit closed-form solutions, and likelihood-based inference extends the capacity for detailed dependence analysis in bivariate cure models. The availability of an R package implementing these tools positions this framework for immediate practical adoption.

Conclusion

The bivariate cure frailty-copula model with zero-inflated gamma frailty described in this paper provides a unified, rigorous method for paired survival analysis with cure fractions, integrating continuous heterogeneity and flexible dependence structures. Closed-form likelihoods, robust estimation, and explicit hypothesis testing are achieved, overcoming prior limitations seen in discrete frailty or coupled cure-frailty approaches. The methodology is validated through simulation and real data analyses, with practical fit and inferential advantages demonstrated. Future developments include extending to broader frailty distributions, enhancing tail dependence analysis, and further clinico-statistical applications.