Estimating equations for causal survival analysis with pooled logistic regression

Technical Summary: Estimating Equations for Causal Survival Analysis with Pooled Logistic Regression

Problem Statement
Pooled logistic regression is a widely used technique in causal survival analysis, particularly as a nuisance model for g-computation and inverse probability weighting (IPW) estimators. The standard implementation requires transforming data into a "long" format, where each row represents a person-period (an individual for a specific time interval). This approach faces three primary computational and practical challenges:

1. Data Volume: The long dataset can become unwieldy with large sample sizes or fine-grained time intervals, often necessitating the coarsening of time (e.g., using weeks instead of days), which results in information loss.
2. Computational Burden: For causal inference (e.g., estimating marginal risks), variance estimation typically relies on the nonparametric bootstrap. Because the pooled logistic model must be re-fitted hundreds or thousands of times, this process is computationally expensive, often taking hours to complete, which threatens reproducibility.
3. Implementation Errors: The manual transformation of data into long formats introduces a risk of coding errors (e.g., misaligned time or miscoded events).

Methodology
The authors propose reframing pooled logistic regression using estimating equations (specifically M-estimation) to bypass the need for a long dataset and eliminate the reliance on bootstrapping for variance estimation.

• Estimating Equations Formulation:

  • Instead of maximizing a likelihood over a long dataset, the method defines a unit-specific estimating function $\psi(O_i; \theta)$ that sums contributions across all time intervals for a single individual.
  • The score equation for the regression parameters $\beta$ is expressed as the sum of residuals weighted by covariates over the risk set, aggregated per unit rather than per person-period.
  • To estimate the marginal risk (the target parameter), a second estimating function is constructed that transforms the predicted discrete-time hazards into cumulative risks.
  • These equations are stacked to jointly estimate nuisance parameters and the target causal effect.

• Variance Estimation:

  • The method utilizes the empirical sandwich variance estimator to consistently estimate the variance of the target parameter. This accounts for the uncertainty in the nuisance model parameters without resampling.
  • This allows for the construction of Wald-type confidence intervals (CIs) and, for time-varying effects, confidence bands (CBs) using the sup-t method.

• Vectorization and Implementation:

  • The authors describe a vectorized implementation of the estimating functions using matrix operations (Hadamard products and cumulative products) to avoid slow iterative loops common in interpreted languages like Python and R.
  • This approach reduces memory requirements compared to the standard long-data approach, particularly when the number of time intervals ($K$) or sample size ($n$) is large.
  • The implementation is available in the open-source Python library delicatessen and the R package geex.

Key Results
The proposed method was evaluated through two applied examples using public data and a Monte Carlo simulation study.

1. Applied Examples:

   • Bladder Cancer Trial: The estimating equation approach reduced run-times from over a minute (standard implementation with disjoint indicators) or several hours (bootstrap) to less than one second, while producing identical risk difference estimates and confidence intervals.
   • HIV Cohort Study (WIHS): In a scenario with 1,164 subjects and 3,653 days of follow-up, the standard implementation with disjoint indicators failed due to memory constraints (requiring ~4.2 million rows). The bootstrap could not be run in parallel due to memory limits. The estimating equation approach successfully computed results in under two minutes (days scale) and seconds (months scale) without data coarsening.

2. Simulation Study:

   • A simulation with 5,000 iterations across varying sample sizes ($n=250$ to $2000$) assessed the performance of the empirical sandwich variance estimator.
   • Coverage: When the functional form for time was correctly specified (e.g., log-linear, splines, or disjoint indicators), the estimator achieved nominal 95% confidence interval coverage.
   • Bias: Misspecification of the time functional form (e.g., using only an intercept) resulted in bias and poor coverage, highlighting the importance of flexible time modeling.
   • Variance Accuracy: The standard error ratio (ASE/ESE) was near 1.0 for correctly specified models, indicating the sandwich estimator accurately captures the variability.

Significance and Claims
The paper claims that this implementation offers a computationally efficient alternative to the standard pooled logistic regression approach for causal survival analysis. Its primary contributions are:

• Elimination of the Long Dataset: It avoids the creation and storage of massive person-period datasets, thereby removing the need to coarsen time intervals and preserving data granularity.
• Efficient Inference: By replacing the nonparametric bootstrap with the empirical sandwich variance estimator, it drastically reduces computation time (from hours to seconds/minutes) while maintaining valid statistical inference.
• Feasibility: It enables the analysis of datasets and time scales that are computationally intractable under the standard bootstrap approach.

The authors note that while the focus is on pooled logistic regression, the framework is generalizable to other pooled regression models (e.g., complementary log-log or multinomial logistic for competing risks). They acknowledge limitations, including the current inability to handle time-varying treatments directly within this specific formulation and the reliance on discretizing time, though they note that theoretical extensions for growing dimensions exist.