Technical Summary: Debiased Causal Mediation Analysis in Ultra-High-Dimensional Settings in the Presence of Interaction Effects

1. Problem Setting and Motivation

This paper addresses the challenge of estimating natural direct effects (NDE) and natural indirect effects (NIE) in causal mediation analysis when both the pre-treatment covariates ($X \in \mathbb{R}^p$) and the post-treatment mediators ($M \in \mathbb{R}^q$) are ultra-high-dimensional (UHD), meaning $p \wedge q \gg n$. Crucially, the proposed framework accommodates interaction effects between the binary treatment ($A$) and both the covariates and the mediators.

Existing literature on high-dimensional mediation analysis has largely relied on heuristic screening, dimension reduction, or pathway selection, often lacking rigorous inferential theory for natural effects in the presence of UHD variables. Furthermore, recent theoretical advances in high-dimensional mediation (e.g., Guo et al., 2022–2024; Lin et al., 2025) typically assume fixed-dimensional covariates or exclude treatment-covariate/mediator interactions. The specific difficulty addressed here is the estimation of the mediation functional $\theta_{10} = E[Y(1, M(0))]$, which involves a "cross-arm" interaction: it combines parameters learned from the treatment arm (the outcome regression coefficients $\gamma_1$) with those learned from the control arm (the mediator regression matrix $B_0$). This structure prevents the direct application of standard debiasing techniques for linear or bilinear functionals, as the target involves a composite term $\bar{X}^\top B_0^\top \gamma_1$ where $B_0 \bar{X}$ may be dense and not consistently estimable in the $\ell_1/\ell_2$ norm.

2. Methodology

The authors propose a novel multi-step debiasing estimator that combines high-dimensional regression with carefully designed balancing steps. The method utilizes a cross-fitting procedure to split the data into two folds ($D_1, D_2$) to avoid overfitting and ensure valid inference.

The Estimation Procedure (Algorithm 1)

The core of the methodology estimates the parameter $\theta_{10} = \bar{X}^\top (\beta_1 + B_0^\top \gamma_1)$ through the following steps:

1. Balancing for Approximation (Step 1): Using control observations ($A=0$), the method constructs a balancing weight vector $\hat{\tau}_2$ by solving an optimization problem. This weight approximates the mean covariate vector $\bar{X}$ in the space of the control covariates. Crucially, the product $\hat{m}_c = M_c^\top \hat{\tau}_2$ serves as a direct approximation of the unknown vector $B_0 \bar{X}$ in the $\ell_\infty$ norm, bypassing the need to estimate the full UHD matrix $B_0$.
2. Outcome Regression (Step 2): Using treatment observations ($A=1$), a Lasso regression estimates the outcome coefficients $\hat{\phi}_1 = (\hat{\beta}_1, \hat{\gamma}_1)$.
3. First-Stage Debiasing (Step 3): The estimator $\hat{\phi}_1$ is debiased in the direction $a = (\bar{X}^\top, \hat{m}_c^\top)^\top$. This step corrects for the regularization bias of the Lasso but introduces a new bias because the direction $\hat{m}_c$ is an approximation of $B_0 \bar{X}$, not the true value.
4. Second-Stage Debiasing (Step 4): To remove the bias introduced by the approximation in Step 3, a second regression is performed on the control data. The pseudo-outcome $M_c^\top \hat{\gamma}_1$ is regressed on $X_c$ to estimate $b = B_0^\top \hat{\gamma}_1$. This estimator is then debiased using the same balancing weight $\hat{\tau}_2$ used in Step 1.
5. Combination (Step 5): The final estimator is constructed by combining the debiased outcome term and the debiased approximation term, subtracting a correction term to account for the overlap.
6. Cross-Fitting (Step 6): The roles of $D_1$ and $D_2$ are swapped, and the results are averaged to obtain the final cross-fitted estimator $\hat{\theta}_{10}^{cf}$.

The estimators for NDE and NIE are constructed by combining $\hat{\theta}_{10}^{cf}$ with analogous estimators for $\theta_{11}$ and $\theta_{00}$ (which are simpler as they do not involve cross-arm interactions).

3. Theoretical Contributions and Results

The paper establishes that the proposed estimator is $\sqrt{n}$-consistent and asymptotically normal under fairly general sparsity, regularity, and restricted eigenvalue conditions.

• Bias Structure: The authors prove that the estimation error decomposes into a leading stochastic term and a remainder term $\Delta_n$. The remainder term exhibits a product form (e.g., $\|\text{approximation error}\|_\infty \times \|\text{estimation error}\|_1$). This structure ensures that the bias is asymptotically negligible ($o_p(n^{-1/2})$) even when individual estimation errors are large, provided the sparsity conditions $(1+s)(1+k)\log(p+q)/\sqrt{n} \to 0$ hold.
• Asymptotic Normality: The paper derives the asymptotic variance and constructs consistent estimators for it, enabling the construction of valid confidence intervals for NDE, NIE, and the total effect. The results hold for both sample-average and population-average targets, and are extended to handle conditionally heteroskedastic errors.
• Comparison to Existing Methods: The authors demonstrate that their approach avoids the strong "beta-min" (minimum signal strength) and support recovery assumptions required by some existing high-dimensional inference methods (e.g., Guo et al., 2023). Their method allows for weaker signals and higher sparsity levels, provided the product of sparsity and dimension grows slowly enough relative to the sample size.

4. Empirical Evaluation

The methodology was evaluated through extensive simulation studies and a real-data application.

• Simulations: The proposed debiased estimator ($\hat{\theta}_{10}^{cf}$) was compared against a standard plug-in estimator and a modified version of the method by Guo et al. (2022).
  • The debiased estimator consistently achieved lower Root Mean Squared Error (RMSE) and coverage probabilities close to the nominal level (95%) across various regimes, including increasing sparsity, vanishing minimum signal strength, and heteroskedastic noise.
  • In contrast, the standard estimator exhibited substantial bias, and the modified Guo et al. method showed instability and poor coverage in ultra-high-dimensional settings, particularly when sparsity was high or signals were weak.
• Real Data Application: The method was applied to the TCGA lung cancer dataset to estimate the effect of smoking on survival time, mediated by DNA methylation.
  • The analysis involved $n=738$ patients with $p=42$ covariates and $q=264,797$ mediators.
  • Results indicated that the total effect of smoking on survival time was primarily mediated through DNA methylation (NIE), with the indirect effect being statistically significant while the direct effect was not.
  • A comparison with a fixed-dimensional OLS approach (using only 6 pre-selected CpG sites) yielded contradictory and counter-intuitive results (suggesting smoking increased survival), highlighting the importance of high-dimensional methodology in complex biological systems.

5. Significance and Claims

The paper claims to be the first work to provide a methodology for estimating natural direct and indirect effects in the presence of simultaneously ultra-high-dimensional covariates, ultra-high-dimensional mediators, and their interactions with the treatment variable.

Key claims regarding the paper's significance include:

1. Structural Innovation: The proposed multi-step debiasing strategy addresses a unique structural challenge (cross-arm nuisance parameters involving UHD matrices and vectors) that existing debiasing techniques for linear or bilinear functionals cannot handle directly.
2. Weaker Assumptions: The theoretical framework operates under weaker signal-strength assumptions than competing high-dimensional mediation inference methods, allowing for the detection of effects even when some active coefficients are weak.
3. Practical Utility: The method provides valid inference for complex biological mechanisms (like epigenetic mediation) where the number of potential mediators vastly exceeds the sample size, a scenario where traditional low-dimensional selection methods may fail or yield biased results.

The authors conclude that while the current work focuses on linear models, the developed framework offers a foundation for future extensions to generalized linear models and other structural equation models with similar cross-arm complexities.