Robust and Flexible Estimation of Stochastic Mediation Effects: A Proposed Method and Example in a Randomized Trial Setting

1 Introduction

Mediation allows for an examination of the mechanisms driving a relationship. Much of epidemiology entails reporting exposure-outcome associations where the exposure may be multiple steps removed from the outcome. For example, risk-factor epidemiology demonstrated that obesity increases the risk of type 2 diabetes, but biochemical mediators linking the two have advanced our understanding of the causal relationship (Kahn, Hull, and Utzschneider 2006). Mediators have been similarly important in understanding how social exposures act to affect health outcomes. In the illustrative example we consider in this paper, the Moving to Opportunity (MTO) experiment randomized families living in public housing to receive a voucher that they could use to rent housing on the private market, which reduced their exposure to neighborhood poverty (Kling, Liebman, and Katz 2007). Ultimately, being randomized to receive a voucher affected subsequent adolescent drug use (Orr et al. 2003; Rudolph et al. 2017). In the illustrative example, we test the extent to which the effect operates through a change in the adolescent’s school environment.

Causal mediation analysis (Imai, Keele, and Tingley 2010; Valeri and VanderWeele 2013; Zheng and van der Laan 2012b) (also called mediation analysis using the counterfactual framework (Valeri and VanderWeele 2013)) shares similar goals with the standard mediation approaches, e.g., structural equation modeling and the widely used Baron and Kenny “product method” approach (Baron and Kenny 1986; Valeri and VanderWeele 2013). They all aim to test mechanisms and estimate direct and indirect effects. Advantages of causal mediation analysis include that estimates have a causal interpretation (under specified identifying assumptions) and some approaches make fewer restrictive parametric modeling assumptions. For example, in contrast to traditional approaches, approaches within the causal mediation framework (1) allow for interaction between the treatment and mediator (VanderWeele 2009), (2) allow for modeling nonlinear relationships between mediators and outcomes (VanderWeele 2009), and (3) allow for incorporation of data-adaptive machine learning methods and double robust estimation (Zheng and van der Laan 2012b; Tchetgen and Shpitser 2014).

However, despite these advantages, the assumptions required to estimate certain causal mediation effects may sometimes be untenable; for example, the assumption that there is no confounder of the mediator-outcome relationship that is affected by treatment (in the literature, such a confounder is referred to as confounding by a causal intermediate (Petersen, Sinisi, and van der Laan 2006), a time-varying confounder affected by prior exposure (VanderWeele and Tchetgen Tchetgen 2017), or time-dependent confounding by an intermediate covariate (van der Laan and Gruber 2012)). For brevity, we will refer to such a variable as an intermediate confounder.

There have been recently proposed causal mediation estimands, called randomized (i.e., stochastic) interventional direct effects and interventional indirect effects that do not require this assumption (Didelez, Dawid, and Geneletti 2006; van der Laan and Petersen 2008; Zheng and van der Laan 2012a; VanderWeele, Vansteelandt, and Robins 2014; Vansteelandt and Daniel 2017; VanderWeele and Tchetgen Tchetgen 2017; Zheng and van der Laan 2017). We build on this work, proposing a robust and flexible estimator for these effects, which we call stochastic direct and indirect effects (SDE and SIE).

This paper is organized as follows. In the following section, we review and compare common causal mediation estimands, providing the assumptions necessary for their identification. Then, we describe our proposed estimator, its motivation, and its implementation in detail. Code to implement this method is provided in the Appendix. We then provide results from a limited simulation study demonstrating the estimator’s finite sample performance and robustness properties. Lastly, we apply the method in a longitudinal, randomized trial setting.

2 Notation and causal mediation estimands

Let observed data: O=(W,A,Z,M,Y) with n i.i.d. copies O1,...,On∼P0, where W is a vector of pre-treatment covariates, A is the treatment, Z is the intermediate confounder affected by A, M is the mediator, and Y is the outcome. For simplicity, we assume that A,Z,M,and Y are binary. In our illustrative example, A is an instrument, so it is reasonable to assume that M and Y are not affected by A except through its effect on Z. Mirroring the structural causal model (SCM) of our illustrative example, we assume that M is affected by {Z,W} but not A, and that Y is affected by {M,Z,W} but not A. We assume exogenous random errors: (UW,UA,UZ,UM,UY). This SCM is represented in Figure 1 and the following causal models: W=f(UW), A=f(UA), Z=f(A,W,UZ), M=f(Z,W,UM), and Y=f(M,Z,W,UY). Note that this SCM (including that UY and UM are not affected by A) puts the following assumptions on the probability distribution: P(Y|M,Z,A,W)=P(Y|M,Z,W) and P(M|Z,A,W)=P(M|Z,W). However, our approach generalizes to scenarios where A also affects M and Y as well as to scenarios where A is not random. We provide details and discuss these generalizations in the Appendix. We can factorize the likelihood for the SCM reflecting our illustrative example as follows: P(O)=P(Y|M,Z,W)P(M|Z,W)P(Z|A,W)P(A)P(W).

Figure 1:

Structural causal model reflecting the illustrative example.

Causal mediation analysis typically involves estimating one of two types of estimands: controlled direct effects (CDE) or natural direct and indirect effects (NDE, NIE). Controlled direct effects involve comparing expected outcomes under different values of the treatment and setting the value of the mediator for everyone in the sample. For example the CDE can be defined: E(Ya,m−Ya∗,m), where Ya,m, Ya∗,m, is the counterfactual outcome setting treatment A equal to a or a∗, respectively (the two treatment values being compared), and setting mediator M equal to m. In contrast, the NDE can be defined: E(Ya,Ma∗−Ya∗,Ma∗), where Ya,Ma∗, Ya∗,Ma∗ is the counterfactual outcome setting A equal to a or a∗ but this time setting Ma∗ to be the counterfactual value of the mediator had A been set to a∗ (possibly contrary to fact). Similarly, the NIE can be defined: E(Ya,Ma−Ya,Ma∗). Natural direct and indirect effects are frequently used in epidemiology and have the appealing property of adding to the total effect (Pearl 2001).

Although the NDE and NIE are popular estimands, their identification assumptions may sometimes be untenable. Broadly, identification of their causal effects relies on the sequential randomization assumption on intervention nodes A and M and positivity. Two specific ignorability assumptions are required to identify CDEs and NDE/NIEs: (1) A⊥Ya,m|W and (2) M⊥Ya,m|W,A (Pearl 2001). The positivity assumptions are: P(M=m|A=a,W)<0 a.e. and P(A=a|W)<0 a.e. Two additional ignorability assumptions are required to identify NDE/NIEs: (3) A⊥Ma|W and (4) Ma∗⊥Ya,m|W (Pearl 2001). This last assumption states that, conditional on W, knowledge of M in the absence of treatment A provides no information of the effect of A on Y (Petersen, Sinisi, and van der Laan 2006). This assumption is violated when there is a confounder of the M−Y relationship that is affected by A (i.e., an intermediate confounder) (Avin, Shpitser, and Pearl 2005; Petersen, Sinisi, and van der Laan 2006; VanderWeele and Tchetgen Tchetgen 2017). This assumption is also problematic because it involves independence of counterfactuals under separate worlds (a and a∗) which can never simultaneously exist.

This last assumption that there is no confounder of the mediator-outcome relationship affected by prior treatment is especially concerning for epidemiology studies where longitudinal cohort data may reflect a data structure in which a treatment affects an individual characteristic measured at follow-up that in turn affects both a mediating variable and the outcome variable (see (Bild, D. E., Bluemke, D. A., Burke, G. L., Detrano, R., Roux, A. V. D., Folsom, A. R., Greenland, P., Jacobs Jr., D. R., Kronmal, R., Liu, K., et al 2002; Eaton and Kessler 2012; Phair et al. 1992) for some examples). It is also problematic for mediation analyses involving instrumental variables such as randomized encouragement-design interventions where an instrument, A, encourages treatment uptake, Z, which then may influence Y potentially through M. Such a design is present in the MTO experiment that we will use as an illustrative example. Randomization to receive a housing voucher (A) was the instrument that “encouraged” the treatment uptake, moving with the voucher out of public housing (Z, which we will call the intermediate confounder). In turn, Z may influence subsequent drug use among adolescent participants at follow-up (Y), possibly through a change the children’s school environment (M). In the illustrative example, our goal was to examine mediation of the effect of receiving a housing voucher (A) on subsequent drug use (Y) by changing school districts (M) in the MTO data.

There has been recent work to relax the assumption of no intermediate confounder, Ma∗⊥Ya,m|W, by using a stochastic intervention on M (Didelez, Dawid, and Geneletti 2006; van der Laan and Petersen 2008; VanderWeele and Tchetgen Tchetgen 2017; VanderWeele, Vansteelandt, and Robins 2014; Vansteelandt and Daniel 2017; Zheng and van der Laan 2017; 2012a). In this paper, we build on the approach described by VanderWeele and Tchetgen Tchetgen (2017) in which they defined the stochastic distribution on M as: gM|a,W or gM|a∗,W , where

(eq. 1 is an example of VanderWeele and Tchetgen Tchetgen (2017)’s work.) In other words, instead of formulating the individual counterfactual values of Ma or Ma∗, values are stochastically drawn from the distribution of M, conditional on covariates W but marginal over intermediate confounder Z, setting A=a or A=a∗, respectively. The corresponding estimands of interest are the SDE=E(Ya,gM|a∗,W)−E(Ya∗,gM|a∗,W), and SIE=E(Ya,gM|a,W)−E(Ya,gM|a∗,W).

Others have taken a similar approach. For example, Zheng and van der Laan (2017) formulate a stochastic intervention on M that is fully conditional on the past:

(note that per our SCM, P(M=1|Z,W)=P(M=1|Z,A,W), so in our case, gM|Z,a∗,W(Z,W)=P(M=1|Z,W).) The corresponding estimands are the stochastic direct and indirect effects fully conditional on the past: CSDE=E(Ya,gM|Z,a∗,W)−E(Ya∗,gM|Z,a∗,W), and CSIE=E(Ya,gM|Z,a,W)−E(Ya,gM|Z,a∗,W). However, Zheng and van der Laan (2017)’s formulation shown in eq. 2 is not useful for understanding mediation under the instrumental variable SCM we consider here, as there is no direct pathway from A to M. Because of the restriction on our statistical model that P(M|Z,A,W)=P(M|Z,W), gM|Z,a∗,W(Z,W)=gM|Z,a,W(Z,W), so CSIE’s under this model would equal 0. Thus, in this scenario, the NDE and CSDE are very different parameters. We note that it is also because of these restrictions on our statistical model stemming from the instrumental variable SCM that the sequential mediation analysis approach proposed by VanderWeele and Vansteelandt (2014) would also result in indirect effects equal to 0.

Because the CSIE and CSDE do not aid in understanding the role of M as a potential mediator in this scenario, we focus instead on VanderWeele and Tchetgen Tchetgen (2017)’s SDE and SIE that condition on W but marginalize over Z, thus completely blocking arrows into M (similar to an NDE and NIE). The SDE and SIE coincide with the NDE and NIE in the absence of intermediate confounders (VanderWeele and Tchetgen Tchetgen 2017).

3 Targeted minimum loss-based estimator

Our proposed estimator uses targeted minimum loss-based estimation (TMLE) (van der Laan and Rubin 2006), targeting the stochastic, counterfactual outcomes that comprise the SDE and SIE. To our knowledge, it is the first such estimator appropriate for instrumental variable scenarios. TMLE is a substitution estimation method that solves the EIC estimating equation. Its robustness properties differ for the fixed and data-dependent parameters. For the data-dependent SDE and SIE, if either the Y model is correct or the A and M models given the past are correct, then one obtains a consistent estimator of the parameter. Robustness to misspecification of the treatment model is relevant under an SCM with nonrandom treatment; we discuss the generalization of our proposed estimator to such an SCM in the Appendix. Note that g^M|a∗,W for the stochastic intervention is not the same as the conditional distribution of M given the past, so the first could be inconsistent while the latter is consistent. For the fixed SDE and SIE, we also need to assume consistent estimation of gM|a∗,W, since it does not target gM|a∗,W (and is therefore not a full TMLE for the fixed parameters).

The estimator integrates two previously developed TMLEs: one for stochastic interventions (Muñoz and van der Laan 2012) and one for multiple time-point interventions (van der Laan and Gruber 2012), which is built on the iterative/recursive g-computation approach (Bang and Robins 2005). This TMLE is not efficient under the SCM considered, because of the restriction on our statistical model that P(Y|M,Z,A,W)=P(Y|M,Z,W). However, it is still a consistent estimator if that restriction on our model does not hold (i.e., P(Y|M,Z,A,W)≠P(Y|M,Z,W)), because the targeting step adds dependence on A.

The TMLE is constructed using the sequential regressions described in the above section with an additional targeting step after each regression. The TMLE solves the EIC for the target parameter that treats gM|a∗,W as given. A similar EIC has been described previously (Bang and Robins 2005; van der Laan and Gruber 2012). The EIC for the parameter Ψ(P)(a,gM|a∗,W) for a given gM|a∗,W is given by:

Substitution of gM|a∗,W=g^M|a∗,W yields the EIC used for the data-dependent parameter Ψ(P)(a,g^M|a∗,W). The EIC for the parameter Ψ(P)(a,gM|a∗,W) in which the stochastic intervention equals the unknown gM|a∗,W is an area of future work.

We now describe how to compute the TMLE. In doing so, we use parametric model/regression language for simplicity but data-adaptive estimation approaches that incorporate machine learning

(e.g.,van der Laan, Polley, and Hubbard 2007) may be substituted and may be preferable (we use such a data-adaptive approach in the illustrative example analysis). We note that survey or censoring weights could be incorporated into this estimator as described previously (Rudolph et al. 2014). We use notation reflecting estimation of the data-dependent parameters, but note that estimation of the fixed parameters would be identical–-in the fixed parameter case, the notation would refer to gM|a∗,W instead of g^M|a∗,W.

First, one estimates g^M|a∗,W(W), which is the estimate of gM|a∗,W(W), defined in eq. 1, using observed data. Consider a binary Z. We estimate gZ|a∗,W(W)=P(Z=1|A=a∗,W). We then estimate gM|z,W(W)=P(M=1|Z=z,W) for z∈{0,1}. We use these quantities to calculate g^M|a∗,W=g^M|z=1,Wg^Z|a∗,W+g^M|z=0,W(1−g^Z|a∗,W). We can obtain g^Z|a∗,W(W) from a logistic regression of Z on A,W setting A=a∗, and g^M|z,W(W) from a logistic regression of M on Z,W, setting Z={0,1}. We will then use this stochastic intervention in the TMLE, whose implementation is described as follows.

1. Let Q¯Y,n(M,Z,W) be an estimate of Q¯Y(M,Z,W)≡E(Y|M,Z,W). To obtain Q¯Y,n(M,Z,W), predict values of Y from a regression of Y on M,Z,W.

2. Estimate the weights to be used for the initial targeting step:h1(a)=I(A=a){I(M=1)g^M|a∗,W+I(M=0)(1−g^M|a∗,W)}P(A=a){I(M=1)gM|Z,W+I(M=0)(1−gM|Z,W)}, where estimates of gM|Z,W are predicted probabilities from a logistic regression of M=m on Z and W. Let h^1,n(a) denote the estimate of h1(a).

3. Target the estimate of Q¯Y,n(M,Z,W) by considering a univariate parametric submodel {Q¯Y,n(M,Z,W)(ϵ):ϵ} defined as: logit(Q¯Y,n(ϵ)(M,Z,W))=logit(Q¯Y,n(M,Z,W))+ϵ. Let ϵn be the MLE fit of ϵ. We obtain ϵn by setting ϵ as the intercept of a weighted logistic regression model of Y with logit(Q¯Y,n(M,Z,W)) as an offset and weights h^1,n(a). (Note that this is just one possible TMLE.)The update is given by Q¯Y,n∗(M,Z,W)=Q¯Y,n(ϵn)(M,Z,W).Y can be bounded to the [0,1] scale as previously recommended (Gruber and van der Laan 2010).

4. Let Q¯M,ng(Z,W) be an estimate of Q¯Mg(Z,W). To obtain Q¯M,ng(Z,W), we integrate out M to from Q¯Y,n∗(M,Z,W). First, we estimate Q¯Y,n∗(M,Z,W) setting m=1 and m=0, giving Q¯Y∗(m=1,z,w) and Q¯Y∗(m=0,z,w). Then, multiply these predicted values by their probabilities under g^M|a∗,W(W) (for a∈{a,a∗}), and add them together (i.e., Q¯M,ng^(Z,W)=Q^Y∗(m=1,z,w)g^M|a∗,W+Q^Y∗(m=0,z,w)(1−g^M|a∗,W)).

5. We now fit a regression of Q¯M,ng^,∗(Z,W) on W among those with A=a. We call the predicted values from this regression Q¯Z,na(W). The empirical mean of these predicted values is the TMLE estimate of Ψ(P)(a,g^M|a∗,W).

6. Repeat the above steps for each of the interventions. For example, for binary A, we would execute these steps a total of three times to estimate: 1) Ψ(P)(1,g^M|1,W), 2) Ψ(P)(1,g^M|0W), and 3) Ψ(P)(0,g^M|0,W).

7. The SDE can then be obtained by substituting estimates of parameters Ψ(P)(a,g^M|a∗,W)−Ψ(P)(a∗,g^M|a∗,W) and the SIE can be obtained by substituting estimates of parameters Ψ(P)(a,g^M|a,W)−Ψ(P)(a,g^M|a∗,W).

8. For the fixed parameters, the variance can be estimated using the bootstrap. For the data-dependent parameters, the variance of each estimate from Step 6 can be estimated as the sample variance of the EIC (defined above, substituting in the targeted fits Q¯Y,n∗(M,Z,W) and Q¯Z,na,∗(W)) divided by n. First, we estimate the EIC for each component of the data-dependent SDE/SIE, which we call EICΨ(P)(a,g^M|a∗,W). Then we estimate the EIC for the estimand of interest by subtracting the EICs corresponding to the components of the estimand. For example EICSDE=EICΨ(P)(a,g^M|a∗,W)−EICΨ(P)(a∗,g^M|a∗,W). The sample variance of this EIC divided by n is the influence curve (IC)-based variance of the data-dependent estimator.

4 Simulation

5 Empirical illustration

5.2 Results

Figure 2 shows the data-dependent SDE and SIE estimates by type of estimator (TMLE, IPTW, and EE) for boys in the Boston MTO site (N=228). SDE and SIE estimates are similar across estimators. We find no evidence that change in school district mediated the effect of being randomized to the voucher group on marijuana use, with null SIE estimates (TMLE risk difference: − 0.003, 95% CI: − 0.032, 0.026). The direct effect of randomization to the housing voucher group on marijuana use suggests that boys who were randomized to this group were 9% more likely to use marijuana than boys in the control group, though this difference is not statistically significant (TMLE risk difference: 0.090, 95% CI: -0.065-0.245).

Figure 2:

Mediated effect estimates and 95% confidence intervals using interim follow-up data from adolescent boys in the Boston site of the Moving to Opportunity experiment. The data-dependent SDE is interpreted as the direct effect of being randomized to receive a housing voucher on risk of marijuana use that is not mediated through a change in school district. The data-dependent SIE is interpreted as the effect of being randomized to receive a housing voucher on marijuana use that is mediated by changing school districts.

6 Discussion

We proposed robust targeted minimum loss-based estimators to estimate fixed and data-dependent stochastic direct and indirect effects that are the first to naturally accommodate instrumental variable scenarios. These estimators build on previous work identifying and estimating the SDE and SIE (VanderWeele and Tchetgen Tchetgen 2017). The SDE and SIE have the appealing properties of (1) relaxing the assumption of no intermediate confounder affected by prior exposure, and (2) utility in studying mediation in the context of instrumental variables that adhere to the exclusion restriction assumption (a common assumption of instrumental variables which states that there is no direct effect between A and Y or between A and M (Angrist, Imbens, and Rubin 1996)) due to completely blocking arrows into the mediator by marginalizing over the intermediate confounder, Z. Given the restrictions that this assumption places on the statistical model, several alternative estimands are not appropriate for understanding mediation in this context as the indirect effect would always equal zero (e.g., Tchetgen Tchetgen 2013; VanderWeele and Vansteelandt 2014; Zheng and van der Laan 2017).

Inference for the fixed SDE and SIE can be obtained from bootstrapping, using parametric models for nuisance parameters. Inference for the data-dependent SDE and SIE can be obtained from the data-dependent EIC that assumes known g^M|a∗,W estimated from the data, and is appropriate for integrating machine learning in modeling nuisance parameters. The ability to incorporate machine learning is a significant strength in this case; if using the parametric alternative, multiple models would need to be correctly specified (VanderWeele and Tchetgen Tchetgen 2017). IC-based variance is possible in estimating the data-dependent SDE and SIE, because the data-dependent EIC has a form that is solvable using existing statistical tools; in contrast, the EIC for the fixed parameters is more complex and is not solvable with current statistical tools.

Our proposed estimator for the fixed and data-dependent parameters is simple to implement in standard statistical software, and we provide R code to lower implementation barriers. Another advantage of our TMLE estimator, which is shared with other estimating equation approaches, is that it is robust to some model misspecification. In estimating the data-dependent SDE and SIE, one could obtain a consistent estimate as long as either the Y model or the A and M models given the past were correctly specified. Obtaining a consistent estimate of the fixed SDE and SIE would also require consistent estimation of gM|a∗,W. In addition, our proposed estimation strategy is less sensitive to positivity violations than weighting-based approaches. First, TMLE is usually less sensitive to these violations than weighting estimators, due in part to it being a substitution estimator, which means that its estimates lie within the global constraints of the statistical model. This is in contrast to alternative estimating equation approaches, which may result in estimates that lie outside the parameter space. Second, we formulate our TMLE such that the targeting is done as a weighted regression, which may smooth highly variable weights (Stitelman, De Gruttola, and van der Laan 2012). In addition, moving the targeting into the weights improves computation time (Stitelman, De Gruttola, and van der Laan 2012).

However, there are also limitations to the proposed approach. We have currently only implemented it for a binary A and M, though extensions to multinomial or continuous versions of those variables are possible (Rosenblum and van der Laan 2010; Diaz and Rosenblum 2015). Extending the estimator to allow for a high-dimensional M is less straightforward, though it is of interest and an area for future work as allowing for high-dimensional M is a strength of other mediation approaches (Tchetgen Tchetgen 2013; Zheng and van der Laan 2017). We also plan to focus future work on developing a full TMLE for the fixed SDE and SIE parameters.