Cox Regression on the Plane

Technical Summary: Cox Regression on the Plane

Problem Statement
The Cox proportional hazards model is the standard for univariate survival analysis, yet extensions to bivariate survival data (correlated failure times) remain limited. Existing approaches face significant trade-offs:

• Shared frailty and copula models: These often rely on parametric assumptions regarding the dependence structure (e.g., specific copula families or frailty distributions). If the dependence structure is misspecified, covariate effects at the joint survival level may be uninterpretable or incorrect.
• Multiplicative hazard extensions: Previous attempts to specify three separate proportional hazards models (for single-failure and double-failure rates) generally fail to yield a valid bivariate survival distribution because the resulting hazard functions do not satisfy necessary structural constraints.

The authors identify a gap for a method that preserves the semiparametric nature of the Cox model, yields a valid bivariate survival distribution, and offers transparent covariate interpretation without requiring the specification of a copula or frailty distribution.

Methodology
The paper proposes two semiparametric extensions based on a generalized Lehmann representation of the survival function.

1. The Models:

   • Simple Lehmann Model: A direct extension where the bivariate survival function is $S(t_1, t_2 | Z) = [S_0(t_1, t_2)]^{\exp(b^T Z)}$. This imposes a common proportional effect on the log-survival of both marginals and the joint dependence, preserving the baseline dependence structure.
   • Generalized Lehmann Model: A more flexible formulation that allows distinct regression parameters for the two marginal hazards ($\alpha, \beta$) and the dependence component ($\gamma$). It is defined as:
     $$S(t_1, t_2 | Z) = [S_{0T1}(t_1)]^{\exp(\alpha^T Z)} [S_{0T2}(t_2)]^{\exp(\beta^T Z)} \left[ \frac{S_0(t_1, t_2)}{S_{0T1}(t_1)S_{0T2}(t_2)} \right]^{\exp(\gamma^T Z)}$$
     This model treats the dependence component as a multiplicative factor on the baseline dependence ratio. The authors derive the corresponding hazard functions and establish sufficient conditions (e.g., constraints on $\gamma$ relative to $\alpha$ and $\beta$) to ensure the resulting functions define valid, non-negative survival distributions.

2. Estimation:

   • The authors utilize a pseudo-observations approach to avoid the complexities of partial likelihood in the bivariate setting, where baseline hazards do not cancel out.
   • Simple Model: Estimated via Generalized Estimating Equations (GEE) using bivariate pseudo-observations as responses in a GLM with a complementary log-log (cloglog) link.
   • Generalized Model: Estimated via a two-step procedure:
     1. Estimate marginal parameters ($\alpha, \beta$) using univariate pseudo-observations and standard GEE.
     2. Construct modified responses for the dependence component (the ratio of joint survival to the product of marginals) and estimate $\gamma$ using a second GEE step.
   • The paper establishes the consistency and asymptotic normality of these estimators under independent censoring.

Key Results

• Simulation Studies: The authors evaluated the methods against data generated from various copulas (Frank, Clayton, Gumbel-Barnett) representing both positive (PQD) and negative (NQD) dependence.
  • The estimators demonstrated low bias and coverage close to nominal levels for sample sizes of $n=800$.
  • The proposed models maintained competitive predictive accuracy (Mean Absolute Error) compared to flexible copula-based regression frameworks (Marra and Radice, 2020), even when the simple Lehmann model was misspecified (i.e., when $\alpha \neq \beta \neq \gamma$).
  • The methods proved robust to model misspecification when applied to data generated from a log-normal distribution, where neither the Lehmann nor copula models were strictly correct.
• Data Application (Global Retinoblastoma Outcome Study - GROS):
  • The method was applied to 977 bilateral retinoblastoma cases to analyze times to enucleation for both eyes.
  • Dependence Structure: Nonparametric estimates indicated negative quadrant dependence (NQD) between the eyes. However, a standard copula-based approach selected a Clayton copula, which implies positive dependence (PQD), leading to conflicting conditional survival probabilities.
  • Interpretability: The Lehmann models correctly captured the NQD structure without forcing a parametric copula. The results showed that if one eye was enucleated early, the probability of the other eye surviving longer was higher than if the first eye survived, a finding consistent with clinical expectations but missed by the copula model.
  • Covariate Effects: The models identified that older age at presentation, higher tumor stage, and income level significantly affected survival probabilities, while sex and family history did not.

Significance and Claims
The paper claims to provide a transparent, fully semiparametric framework for assessing covariate effects on both marginal survival probabilities and their dependence. Its primary contributions are:

1. Avoidance of Copula Misspecification: By modeling covariate effects directly on the bivariate survival function rather than through a parametric copula family, the approach avoids the risk of incorrect dependence structure assumptions, which is critical for interpreting conditional probabilities.
2. Validity and Interpretability: The proposed models are proven to define valid bivariate survival distributions under specific conditions, and the regression coefficients have direct interpretations on the survival scale (e.g., $\gamma$ governs how covariates amplify or attenuate baseline dependence).
3. Novel Estimation: The authors present the first two-step pseudo-observations-based method for bivariate survival regression, establishing its theoretical properties.
4. Clinical Relevance: The application to retinoblastoma data demonstrates that standard model selection criteria (like BIC) can favor copulas with incorrect dependence structures (PQD vs. NQD), whereas the proposed framework yields clinically interpretable results that align with empirical evidence and domain knowledge.

The authors note limitations, including the computational cost of recomputing survival estimators for pseudo-observations and the assumption of independent censoring, but emphasize that the approach offers a robust alternative for bivariate survival analysis where dependence structure is unknown or complex.