Data-Adaptive Identification of Effect Modifiers through Stochastic Shift Interventions and Cross-Validated Targeted Learning

To estimate the causal effect of stochastic interventions effectively, we aim to construct a semiparametric efficient estimator. Such an estimator achieves the smallest possible asymptotic variance among regular asymptotically linear (RAL) estimators, satisfying:

where $\phi(O)$ is the efficient influence function, and $\mathbb{E}_{n}$ denotes the empirical mean. By the Central Limit Theorem, RAL estimators achieve the asymptotic variance bound $\mathbb{E}[\phi^{2}(O)]$, facilitating valid statistical inference.

In the context of our stochastic shift intervention, the efficient influence function $\phi_{P}(O)$ of $\Psi$ at distribution $P$ is given by:

Several estimators can be employed for estimating our stochastic shift parameter, including inverse probability weighting (IPW), augmented inverse probability weighting (AIPW), and targeted maximum likelihood estimation (TMLE). However, IPW estimators can be inefficient and sensitive to model misspecification, and AIPW estimators may not respect the bounds of the outcome, potentially resulting in greater variability. Therefore, we choose TMLE to estimate our parameter due to its efficiency and robustness in finite samples Li et al. (2022).

Targeted learning begins with initial estimates of the outcome regression $\hat{\overline{Q}}(A,W)$ and the exposure density $\hat{g}(A\mid W)$, which can be obtained through flexible machine learning techniques. These initial estimates are then updated using a clever covariate based on the efficient influence function. The refinement process involves fitting a parametric fluctuation model that targets the parameter of interest. This results in a final estimate $\hat{\overline{Q}}^{\star}(A-\delta,W)$ that is efficient and less biased. The updated estimate is then used to compute the TMLE of $\Psi$, and the influence function is utilized to derive confidence intervals.

The oracle parameter, denoted as $\psi^{*}$, represents the maximal differential impact of an intervention in a specific subregion of the covariate space. It is defined as follows:

where:

• •

  $\boldsymbol{A}=(A_{1},A_{2},\ldots,A_{p})$ denotes the vector of exposures in the mixture, with $A_{i}$ representing the $i$-th exposure and $\mathbf{A}_{-i}$ denoting all exposures excluding $A_{i}$.

• •

  $\mathcal{V}$ and $\mathcal{V}^{c}$ denote a subregion of the covariate space $\mathcal{W}$ and its complement, respectively.

• •

  $Y_{A_{i}-\delta_{i},\mathbf{A}_{-i}}$ represents the potential outcome under a stochastic shift intervention where exposure $A_{i}$ is shifted by $\delta_{i}$, while other exposures $\mathbf{A}_{-i}$ remain unchanged.

$\psi^{*}$ identifies the combination of exposure $A_{i}$ and subregion $\mathcal{V}$ that maximizes the difference in expected outcomes between $\mathcal{V}$ and its complement $\mathcal{V}^{c}$. This parameter captures effect modification, highlighting where the intervention has the most substantial and relevant impact.

In general, we are interested in the causal exposure-covariate region combination where the effect of the intervention, compared to the complementary covariate space, is maximally different, controlling for other exposures. From now on, we refer to these regions as $\mathcal{V}$ and $\mathcal{V}^{c}$.

This method allows us to identify and focus on subpopulations where the intervention has the most substantial and relevant impact, thereby informing targeted public health strategies and interventions.

This section reformulates the causal parameter of interest in terms of observable quantities, thereby defining a statistical parameter that can be estimated from the data. We consider a fixed subpopulation $\mathcal{V}\subseteq\mathcal{W}$, but emphasize that subpopulation-specific parameters are interpreted as features of the full population distribution. Let $\pi_{\mathcal{V}}=P(W\in\mathcal{V})$ denote the probability that a randomly selected unit from the full population belongs to $\mathcal{V}$.

Under standard identification assumptions, including no unmeasured confounding and positivity, the causal parameter $\psi^{*}$ is identified as:

where $\overline{Q}_{0}(a_{i},\mathbf{a}_{-i},w)=\mathbb{E}[Y\mid A_{i}=a_{i},\mathbf{A}_{-i}=\mathbf{a}_{-i},W=w]$.

Fixing $\mathcal{V}$, we define the statistical parameter $\psi_{\mathcal{V}}$ as:

where

This formulation ensures that $\Psi_{\mathcal{V}}(P)$ and $\Psi_{\mathcal{V}^{c}}(P)$ are understood as expectations taken over the full population distribution, normalized by $\pi_{\mathcal{V}}$ and $1-\pi_{\mathcal{V}}$ respectively. The parameter $\psi_{\mathcal{V}}$ thus represents a feature of the entire population, focusing attention on differences between a particular subpopulation $\mathcal{V}$ and its complement.

Building upon the statistical parameters and influence functions defined in the previous section, we now describe a procedure for identifying subpopulations with maximal effect modification for a given exposure. The estimation involves a two-stage approach:

First Stage: We obtain efficient estimates of the Individual Intervention Effects (IIEs) and the corresponding Influence Curve Estimates (ICEs) for each exposure $A_{k}$ in the mixture $\boldsymbol{A}=(A_{1},A_{2},\ldots,A_{p})$.

Second Stage: Utilizing these estimates, we apply a greedy partitioning algorithm to find covariate regions $\mathcal{V}$ where the effect modification $\psi_{\mathcal{V}}$ is maximized and statistically significant.

As discussed, sample splitting is required to avoid bias incurred by discovering subregions of the covariate space and estimating the effect of subpopulation intervention using the same data. To achieve this, we employ $K$-fold cross-validation, where the sample data are partitioned into $K$ folds of approximately equal size. For example, consider $K=10$. In each iteration $k=1,\ldots,K$:

1. 1.

   Training and Validation Split:

   • •

     The training set $T_{n,k}\subset\{1,\ldots,n\}$ consists of the indices from all folds except the $k$-th fold, comprising $\frac{K-1}{K}$ of the data (e.g., 90% for $K=10$).

   • •

     The validation set $V_{n,k}\subset\{1,\ldots,n\}$ consists of the indices in the $k$-th fold, comprising $\frac{1}{K}$ of the data (e.g., 10% for $K=10$).

2. 2.

   Stage 1: Subpopulation Identification and Parameter Estimation on Training Set:

   • •

     Use the training data $T_{n,k}$ to identify the maximal exposure-effect modifier sets within the mixture using the recursive partitioning algorithm described in Section 3.3.

   • •

     Estimate the nuisance parameters $\hat{\overline{Q}}$ and $\hat{g}$ using appropriate machine learning methods.

   • •

     Compute the individual intervention effects (IIE) and influence curve estimates (ICE) for each observation in $T_{n,k}$.

3. 3.

   Stage 2: Parameter Estimation on Validation Set:

   • •

     Apply the subpopulation rule $\mathcal{V}_{k}$ identified from the training set to the validation data $V_{n,k}$, assigning observations to $\mathcal{V}_{k}$ or $\mathcal{V}_{k}^{c}$.

   • •

     Estimate the subpopulation stochastic shift intervention parameter $\Psi_{\mathcal{V}_{k}}(P)$ using the CV-TMLE approach with the nuisance parameters estimated from $T_{n,k}$.

After completing all $K$ folds, we aggregate the estimates across folds to obtain the final estimate of the effect modification parameter:

If the identified subpopulation $\mathcal{V}_{k}$ is consistent across folds, this pooled estimate leverages the full data for increased efficiency.

Upon completion of the cross-estimation procedure, a pooled TMLE update provides a summary measure of the oracle maximum effect modification parameter across $K$-folds. In each fold $k$, we obtain estimates for $\overline{Q}$ and $g$ for the respective regions identified as the maximum effect-modifying regions and their complements. These regions may vary across different folds.

For the pooled TMLE update, we perform the following steps:

1. 1.

   Stacking Nuisance Parameters: Stack the estimates of the nuisance parameters $\hat{\overline{Q}}$ and $\hat{g}$ from each training sample $T_{n,k}$.

2. 2.

   Pooled TMLE Update: Apply a pooled TMLE update on the cumulative initial estimates. This involves:

   • •

     Combining the initial outcome regression estimates and propensity score estimates across all folds.

   • •

     Performing a TMLE update on this pooled data to refine the estimates, ensuring they are both efficient and unbiased.

3. 3.

   Aggregating Estimates: Compute the pooled estimate of the effect modification parameter by averaging the estimates across all folds:

   $$\Psi_{\text{pooled}}=\frac{1}{K}\sum_{k=1}^{K}\Psi_{F_{T_{n,k}}}\left(V_{n,k}\right),$$

   (15)

   where $\Psi_{F_{T_{n,k}}}$ identifies the optimal exposure-covariate modifier set and formulates a plug-in estimator using the training data $T_{n,k}$. $V_{n,k}$ refers to the validation sample for deriving the estimates.

4. 4.

   Interpretation of Pooled Estimates: The pooled estimate $\Psi_{\text{pooled}}$ represents an average across folds. If the identified subpopulations $\mathcal{V}_{k}$ and $\mathcal{V}_{k}^{c}$ are consistent across all folds (e.g., consistently identifying age $<22$ for $\mathcal{V}$ and age $\geq 22$ for $\mathcal{V}^{c}$), the pooled estimates are interpretable and provide increased statistical power. This leverages the full data, resulting in tighter confidence intervals and more reliable inference.

When analyzing a mixture of multiple exposures, it is often important to assess the relative importance of each exposure component. Our approach facilitates this by estimating a stochastic shift for each exposure and deriving the corresponding influence curve estimates (ICE) to identify regions with maximal effect modification.

Since our target parameter requires stochastic shift interventions for each exposure as a preliminary step before recursive partitioning, we can efficiently leverage these estimates to obtain CV-TMLE for each exposure. The procedure is as follows:

Utilizing Existing Estimates: In the first stage of our methodology, we estimate stochastic shift interventions for each exposure $A_{i}$ in the mixture $\boldsymbol{A}$ and derive the corresponding ICEs. These estimates are obtained through the training phase of our cross-validation process and are applied to the validation folds.

Pooled CV-TMLE Estimates: After estimating the stochastic shifts and ICEs for each exposure across all folds, we aggregate these results to compute CV-TMLE estimates for each exposure shift. This pooling process involves:

• •

  Saving Validation Estimates: We store the TMLE-updated estimates of the outcome changes due to shifts in each exposure from the validation folds.

• •

  Aggregating Estimates: By averaging these estimates across all folds, we obtain pooled CV-TMLE estimates for each exposure $A_{i}$. This aggregation enhances the reliability and statistical power of our estimates.

Variable Importance Measures: The pooled CV-TMLE estimates for each exposure shift allow us to assess the relative importance of each exposure within the mixture. Specifically, these estimates quantify the impact of shifting each exposure $\delta_{i}$ on the outcome, providing a clear measure of variable importance. This comprehensive estimation framework not only identifies regions with maximal effect modification but also offers insights into the marginal effects of each exposure, thereby supporting informed and targeted public health interventions.

In summary, because our target parameter inherently requires stochastic shift interventions for each exposure as part of the initial estimation process, we can directly utilize these results to obtain CV-TMLE estimates. This approach simplifies the workflow and enables the derivation of variable importance measures for the entire exposure mixture.

Accurate estimation of the conditional exposure densities $g(A\mid W)$ is crucial for evaluating the effects of stochastic interventions within a causal inference framework. We employ two approaches for density estimation: direct conditional density estimation and a reparameterization strategy that transforms density ratio estimation into a classification problem. Both methods are implemented using the Super Learner ensemble van der Laan et al. (2007), which optimally combines multiple algorithms to produce robust and flexible estimates.

To evaluate the performance of our proposed method, we conducted a simulation study using synthetic data that mirrors a realistic environmental epidemiology scenario. We generated a large population with 10000 observations to estimate ground truth of the subpopulation intervention effect. Three binary confounders, $W_{1}$, $W_{2}$, and $W_{3}$, were independently generated from Bernoulli distributions with a success probability of 0.5.

The exposures ($A_{1},A_{2},A_{3}$) were generated as follows:

The outcome $Y$ was then simulated based on the exposures and covariates, incorporating effect modification effects:

where $\epsilon\sim\mathcal{N}(0,1)$ represents random noise.

To make detection of the modifying region more difficult and arguably more realistic compared to the binary case, we also create a DGP with a continuous covariate where the impact on the outcome jumps at a threshold. For this, we change $W_{3}$ to a continuous $\mathcal{N}(40,10)$ and create a binary indicator for when $W_{3}>55$ and use this indicator in place of $W_{3}$ in the above outcome generation equation, but provide the EffectXshift algorithm the continuous covariate to determine whether it detects the correct threshold that leads to modification.

To establish the ground truth for the differential impacts of exposures, we applied a shift intervention to reduce $A_{1}$ by 0.5 units, i.e., we set $A_{1}\leftarrow A_{1}-0.5$, and recalculated the outcome $Y$ to obtain $Y_{\delta}$. Other exposures remained unchanged. The difference between the shifted and unshifted outcomes ($\Psi=Y_{\delta}-Y$) was aggregated to estimate the true average effect within each level of the covariate $W_{3}$, where the true effect modification was created. Here the true average in level 0 for $W_{3}$ is region $\mathcal{V}$ and level 1 is $\mathcal{V}^{c}$.

We assessed the bias, variance, mean squared error (MSE), and coverage of our method across different sample sizes ($n=300,500,1000,2000,5000$). For each sample size, we ran the simulation 100 times to evaluate the stability and reliability of our estimator. Bias: The average difference between the estimated and true effect. Variance: The variability of the estimator across simulations. MSE: Sum of the squared bias and variance, providing a comprehensive measure of estimator accuracy. Coverage: The proportion of times the true effect lies within the estimated 95% confidence interval.

These metrics are essential for ensuring the $\sqrt{n}$-consistent behavior of our estimator, meaning that as the sample size increases, the estimator converges to the true parameter value at a rate proportional to the square root of the sample size. This property is crucial for an efficient estimator, as it guarantees precise and reliable estimates in large samples.

For our estimator to be asympotically unbaised we need to accurately detect the correct region with increasing sample size. As such we report metrics for Accuracy, Precision, Recall and F1 statistics as well.

The simulation process involved the following steps: 1. Data Generation: For each sample size, we generated datasets of the above data. 2. Estimator Application: We applied our method to each dataset, to estimate the intervention effects within discovered subregions using 10-fold CV and default machine learning algorithms in the Super Learners (xgboost, random forest, glm) for the outcome regression and density estimation. 3. Metric Calculation: For each sample size and iteration, we calculated the bias, variance, MSE, and coverage of the estimated effects. 4. Result Aggregation: The results were aggregated to evaluate the overall performance of our method across different sample sizes.

For the binary effect modifier, for all iterations and sample sizes, we consistently identified the correct modifying variable, $W_{3}$ for the correct region, therefore accuracy, precision, recall and f1 were all 1.0.

The performance of our method in identifying the true subpopulation in continuous data improves as sample size increases. At a sample size of 300, the accuracy is 59.1%, precision is 18.9%, recall is 100%, and the F1 score is 29.9%. The high recall indicates perfect sensitivity, but the low precision suggests a high number of false positives. With a sample size of 500, accuracy improves to 94.2%, precision to 56.5%, recall remains at 100%, and the F1 score rises to 71.2%. The precision increase indicates fewer false positives. At 1000 samples, accuracy is 96.9%, precision is 81.8%, recall is 90.9%, and the F1 score is 82.9%. The balance between precision and recall reflects improved detection. For 2000 samples, accuracy is 99.7%, precision is 96.4%, recall is 99.3%, and the F1 score is 97.8%, showing minimal false positives and false negatives. At 5000 samples, the method nearly perfects performance with an accuracy of 99.9%, precision of 99.6%, recall of 99.1%, and an F1 score of 99.3%. Figure 1 shows these results.

The results are summarized below based on bias, variance, mean squared error (MSE), and coverage probability.

• •

  Bias and Variance: The bias and variance generally decreased as the sample size increased. For smaller samples (e.g., $n=300$), the bias was more pronounced in the measure in the $v^{c}$ region compared to the $v$ region, suggesting potential overfitting when sample sizes were limited. However, as the sample size increased, the bias significantly reduced, particularly evident in the $v^{c}$ region at $n=5000$ where the bias was minimal ($-0.007$).

• •

  Mean Squared Error (MSE): The MSE consistently decreased as the sample size increased, reflecting improvements in estimator accuracy. For instance, the MSE in the $v$ region decreased from 0.017 at $n=300$ to 0.0001 at $n=5000$. Similarly, in the $v^{c}$ region, the MSE decreased from 0.010 at $n=300$ to 0.0001 at $n=5000$.

• •

  Coverage Probability: The coverage probability was stablazlied to 0.95 for the $v$ region across sample sizes but $v^{c}$ region maintained a coverage probability of 1, indicating overly conservative estimates which may be due to limited iterations or our DGP.

The results of the continuous effect modifier simulation are summarized below based on bias, variance, mean squared error (MSE), and coverage probability.

• •

  Bias and Variance: As the sample size increased, both bias and variance generally decreased. For the $v$ region, the bias reduced from 0.155 at $n=300$ to a minimal 0.001 at $n=5000$. Similarly, in the $v^{c}$ region, the bias decreased from 1.025 at $n=300$ to 0.112 at $n=5000$, indicating significant improvement in accuracy with larger sample sizes.

• •

  Mean Squared Error (MSE): The MSE consistently decreased as the sample size increased, reflecting improvements in estimator accuracy. For the $v$ region, the MSE dropped from 0.026 at $n=300$ to 0.0007 at $n=5000$. In the $v^{c}$ region, the MSE decreased from 1.087 at $n=300$ to 0.014 at $n=5000$, demonstrating a substantial reduction in prediction error with larger samples.

• •

  Coverage Probability: The coverage probability showed consistent improvement with increasing sample size. For the $v$ region, coverage increased from 0.80 at $n=300$ to 0.95 at $n=5000$, indicating better confidence interval reliability. In the $v^{c}$ region, coverage improved from 0.72 at $n=300$ to 0.95 at $n=5000$, suggesting that the true parameter value was increasingly captured within the calculated confidence intervals as sample size grew. Variability in coverage may be due to limited iterations in our simulation.

Overall, the results suggest that our estimator’s mean squared error decreases at a rate proportional to $\sqrt{n}$, aligning with theoretical expectations.

We evaluated the impact of 16 persistent organic pollutants (POPs) on leukocyte telomere length (LTL). Our primary objective was to identify the baseline covariate region and the exposure with the maximum differential impact due to a reduction in exposure.

The shift interventions were applied to a subset of uncorrelated exposures: PCB74, 1,2,3,6,7,8-hxcdd, 1,2,3,4,6,7,8-hpcdd, 2,3,4,7,8-pncdf Lipid Adj, 1,2,3,4,6,7,8-hxcdf, and 3,3’,4,4’,5-pcnb. Age was consistently identified as the maximum effect modifier across all analysis folds. The exposure 3,3’,4,4’,5-pcnb (LBXPCBLA) was found in 100% of the folds, indicating a robust association. Table 1 shows the K-fold results. The change in outcome due to a 1 standard deviation reduction in 3,3’,4,4’,5-pcnb is shown under the column ”Effect”, inference is given in the standard error (SE) and confidence interval columns. ”Modifier” is the region in the covariate space and its complement and Fold is the fold number.

Age consistently emerged as the maximum effect modifier across all folds. The exposure 3,3’,4,4’,5-pcnb (LBXPCBLA) was identified in 100% of the folds. Our findings indicate that a reduction in exposure to 3,3’,4,4’,5-pcnb by one standard deviation (52.7907) leads to an increase in LTL. This effect is more pronounced in younger populations, with central ages of 11, 14, and 17 consistently identified as split points. All these results were statistically significant based on the confidence intervals, underscoring the differential impact of this exposure across age.

For each region and its complement, we provide the pooled oracle estimates. Here, $v$ denotes the region for the lower age group and $v^{c}$ for the higher age group. Below is Table 2 of these pooled results.

These estimates indicate that shifting 3,3’,4,4’,5-pcnb by one standard deviation (approximately -52.8) leads to an estimated increase in telomere length of about 0.488 in the lower age group ($v$) and about 0.316 in the higher age group ($v^{c}$).

To formally assess whether the difference between these two subgroups is statistically significant, we can compare the estimates and their variances directly. Letting $\hat{\psi}_{\mathcal{V}}=\hat{\Psi}_{v}-\hat{\Psi}_{v^{c}}=0.4879-0.3156=0.1723$, and assuming negligible covariance between the estimates, the variance of the difference is $\widehat{\mathrm{Var}}(\hat{\psi}_{\mathcal{V}})=0.0001027+0.0010592=0.0011619$. Thus, $\text{SE}(\hat{\psi}_{\mathcal{V}})=\sqrt{0.0011619}\approx 0.0341$.

The corresponding z-statistic is:

which yields a p-value well below 0.001. This confirms that the difference in the effect estimates between the two subpopulations is highly statistically significant.