Understanding the pathways whereby an intervention has an effect on an outcome is a common scientific goal. A rich body of literature provides various decompositions of the total intervention effect into pathway-specific effects. Interventional direct and indirect effects provide one such decomposition. Existing estimators of these effects are based on parametric models with confidence interval estimation facilitated via the nonparametric bootstrap. We provide theory that allows for more flexible, possibly machine learning-based, estimation techniques to be considered. In particular, we establish weak convergence results that facilitate the construction of closed-form confidence intervals and hypothesis tests and prove multiple robustness properties of the proposed estimators. Simulations show that inference based on large-sample theory has adequate small-sample performance. Our work thus provides a means of leveraging modern statistical learning techniques in estimation of interventional mediation effects.
Recent advances in causal inference have provided rich frameworks for posing interesting scientific questions pertaining to the mediation of effects through specific biologic pathways (Yuan and MacKinnon , Imai et al. , Valeri and VanderWeele , Pearl , Naimi et al. , Zheng and van der Laan , VanderWeele and Tchetgen Tchetgen , among others). Foremost among these advances is the provision of model-free definitions of mediation parameters, which enables researchers to develop robust estimators of these quantities. A debate in this literature has emerged pertaining to the reliance of methodology on cross-world independence assumptions that are fundamentally untestable even in randomized controlled experiments [8,9,10]. One approach to this problem is to utilize methods that attempt to estimate bounds on effects (Robins and Richardson , Tchetgen Tchetgen and Phiri , among others). A second approach considers seeking alternative definitions of mediation parameters that do not require such cross-world assumptions (VanderWeele et al. , Rudolph et al. , among others). Rather than considering deterministic interventions on mediators (i.e., a hypothetical intervention that fixes every individual mediator to a particular value), these approaches consider stochastic interventions on mediators (i.e., hypothetical interventions where the mediator is drawn from a particular conditional distribution). In this class of approaches, that of Vansteelandt and Daniel  is particularly appealing. Building on the prior work of VanderWeele et al. , the authors provide a simple decomposition of the total effect into direct effects and pathway-specific effects via multiple mediators. Interestingly, their decompositions hold even when the structural dependence between mediators is unknown.
Vansteelandt and Daniel  described two approaches to estimation of the effects using parametric working models for relevant nuisance parameters. In both cases, the nonparametric bootstrap was recommended for inference. A potential limitation of the proposal is that correctly specifying a parametric working model may be difficult in many settings. In these instances, we may rely on flexible estimators of nuisance parameters, for example, based on machine learning. When such techniques are employed, the nonparametric bootstrap does not generally guarantee valid inference . This fact motivates the present work, where we develop nonparametric efficiency theory for the interventional mediation effect parameters. This theory allows us to utilize frameworks for nonparametric efficient inference to develop estimators of the quantities of interest. We propose a one-step and a targeted minimum loss-based estimator and demonstrate that under suitable regularity conditions, both estimators are nonparametric efficient among the class of regular asymptotically linear estimators. The estimators also enjoy a multiple robustness property, which ensures consistency of effect estimates if at least some combinations of nuisance parameters are consistently estimated. Another benefit enjoyed by our estimators is the availability of closed-form confidence intervals and hypothesis tests.
2 Interventional effects
Adopting the notation of Vansteelandt and Daniel , suppose the observed data are represented as independent copies of the random variable , where is a vector of confounders, is a binary intervention, and are mediators, and is a relevant outcome. Our developments pertain to both discrete and real-valued mediators, while without loss of generality, we assume . We assume ; that is, any subgroup defined by covariates that is observed with positive probability should have some chance of receiving both interventions. We also assume that for , the probability distribution of given has density with respect to some dominating measures and this density satisfies , where the infimum is taken over . Similarly, we assume that for all . Beyond these conditions, encodes no assumptions about ; however, the efficiency theory that we develop still holds under a model that makes assumptions about , including the possibility that this quantity is known exactly, as in a stratified randomized trial.
To define interventional mediation effects, notation for counterfactual random variables is required. For , and , let denote the counterfactual value for the th mediator when is set to . Similarly, let denote the counterfactual outcome under an intervention that sets and . As a point of notation, when introducing quantities whose definition depends on particular components of the random variable , we will use lower case letters to denote the particular value and assume that the definition at hand applies for all values in the support of that random variable.
The total effect of intervening to set versus is , where we use to emphasize that we are taking an expectation with respect to a distribution of a counterfactual random variable. The total effect describes the difference in counterfactual outcome considering an intervention where we set and allow the mediators to naturally assume the value that they would under intervention versus an intervention where we set and allow the mediators to vary accordingly. To contrast with forthcoming effects, it is useful to write the total effect in integral form. Specifically, we use to denote the covariate-conditional mean of the counterfactual outcome , to denote the covariate-conditional bivariate cumulative distribution function of , and to denote the marginal distribution of . The total effect can be written as
The total effect can be decomposed into interventional direct and indirect effects. The interventional direct effect is the difference in average counterfactual outcome under two population-level interventions. The first intervention sets , and subsequently for individuals with draws mediators from . Thus, on a population level the covariate conditional distribution of mediators in this counterfactual world is the same as it would be in a population where everyone received intervention . This is an example of a stochastic intervention . The second intervention sets and subsequently allows the mediators to naturally assume the value that they would under intervention , so that the population level mediator distribution is again . The interventional direct effect compares the average outcome under these two interventions,
For interventional indirect effects, we require definitions for the covariate-conditional distribution of each mediator, which we denote for by . The interventional indirect effect through is
As with the direct effect, this effect considers two interventions. Both interventions set . The first intervention draws mediator values independently from the marginal mediator distributions and , while the second intervention draws mediator values independently from the marginal mediator distributions and . The effect thus describes the average impact of shifting the population level distribution of , while holding the population level distribution of fixed. The interventional indirect effect on the outcome through is similarly defined as
Note that when defining interventional indirect effects, mediators are drawn independently from marginal mediator distributions. The final effect in the decomposition essentially describes the impact of drawing the mediators from marginal rather than joint distributions. Thus, we term this effect the covariant mediator effect, defined as
where . Vansteelandt and Daniel  discussed situations where these effects are of primary interest.
From the aforementioned definitions, we have the following effect decomposition . These component effects can be identified using the observed data under the following assumptions:
the effect of on is unconfounded given , ;
the effect of and on is unconfounded given and , ;
the effect of on is unconfounded given , .
The identifying formula for each effect can now be written as a statistical functional of the observed data distribution by substituting the outcome regression for and the observed-data mediator distributions for the respective counterfactual distributions in the aforementioned integral expressions.
We note that the aforementioned assumptions preclude the existence of treatment-induced confounding of the mediator-outcome association. In the Supplementary material, we provide relevant extensions to this setting.
3.1 Efficiency theory
In this section, we develop efficiency theory for nonparametric estimation of interventional effects. This theory centers around the efficient influence function of each parameter. The efficient influence function is important for several reasons. First, it allows us to utilize two existing estimation frameworks, one-step estimation [17,18] and targeted minimum loss-based estimation [19,20], to generate estimators that are nonparametric efficient. That is, under suitable regularity conditions, they achieve the smallest asymptotic variance among all regular estimators that, when scaled by , have an asymptotic normal distribution. We discuss how these estimators can be implemented in Section 3.2. The second important feature of the efficient influence function is that its variance equals the variance of the limit distribution of the scaled estimators. Thus, an estimate of the variance of the efficient influence function is a natural standard error estimate, which affords closed-form Wald-style confidence intervals and hypothesis tests (Section 3.3). Finally, the efficient influence function also characterizes robustness properties of our proposed estimators (Section 3.4).
To introduce the efficient influence function, several additional definitions are required. For a given distribution , we define , commonly referred to as a propensity score. For and , we introduce the following partially marginalized outcome regressions, . We also introduce notation for the indicator function defined by if and zero otherwise. is similarly defined.
Under sampling from , the efficient influence function evaluated on a given observation for the total effect is
The efficient influence function for the interventional direct effect is
The efficient influence function for the interventional indirect effect through is
The efficient influence function for the interventional indirect effect through is
The efficient influence function for the covariant interventional effect is .
A proof of Theorem 1 is provided in the Supplementary material.
We propose estimators of each interventional effect using one-step and targeted minimum loss-based estimation. Both techniques develop along a similar path. We first obtain estimates of the propensity score, outcome regression, and joint mediator distribution; we collectively refer to these quantities as nuisance parameters. With estimated nuisance parameters in hand, we subsequently apply a correction based on the efficient influence function to the nuisance estimates.
To estimate the propensity score, we can use any suitable technique for mean regression of the binary outcome onto confounders . Working logistic regression models are commonly used for this purpose, though semi- and nonparametric alternatives would be more in line with our choice of model. We denote by the chosen estimate of . Similarly, the outcome regression can be estimated using mean regression of the outcome onto and . For example, if the study outcome is binary, logistic regression could again be used, though more flexible regression estimators may be preferred. As above, we denote by the estimated outcome regression evaluated under , with providing an estimate of . To estimate the marginal cumulative distribution of , we will use the empirical cumulative distribution function, which we denote by .
Estimation of the conditional joint distribution of the mediators is a more challenging proposition, as fewer tools are available for flexible estimation of conditional multivariate distribution functions. We hence focus our developments on the development of approaches for discrete-valued mediators. The approach we adopt could be extended to continuous-valued mediators by considering a fine partitioning of the mediator values. We examine this approach via simulation in Section 4. To develop our density estimators, we use the approach of Dìaz Muñoz and van der Laan , which considers estimation of a conditional density via estimation of discrete conditional hazards. Briefly, consider estimation of the distribution of given and , and, for simplicity, suppose that the support of is . We create a long-form data set, where the number of rows contributed by each individual contribute is equal to their observed value of . An example is illustrated in Table 1. We see that the long-form data set includes an integer-valued column named “bin” that indicates to which value of each row corresponds, as well as a binary column indicating whether the observed value of corresponds to each bin. These long-form data can be used to fit a regression of the binary outcome onto , , and bin. This naturally estimates , the conditional discrete hazard of given and . Let denote the estimated hazard obtained from fitting this regression. An estimate of the density at is
Similarly, an estimate of the conditional distribution of given can be obtained. An estimate of the joint conditional density is implied by these estimates, , while an estimate of the marginal distribution of is .
An ID is uniquely assigned to each independent data unit and a single confounder is included in the mock data set.
In principle, one could reverse the roles of and in the above procedure. That is, we could instead estimate the distribution of given and of given . Cross-validation could be used to pick between the two potential estimators of the joint distribution. Other approaches to conditional density estimation are permitted by our procedure as well. For example, approaches based on working copula models may be particularly appealing in this context, as they allow separate specification of marginal vs joint distributions of the mediators.
Given estimates of nuisance parameters, we now illustrate one-step estimation for the interventional direct effect. One-step estimators of other effects can be generated similarly. A plug-in estimate of the conditional interventional direct effect given is the difference between
To obtain a plug-in estimate of , we standardize the conditional effect estimate with respect to , the empirical distribution of . Thus, the plug-in estimator of is .
The one-step estimator is constructed by adding an efficient influence function-based correction to an initial plug-in estimate. Suppose we are given estimates of all relevant nuisance quantities and let denote any probability distribution in that is compatible with these estimates. The efficient influence function for under sampling from is , and the one-step estimator is . All other effect estimates are generated in this vein: estimated nuisance parameters are plugged in to the efficient influence function, the resultant function is evaluated on each observation, and the empirical average of this quantity is added to the plug-in estimator.
While one-step estimators are appealing in their simplicity, the estimators may not obey bounds on the parameter space in finite samples. For example, if the study outcome is binary, then the interventional effects each represent a difference in two probabilities and thus are bounded between and 1. However, one-step estimators may fall outside of this range. This motivates estimation of these quantities using targeted minimum loss-based estimation, a framework for generating plug-in estimators. The implementation of such estimators is generally more involved than that of one-step estimators. In this approach, a second-stage model fitting is used to ensure that nuisance parameter estimates satisfy efficient influence function estimating equations. The approach for this second-stage fitting is dependent on the specific effect parameter considered and the procedure differs subtly for the various effect measures presented here. Supplementary material includes a detailed exposition of how such estimators can be implemented.
3.3 Large sample inference
We now present a theorem establishing the joint weak convergence of the proposed estimators to a random variable with a multivariate normal distribution. Because the asymptotic behavior of the one-step and targeted minimum loss estimators (TMLEs) is equivalent, we present a single theorem. A discussion of the differences in regularity conditions required to prove the theorem for one-step versus targeted minimum loss estimation is provided in the Supplementary material. Let denote the vector of (one-step or targeted minimum loss) estimates of and let denote the vector of efficient influence functions defined by
In the theorem, we use to denote the -norm, define for any -measurable as .
Under sampling from , if for ,
in probability as ,
in probability as ,
, , and
in probability as and falls in a -Donsker class with probability tending to 1,
The regularity conditions required for Theorem 2 are typical of many problems in semiparametric efficiency theory. We provide conditions in terms of -norm convergence, as this is typical of this literature; however, alternative and potentially weaker conditions are possible to derive. For further discussion, see the supplementary material. As with any nonparametric procedure, there is a concern related to the dimensionality , particularly in situations with real-valued mediators. Minimum loss estimators (MLEs) in certain function classes can attain the requisite convergence rates. For example, an MLE in the class of functions that are right-continuous with left limits (i.e., càdlàg) with variation norm bounded by a constant achieves an convergence rate faster than irrespective of the dimension of the conditioning set . However, this may not allay all concerns pertaining to the curse of dimensionality due to the fact that in moderately high dimensions, these function classes can be restrictive and thus the true function may fall outside this class. Nevertheless, we suggest (and our simulations show) that in spite of concerns pertaining to the curse of dimensionality our procedure will enjoy reasonable finite-sample performance in many settings.
The covariance matrix may be estimated by the empirical covariance matrix of the vector applied to the observed data, where is any distribution in the model that is compatible with the estimated nuisance parameters. With the estimated covariance matrix, it is straightforward to construct Wald confidence intervals and hypothesis tests about the individual interventional effects or comparisons between them. For example, a straightforward application of the delta method would allow for a test of the null hypothesis that .
3.4 Robustness properties
As with many problems in causal inference, consistent estimation of interventional effects requires consistent estimation only of certain combinations of nuisance parameters. To determine these combinations, we may study the stochastic properties of the efficient influence function. In particular, consider a parameter whose value under is and whose efficient influence function under sample from can be written , where is the value of the parameter of interest under . Then we may study the circumstances under which . This generally entails understanding which parameters of must align with those parameters of to ensure that the influence function has mean zero under sampling from . We present the results of this analysis in a theorem below and refer readers to the Supplementary material for the proof.
Locally efficient estimators of the total effect and the intervention direct, indirect, and covariant effects are consistent for their respective target parameters if the following combinations of nuisance parameters are consistently estimated:
Total effect: or
Interventional direct effect: or or ;
Interventional indirect effect through : or or or ;
Interventional indirect effect through : or or or ;
Interventional covariant effect: or or .
The most interesting robustness result is perhaps that pertaining to the indirect effects. The first condition for consistent estimation is expected, as the propensity score plays no role in the definition of the indirect effect. The second condition shows that the joint mediator distribution and propensity score together can compensate for inconsistent estimation of the outcome regression, while the relevant marginal mediator distributions are required to properly marginalize the resultant quantity. The third and fourth conditions show that inconsistent estimation of the marginal distribution one, but not both, of the mediators can be corrected for via the propensity score.
We note that Theorem 3 provides sufficient, but not necessary, conditions for consistent estimation of each effect. For example, a consistent estimate of the total effect is implied by a consistent estimate of and , a condition that is generally weaker than requiring consistent estimation of the outcome regression and joint mediator distribution. Because our estimation strategy relies on estimation of the joint mediator distribution, we have described robustness properties in terms of the large sample behavior of estimators of those quantities.
In the Supplementary material, we provide relevant extensions to the setting where the mediator–outcome relationship is confounded by measured covariates whose distributions are affected by the treatment. In this case, both the effects of interest and their efficient influence functions involve the conditional distribution of the confounding covariates. We discuss the relevant modifications to the estimation procedures to accommodate this setting in the supplement.
Generalization to other effect scales requires only minor modifications. First, we determine the portions of the efficient influence function that pertain to each component of the additive effect. For example, considering , we identify the portions of the efficient influence function that pertain to the mean counterfactual under draws of from and of from versus those portions that pertain to the mean counterfactual under draws of from and of from . We then develop a one-step or TMLE for each of these components separately. Finally, we use the delta method to derive the resulting influence function. In the Supplementary material, we illustrate an extension to a multiplicative scale.
Our results can also be extended to estimation of interventional effects for more than two mediators. As discussed in Vansteelandt and Daniel , when there are more than two mediators, say , there are many possible path-specific effects. However, our scientific interest is usually restricted to learning effects that are mediated through each of the mediators, rather than all possible path-specific effects. Moreover, strong untestable assumptions are required to infer all path-specific effects, including assumptions about the direction of the causal effects between mediators. Therefore, it may be of greatest interest to evaluate direct effects such as
which describes the effect of setting versus , while drawing all mediators from the joint conditional distribution given , and for , indirect effects such as
which describes the effect of setting to the value it would assume under versus while drawing from their respective marginal distributions given and drawing from their marginal distribution given . We provide relevant efficiency theory for these parameters in the Supplementary material.
4.1 Discrete mediators
We evaluated the small sample performance of our estimators via Monte Carlo simulation. Data were generated as follows. We simulated by drawing independently from Uniform(0,1) distributions, and independently from Bernoulli distributions with success probability of 0.25 and 0.5, respectively. The treatment variable , given , was drawn from a Bernoulli distribution with and . Here, we consider and . Given , the first mediator was generated by taking draws from a geometric distribution with success probability . Any draw of six or greater was set equal to six. The second mediator was generated from a similarly truncated geometric distribution with success probability . Given , the outcome was drawn from a Bernoulli distribution with success probability . The mediator distribution is visualized for combinations of and in Figure 1. The true total effect is approximately 0.06, which decomposes into a direct effect of 0.05, an indirect effect through of , an indirect effect through of 0.02, and a covariant effect of 0.
The nuisance parameters were estimated using regression stacking [23,24], also known as super learning  using the SuperLearner package for the R language . We used this package to generate an ensemble of a main-terms logistic regression (as implemented in the SL.glm function in SuperLearner), polynomial multivariate adaptive regression splines (SL.earth), and a random forest (SL.ranger). The ensemble was built by selecting the convex combination of these three estimators that minimized tenfold cross-validated deviance.
We evaluated our proposed estimators under this data generating process at sample sizes of 250, 500, 1,000, and 2,000. At each sample size, we simulated 1,000 data sets. Point estimates were compared in terms of their Monte Carlo bias, standard deviation, and mean squared error. We evaluated weak convergence by visualizing the sampling distribution of the estimators after centering at the true parameter value and scaling by an oracle standard error, computed as the Monte Carlo standard deviation of the estimates, as well as scaling by an estimated standard error based on the estimated variance of the efficient influence function. Similarly, we evaluated the coverage probability of a nominal 95% Wald-style confidence interval based on the oracle and estimated standard errors.
In terms of estimation, one-step and TMLE behave as expected in large samples (Figure 2). The estimators are approximately unbiased in large samples and have mean squared error appropriately decreasing with sample size. Comparing the two estimation strategies, we see that one-step and TMLEs had comparable performance for the interventional direct effect, while the TMLE had better performance for the indirect effects. However, the one-step was uniformly better for estimating the covariant effect owing to large variability of the TMLE of this quantity. Further examination of the results revealed that the second-stage model fitting required by the targeted minimum loss approach was could be unstable in small samples, leading to extreme results in several data sets.
The sampling distributions of the centered and scaled estimators were approximately a standard normal distribution (Figures 3 and 4), excepting the TMLE scaled by an estimated standard error. Confidence intervals based on an oracle standard error came close to nominal coverage in all sample sizes, while those based on an estimated standard error tended to have under-coverage in small samples.
4.2 Continuous mediators
We examined the impact of discretization of the mediator distributions when in fact the mediators are continuous valued. To that end, we simulated data as follows. Covariates were simulated as above. The treatment variable given was drawn from a Bernoulli distribution with . Given , and were, respectively, drawn from normal distributions with unit variance and mean values and . As above, Super Learner was used to estimate all nuisance parameters. To accommodate appropriate modeling of the interactions, we replaced the main terms GLM (SL.glm) with a forward stepwise GLM algorithm that included all two-way interactions (SL.step.interaction). The true effect sizes were approximately the same as in the first simulation. We evaluated discretization of each continuous mediator distribution into 5 and 10 evenly spaced bins. For the sake of space, we focus results on the one-step estimator; results for TMLE are included in the supplement.
Overall, discretization of the continuous mediator distribution had a greater impact on the performance of indirect effect estimators compared to direct effects (Figures 5 and 6). For the latter effects, oracle confidence intervals for both levels of discretization achieved nominal coverage for all sample sizes considered. For the indirect effects, we found that there was non-negligible bias in the estimates due to the discretization. The impacts in terms of confidence interval coverage were minimal in small sample sizes, but lead to under-coverage in larger sample sizes. Including more bins generally lead to better performance, but these estimates still exhibited bias in the largest sample sizes that impacted coverage. Nevertheless the performance of the indirect effect estimators was reasonable with oracle coverage 90% for all sample sizes.
4.3 Additional simulations
In the Supplementary material, we include several additional simulations studying the impact of the number of levels of the discrete mediator, as well as the impact of inconsistent estimation of the various nuisance parameters. For the former, we found that the results of the simulation were robust to number of mediator levels in the setting considered. For the latter, we confirmed the multiple robustness properties of the indirect effect estimators by studying the bias and standard deviation of the estimators in large sample sizes under the various patterns of misspecification given in our theorem.
Our simulations demonstrate adequate performance of the proposed nonparametric estimators of interventional mediation effects in settings with relatively low-dimensional covariates (five, in our simulation). In certain settings, it may only be necessary to adjust for a limited number of covariates to adequately control confounding. For example, in the study of the mediating mechanisms of preventive vaccines using data from randomized trials, we need to only adjust for confounders of the mediator/outcome relationship, since other forms of confounding are addressed by the randomized design. Generally, there are few known factors that are likely to impact vaccine-induced immune responses and so nonparametric analyses may be quite feasible in this case. For example, Cowling et al.  studied mediating effects of influenza vaccines, adjusting only for age. Thus, we suggest that interventional mediation estimands and nonparametric estimators thereof may be of interest for studying mediating pathways of vaccines. However, in other scenarios, it may be necessary to adjust for a high-dimensional set of confounders. For example, in observational studies of treatments (e.g., through an electronic health records system), we may require control for a high-dimensional set of putative confounders of treatment and outcome. This may raise concerns related to the curse-of-dimensionality when utilizing nonparametric estimators. Studying tradeoffs between the selection of various estimation strategies in this context will be an important area for future research.
We have developed an R package intermed with implementations of the proposed methods, which is included in the Supplementary material. The package focuses on implementations for discrete mediators. However, our simulations demonstrate a clear need to extend the software to accommodate adaptive selection of the number of bins in the mediator density estimation procedure for continuous mediators. In small sample sizes, we found that course binning leads to adequate results, but as sample size increased, unsurprisingly there was a need for finer partitioning to reduce bias. In future versions of the software, we will include such adaptive binning strategies, as well as other methods for estimating continuous mediator densities.
The behavior of the TMLE of the covariant effect in the simulation is surprising as generally we see comparable or better performance of such estimators relative to one-step estimators. This can likely be attributed to the fact that the targeted minimum loss procedure does not yield a compatible plug-in estimator of the vector , in the sense that there is likely no distribution that is compatible with all of the various nuisance estimators after the second-stage model fitting. A more parsimonious approach could consider either an iterative targeting procedure or a uniformly least favorable submodel that simultaneously targets the joint mediator density and outcome regression. The former is implemented in a concurrent proposal , where one-step and TMLEs of interventional effects are developed for a single mediator when the mediator–outcome relationship is subject to treatment-induced confounding. In their set up, if one can treat the treatment-induced confounder as a second mediator, then their proposal results in an estimate of one component of our indirect effect. In their simulations, they find superior finite-sample performance of the TMLE relative to the one step, suggesting that targeting the mediator densities may be a more robust approach. However, their simulation involved only binary-valued mediators, so further comparison of these approaches is warranted in settings similar to our simulation, where mediators can take many values. We leave these developments to future work.
The Donsker class assumptions of our theorem could be removed by considering cross-validated nuisance parameter estimates (also known as cross-fitting) [29,30]. This technique is implemented in our R package, but we leave to future research the examination of its impact of estimation and inference. We hypothesize that this approach will generally improve the anti-conservative confidence intervals in small samples, but will have little impact on performance of point estimates in terms of bias and variance.
Funding information: D. Benkeser was funded by National Institutes of Health award R01AHL137808 and National Science Foundation award 2015540. Code to reproduce simulation results is available at https://github.com/benkeser/intermed/tree/master/simulations.
Conflict of interest: Prof. David Benkeser is a member of the Editorial Board in the Journal of Causal Inference but was not involved in the review process of this article.
 Yuan Y, MacKinnon DP. Bayesian mediation analysis. Psychol Methods. 2009;14(4):301. Search in Google Scholar
 Imai K, Keele L, Tingley D. A general approach to causal mediation analysis. Psychol Methods. 2010;15(4):309. Search in Google Scholar
 Valeri L, VanderWeele TJ. Mediation analysis allowing for exposure-mediator interactions and causal interpretation: theoretical assumptions and implementation with SAS and SPSS macros. Psychol Methods. 2013;18(2):137. Search in Google Scholar
 Pearl J. Interpretation and identification of causal mediation. Psychol Methods. 2014;19(4):459. Search in Google Scholar
 Naimi AI, Schnitzer ME, Moodie EE, Bodnar LM. Mediation analysis for health disparities research. Am J Epidemiol. 2016;184(4):315–24. Search in Google Scholar
 Zheng W, van der Laan MJ. Longitudinal mediation analysis with time-varying mediators and exposures, with application to survival outcomes. J Causal Infer. 2017;5(2):20160006. Search in Google Scholar
 VanderWeele TJ, Tchetgen Tchetgen EJ. Mediation analysis with time varying exposures and mediators. J R Stat Soc B (Statistical Methodology). 2017;79(3):917–38. Search in Google Scholar
 Dawid AP. Causal inference without counterfactuals. J Am Stat Assoc. 2000;95(450):407–24. Search in Google Scholar
 Pearl J. Direct and indirect effects. In: Proceedings of the 17th Conference in Uncertainty in Artificial Intelligence. San Francisco: Morgan Kaufmann; 2001. Search in Google Scholar
 Robins JM, Richardson TS. Alternative graphical causal models and the identification of direct effects. Causality and psychopathology: Finding the determinants of disorders and their cures. Oxford, New York: Oxford University Press; 2011. p. 103–58. Search in Google Scholar
 Tchetgen Tchetgen EJ, Phiri K. Bounds for pure direct effect. Epidemiology (Cambridge, Mass.). 2014;25(5):775. Search in Google Scholar
 VanderWeele TJ, Vansteelandt S, Robins JM. Effect decomposition in the presence of an exposure-induced mediator-outcome confounder. Epidemiology. 2014;25(2):300. Search in Google Scholar
 Rudolph KE, Sofrygin O, Zheng W, van der Laan MJ. Robust and flexible estimation of stochastic mediation effects: a proposed method and example in a randomized trial setting. Epidemiol Methods. 2018;7(1):2017007. Search in Google Scholar
 Vansteelandt S, Daniel RM. Interventional effects for mediation analysis with multiple mediators. Epidemiology. 2017;28(2):258. Search in Google Scholar
 Coyle J, van der Laan MJ. Targeted bootstrap. In Targeted learning for data science. Cham: Springer International Publishing; 2018. p. 523–39. Ch. 28. Search in Google Scholar
 Muñoz ID, van der Laan MJ. Population intervention causal effects based on stochastic interventions. Biometrics. 2012;68(2):541–9. Search in Google Scholar
 Ibragimov I, Khasminskii R. Statistical estimation: asymptotic theory. New York: Springer-Verlag; 1981. Search in Google Scholar
 Bickel P, Klaassen C, Ritov Y, Wellner J. Efficient and adaptive estimation for semiparametric models. Berlin Heidelberg New York: Springer; 1997. Search in Google Scholar
 van der Laan M, Rubin DB. Targeted maximum likelihood learning. Int J Biostat. 2006;2(1):11. Search in Google Scholar
 van der Laan M, Rose S. Targeted learning: causal inference for observational and experimental data. Berlin Heidelberg New York: Springer; 2011. Search in Google Scholar
 Dìaz Muñoz I, van der Laan MJ. Super learner based conditional density estimation with application to marginal structural models. Int J Biostat. 2011;7(1):1–20. Search in Google Scholar
 Benkeser D, van der Laan MJ. The highly adaptive lasso estimator. In: 2016 IEEE International Conference on Data Science and Advanced Analytics (DSAA). IEEE; 2016. p. 689–96. Search in Google Scholar
 Wolpert DH. Stacked generalization. Neural Netw. 1992;5:241–59. Search in Google Scholar
 Breiman L. Stacked regressions. Mach Learn. 1996;24:49–64. Search in Google Scholar
 van der Laan M, Polley E, Hubbard A. Super learner. Stat Appl Genet Mol. 2007;6(1):25. Search in Google Scholar
 Polley E, LeDell E, Kennedy C, van der Laan MJ. SuperLearner: Super Learner Prediction. R package version 2.0-28; 2013. https://CRAN.R-project.org/package=SuperLearnerSearch in Google Scholar
 Cowling BJ, Lim WW, Perera RA, Fang VJ, Leung GM, Peiris JM, et al. Influenza hemagglutination-inhibition antibody titer as a mediator of vaccine-induced protection for influenza B. Clin Infect Dis. 2019;68(10):1713–17. Search in Google Scholar
 Zheng W, van der Laan MJ. Asymptotic theory for cross-validated targeted maximum likelihood estimation. Technical Report 273. Berkeley: Division of Biostatistics, University of California, Berkeley; 2010. Search in Google Scholar
 Chernozhukov V, Chetverikov D, Demirer M, Duflo E, Hansen C, Newey W, et al. Double/debiased machine learning for treatment and structural parameters. Econom J. 2018;21(1):C1–C68. Search in Google Scholar
© 2021 David Benkeser and Jialu Ran, published by De Gruyter
This work is licensed under the Creative Commons Attribution 4.0 International License.