A Statistical Framework for Understanding Causal Effects that Vary by Treatment Initiation Time in EHR-based Studies

1. Problem Statement

Electronic Health Record (EHR) studies often span long periods (e.g., a decade or more) to accumulate sufficient sample sizes for rare treatments like bariatric surgery. Standard practice typically estimates a single, constant Average Treatment Effect (ATE) across the entire study period. However, this approach is flawed for two main reasons:

1. Evolving Treatments: In real-world settings, treatment protocols (e.g., surgical techniques for Sleeve Gastrectomy vs. Roux-en-Y Gastric Bypass) and standards of care evolve over time.
2. Covariate Shift: The population receiving treatment changes over time (e.g., shifts in BMI, age, or comorbidity distributions), which can confound the estimation of temporal trends.

Existing methods often treat calendar time as a nuisance variable to be adjusted for, rather than a source of heterogeneity to be characterized. There is a lack of statistical frameworks that can simultaneously:

• Estimate how treatment effects change over calendar time.
• Distinguish whether these changes are due to changes in treatment efficacy (the treatment itself is working differently) or covariate shift (the patient population receiving the treatment has changed).

2. Methodology

The authors propose a comprehensive statistical framework based on Sequential Target Trial Emulation and Semiparametric Efficiency Theory.

A. Study Design: Sequential Target Trials

Instead of one long observational study, the authors partition the study entry period into a sequence of $M$ "target trials" (e.g., monthly intervals).

• Eligibility ($E_m$): Patients are assessed for eligibility at the baseline of each trial $m$.
• Treatment ($A_m$): Binary exposure (Surgery vs. No Surgery).
• Outcome ($Y_m$): Continuous outcome (e.g., relative weight change).
• Covariates ($L_m$): Baseline characteristics sufficient for confounding control.

B. Estimation of Calendar Time-Specific ATEs

The goal is to estimate $\chi_m(P) = E[Y_m(1) - Y_m(0) | E_m = 1]$, the ATE for trial $m$.

• Double Robustness: The authors use influence function-based estimators that rely on pooled nuisance models (outcome regression $\mu_m$ and propensity score $\pi_m$) estimated across all trials.
• Machine Learning: To handle the complexity of EHR data and potential covariate shifts, they employ flexible machine learning methods (e.g., Random Forests) within a cross-fitting (sample splitting) framework to avoid overfitting and ensure valid inference.
• Projection to Marginal Structural Models (MSM): Since estimating a non-parametric curve for every single month is inefficient (especially with small monthly surgical volumes), the trial-specific ATEs are projected onto candidate parametric MSMs:
  $$E[Y_m(a) | E_m = 1] = \psi_0(m) + \psi(m; \beta)a$$
  A model selection procedure (minimizing a pseudorisk loss function with a penalty for complexity) is used to select the best functional form (e.g., constant, linear, cubic, or spline) that describes the trend $\psi(m; \beta)$.

C. Decomposition of Variation (The "Why")

To determine why effects vary, the authors introduce Cross-Trial Contrasts ($\chi_{j,m}$).

• Definition: $\chi_{j,m}$ represents the treatment effect in trial $m$ if the patient population had the covariate distribution of trial $j$.
• Standardization: This is achieved by reweighting the population of trial $m$ to match trial $j$ using a density ratio (transportability weight) $\xi_{j,m}$.
• Variance Decomposition: The authors define a scalar metric $\theta$ to quantify the proportion of variation in treatment effects attributable to covariate shift versus changes in underlying efficacy:
  $$\theta = \frac{1}{M} \sum_{m=1}^M \frac{\sigma_m^2}{\sigma_m^2 + \gamma_m^2}$$
  Where:
  • $\sigma_m^2$: Variance due to changing populations (covariate shift) when the treatment time is fixed.
  • $\gamma_m^2$: Variance due to changing treatment time (efficacy) when the population is fixed.
  • Interpretation:
    • $\theta \approx 0$: Variation is entirely due to covariate shift (population changes).
    • $\theta \approx 1$: Variation is entirely due to changes in treatment efficacy.
    • $\theta \in (0, 1)$: A mix of both.

D. Inference

The framework provides asymptotic normality for the estimators and uses a parametric bootstrap (approximated via asymptotic normality of the influence functions) to construct confidence intervals for $\theta$, handling boundary issues where $\theta$ might be 0 or 1.

3. Key Contributions

1. Novel Framework for Temporal Heterogeneity: The first application of semiparametric theory to explicitly model and estimate calendar time-varying treatment effects in EHR studies, moving beyond simple adjustment for time.
2. Decomposition Metric ($\theta$): Introduction of a novel summary metric to disentangle effect modification (due to population shifts) from changes in treatment efficacy. This allows researchers to answer not just if effects change, but why.
3. Projection-Based MSM Selection: A data-adaptive procedure to select the appropriate complexity of the temporal trend (e.g., constant vs. spline) while penalizing overfitting, ensuring clinically meaningful summaries.
4. Double Robustness with Cross-Fitting: Implementation of efficient, doubly robust estimators using machine learning and sample splitting, ensuring valid inference even when nuisance models are complex.

4. Results (Application to Bariatric Surgery)

The framework was applied to the DURABLE dataset (Kaiser Permanente, 2005–2011), comparing bariatric surgery (RYGB and SG) to non-surgical care.

• Temporal Trends:
  • For Surgery vs. No Surgery, a non-constant (cubic/spline) MSM was selected. The estimated 3-year weight loss effect changed from -27.4% in 2005 to -18.3% in 2011. A constant model would have masked this significant decline, reporting an average of -19.6%.
  • For RYGB vs. No Surgery, the effect also declined but stabilized after 2007.
  • For SG vs. No Surgery, a constant model was selected for 3-year outcomes, suggesting the average effect (-16.3%) was stable, though short-term (6-month) effects showed improvement.
• Decomposition ($\theta$):
  • For comparisons involving Surgery vs. No Surgery and RYGB vs. No Surgery, $\theta$ was estimated between 0.65 and 0.80. This indicates that while covariate shift played a role, a significant portion (approx. 70%) of the observed decline in efficacy was due to changes in underlying treatment efficacy (likely due to shifts in the mix of procedures, as SG became more common and RYGB less common, and potential changes in surgical technique).
  • For SG vs. No Surgery, $\theta$ was near 1.0 (0.97–0.99), indicating that observed variations were almost entirely driven by changes in the patient population (covariate shift) rather than changes in the efficacy of SG itself.
• Simulation: Extensive simulations confirmed that the method correctly identifies the functional form of the trend and accurately estimates $\theta$, even in the presence of complex covariate shifts and effect modification.

5. Significance

• Clinical Relevance: The study demonstrates that reporting a single "average" treatment effect for bariatric surgery over a decade is misleading. The efficacy of the "aggregated" surgery label has changed significantly over time, largely due to the shift from RYGB to SG and evolving patient populations.
• Methodological Advancement: This framework provides a rigorous tool for EHR researchers to move beyond static causal estimates. It allows for the characterization of the "real-world history" of a treatment, informing whether current clinical practice differs from historical data.
• Policy and Practice: By quantifying the role of covariate shift, the method helps distinguish between a treatment becoming less effective versus the treatment being applied to a different (perhaps harder-to-treat) population. This is crucial for updating clinical guidelines and interpreting historical trial data in the context of modern practice.

In summary, this paper provides a robust statistical toolkit to detect, quantify, and explain temporal heterogeneity in causal effects within large-scale observational studies, offering a more nuanced understanding of treatment effectiveness over time.