Estimating the effect of lymphovascular invasion on 2-year survival probability under endogeneity: a recursive copula-based approach

TL;DR

A semiparametric recursive copula framework is proposed that jointly specifies marginal models for both LVI, treated as an endogenous exposure, and a binary 2-year survival outcome, and links them through a flexible copula to account for latent confounding and accommodate censoring without requiring strong instruments.

Abstract

Lymphovascular invasion (LVI) is an important prognostic marker for head and neck squamous cell carcinoma (HNSC), but the true effect of LVI on survival may be distorted by endogeneity arising from unmeasured confounding. Conventional one-stage conditional models and instrument-based two-stage estimators are prone to bias under endogeneity, and sufficiently strong instruments are often unavailable in practice. To address these challenges, we propose a semiparametric recursive copula framework that jointly specifies marginal models for both LVI, treated as an endogenous exposure, and a binary 2-year survival outcome, and links them through a flexible copula to account for latent confounding and accommodate censoring without requiring strong instruments. In two simulation studies, we systematically varied sample sizes, censoring rates from 0% to 60%, and endogeneity strengths, and assessed robustness under moderate model misspecification. The proposed copula framework exhibited reduced bias and improved interval coverage compared with both one-stage and two-stage approaches while maintaining robustness to moderate misspecification. We applied the method to HNSC cases with associated clinical and microRNA data from The Cancer Genome Atlas (n = 215), and found that LVI significantly reduced 2-year survival probability by approximately 47%, with a 95% confidence interval of -0.61 to -0.29 on the probability scale. The estimated positive dependence parameter indicates that the attenuation is driven by residual dependence between unobserved components of LVI and survival. Overall, the proposed copula framework yields more credible effect estimates for survival outcomes in the absence of strong instruments, mitigating biases due to endogeneity and censoring and strengthening quantitative evidence for HNSC research.

Figures (4)

• Figure 1: SATE estimates across various censoring rates (All Coefficients Non-Zero). A) No censoring. B) 10% censoring. C) 40% censoring. D) 60% censoring. Datasets are simulated from a bivariate, Gaussian, copula, probit–probit marginal link function, assuming $\rho$=0.5, and all covariates have non-zero effects. For all simulations, the sample size is n=200 for 200 independent replicates. The x-axis denotes survival time quantiles (0.1, $...$, 0.9) based on the observed survival time. Each value corresponds to a cutoff point used to define binary survival outcomes (e.g., alive vs. dead) at that specific time threshold. Orange triangles represent the true average treatment effect. Blue circles represent the values obtained using the generalized joint regression model (GJRM). SATE, sample average treatment effect.
• Figure 2: SATE estimates across various censoring rates (Half Coefficients Zero). A) No censoring. B) 10% censoring. C) 40% censoring. D) 60% censoring. Datasets are simulated from a bivariate, Gaussian, copula, probit–probit marginal link function, assuming $\rho$=0.5, and all covariates have non-zero effects. For all simulations, the sample size is n=200 for 200 independent replicates. The x-axis denotes survival time quantiles (0.1, $...$,0.9) based on the observed survival time. Each value corresponds to a cutoff point used to define binary survival outcomes (e.g., alive vs. dead) at that specific time threshold. Orange triangles represent the true average treatment effect. Blue circles represent the values obtained using the generalized joint regression model (GJRM). SATE, sample average treatment effect.
• Figure 3: Scatter plots from simulated copula models at the same Kendall’s $\tau$$(0.5)$. The copula parameter, $\theta$, differs by family: Gaussian $\rho=0.7071$; Joe $\theta = 2.8562$; Clayton-180 $\theta=2$. See Supplementary Table S21 for the $\tau-\theta$ relationships.: (a) Gaussian copula, (b) Joe copula, and (c) $180^{\circ}$ rotated Clayton copula. The Gaussian copula exhibits symmetric dependence, whereas the Joe and rotated Clayton copulas both show pronounced upper-tail clustering but with different patterns and intensities. Axes are labeled as $q_1=\Phi^{-1}(u_1)$ and $q_2=\Phi^{-1}(u_2)$, where $(u_1,u_2)$ are samples from the respective copula families and $\Phi(\cdot)$ is the standard normal cumulative distribution function.
• Figure 4: Estimated smooth effects from the bivariate Gaussian copula model with probit–probit link functions. Panel A shows the treatment equation results, modeling the probability of LVI as a function of selected miRNA expression levels and clinical covariates. Panel B displays the outcome equation results, modeling the probability of 2-year survival. Shaded areas represent 95% confidence intervals.