Robust and Flexible Estimation of Stochastic Mediation Effects: A Proposed Method and Example in a Randomized Trial Setting
Kara E. Rudolph Oleg Sofrygin Wenjing Zheng Mark J. van der Laan
stat.AP stat.ME 62P25 G.3

Dynamic band structure tuning of graphene moiré superlattices with pressure
    Cory R. Dean^1
    December 30, 2023
===========================================================================






  §.§.§ Abstract

Background: Causal mediation analysis can improve understanding of the mechanisms underlying epidemiologic associations. However, the utility of natural direct and indirect effect estimation has been limited by the assumption of no confounder of the mediator-outcome relationship that is affected by prior exposure—an assumption frequently violated in practice. 


Methods: We build on recent work that identified alternative estimands that do not require this assumption and propose a flexible and double robust semiparametric targeted minimum loss-based estimator for data-dependent stochastic direct and indirect effects. The proposed method 
 treats the intermediate confounder affected by prior exposure as a time-varying confounder and intervenes stochastically on the mediator using a distribution which conditions on baseline covariates and marginalizes over the intermediate confounder. 
  In addition, we assume the stochastic intervention is given, conditional on observed data, which results in a simpler estimator and weaker identification assumptions.
 

Results: We demonstrate the estimator's finite sample and robustness properties in a simple simulation study. We apply the method to an example from the Moving to Opportunity experiment. In this application, randomization to receive a housing voucher is the treatment/instrument that influenced moving to a low-poverty neighborhood, which is the intermediate confounder. We estimate the data-dependent stochastic direct effect of randomization to the voucher group on adolescent marijuana use not mediated by change in school district and the stochastic indirect effect mediated by change in school district. We find no evidence of mediation. 


Conclusions: Our estimator is easy to implement in standard statistical software, and we provide annotated R code to further lower implementation barriers. 




Keywords: mediation; direct effect; indirect effect; double robust; targeted minimum loss-based estimation; targeted maximum likelihood estimation; data-dependent 











§ 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 <cit.>. 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 <cit.>. Ultimately, being randomized to receive a voucher affected subsequent adolescent drug use <cit.>. 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 <cit.> (also called mediation analysis using the counterfactual framework <cit.>) shares similar goals with the standard mediation approaches, e.g., structural equation modeling and the widely used Baron and Kenny “product method” approach <cit.>. 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 
 <cit.>, 
 2) allow for modeling nonlinear relationships between mediators and outcomes <cit.>, 
 and 3) allow for incorporation of data-adaptive machine learning methods and double robust estimation <cit.>. 

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 <cit.>, a time-varying confounder affected by prior exposure <cit.>, or time-dependent confounding by an intermediate covariate <cit.>). 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 <cit.>. 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.



§ NOTATION AND CAUSAL MEDIATION ESTIMANDS


Let observed data: O=(W, A, Z, M, Y) with n i.i.d. copies O_1,...,O_n ∼ P_0, 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: (U_W, U_A, U_Z, U_M, U_Y). This SCM is represented in Figure <ref> and the following causal models: W=f(U_W), A=f(U_A), Z=f(A, W, U_Z), M=f(Z,W,U_M), and Y=f(M,Z,W,U_Y). Note that this SCM (including that U_Y and U_M 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).




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(Y_a, m - Y_a*, m), where Y_a,m, Y_a^*,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(Y_a, M_a^* - Y_a*, M_a^*), where Y_a,M_a^*, Y_a^*,M_a^* is the counterfactual outcome setting A equal to a or a^* but this time setting M_a^* 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(Y_a, M_a - Y_a, M_a*). Natural direct and indirect effects are frequently used in epidemiology and have the appealing property of adding to the total effect <cit.>.

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 ⊥ Y_a,m | W and  2) M ⊥ Y_a,m | W, A <cit.>. 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 ⊥ M_a | W and 4) M_a*⊥ Y_a,m | W <cit.>. 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 <cit.>. This assumption is violated when there is a confounder of the M-Y relationship that is affected by A (i.e., an intermediate confounder) <cit.>. 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 <cit.> 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, M_a*⊥ Y_a,m | W, by using a stochastic intervention on M <cit.>. 
In this paper, we build on the approach described by 
 <cit.> in which they defined the stochastic distribution on M as: g_M|a,W or g_M|a^*,W, where 
    g_M|A,W(m,a^*,W) ≡ g_M|a^*,W(W) = ∑_z=0^1 P(M=1|Z=z,W)P(Z=z | A=a^*, W).
 (Equation 1 is an example of <cit.>'s work.) 

 In other words, instead of formulating the individual counterfactual values of M_a or M_a^*, 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(Y_a,g_M|a^*,W) - E(Y_a*,g_M|a^*,W), and SIE = E(Y_a,g_M|a,W) - E(Y_a,g_M|a^*,W). 
 
 Others have taken a similar approach. For example, <cit.> formulate a stochastic intervention on M that is fully conditional on the past: 
    g_M|Z,A,W(m,Z,a^*,W) ≡ g_M|Z,a^*,W (Z,W) = P(M=1|Z, A=a^*, W)
 (note that per our SCM, P(M=1 | Z,W) = P(M=1 | Z,A,W), so in our case, g_M|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(Y_a,g_M|Z,a^*,W) - E(Y_a*,g_M|Z,a^*,W), and CSIE = E(Y_a,g_M|Z,a,W) - E(Y_a,g_M|Z,a^*,W). However, <cit.>'s formulation shown in Equation <ref> is not useful for understanding mediation under the instrumental variable SCM we consider here, as there is no direct pathway from A to M or from A to Y. Because of the restriction on our statistical model that P(M | Z, A, W) = P(M | Z, W), g_M|Z,a^*,W (Z,W) = g_M|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 <cit.> 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 <cit.>'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 <cit.>. 



 §.§ SDE and SIE Estimands and Identification
 Our proposed estimator can be used to estimate two versions of the SDE and SIE: 1) fixed parameters that assume an unknown, true g_M|a^*,W; and 2)  data-dependent parameters that assume known g_M|a^*,W, estimated from the observed data, which we call ĝ_M|a^*,W. Researchers may have defensible reasons for choosing one version over the other, which we explain further below. The fixed SDE and SIE can be identified from the observed data distribution using the g-computation formula as discussed by <cit.>, assuming the sequential randomization assumption on intervention nodes A and M: 1) A ⊥ Y_a,m | W, 2) M ⊥ Y_a,m | W, A=a, Z, and 3) A ⊥ M_a | W, for a particular a and g_M|a^*,W. The data-dependent SDE And SIE can be identified similarly but need only the first two assumptions, 1) A ⊥ Y_a,m | W and 2) M ⊥ Y_a,m | W, A=a, Z, because ĝ_M|a^*,W is assumed known. If any of the above identification assumptions are violated, then the statistical estimands will not converge to their true causal quantities.

The estimand E(Y_a,ĝ_M|a^*W) can be identified via sequential regression, which provides the framework for our proposed estimator that follows. For intervention (A=a, M=ĝ_M|a^*,W), we have Q̅_M^ĝ(Z,W)≡∫_m ∈ℳ∫_y ∈𝒴 p(y | m, Z, W) dμ_Y(y) ĝ_m|a^*,W dμ_M(m), 
  where we integrate out M under our stochastic intervention ĝ_M|a^*,W, and where M has support ℳ and Y has support 𝒴 and where dμ_Y(y) and dμ_M(m) are some dominating measures. This is accomplished by evaluating E(Y|M=m,Z=z,W) at each m and multiplying it by the probability that M=m under ĝ_M|a^*,W, summing over all m. We then integrate out Z and set A=a: Q̅_Z^a(W)≡∫_z ∈𝒵Q̅_M^ĝ(z,W|A=a,W) p(z | A=a, W) dμ_Z(z), where 𝒵 denotes the support of random variable Z and dμ_Z(z) is some dominating measure. Marginalizing over the distribution of W gives the statistical parameter: Ψ(P)(a,ĝ_M|a^*,W)= ∫_w ∈𝒲Q̅_Z^a(w) p(w) dμ_W(w), where 𝒲 denotes the support of random variable W and dμ_W(w) is some dominating measure.

In the next section, we propose a novel, robust substitution estimator that can be used to estimate both the fixed and parametric versions of the SDE and SIE. Inference for the fixed SDE and SIE can be obtained by using the bootstrapped variance, which requires parametric models for the nuisance parameters P(A) and P(M | Z,A,W). This is the same inference strategy as proposed by <cit.>. However, researchers may encounter scenarios for which fitting parametric models is unappealing and using machine learning approaches is preferred. The data-dependent SDE and SIE with inference based on the efficient influence curve (EIC) may be preferable in such scenarios. In contrast to the EIC for an assumed known g_M|a^*,W, the EIC an unknown g_M|a^*,W is complicated due to the dependence of the unknown marginal stochastic intervention 
 on the data distribution. Such an EIC would include an M component, the form of which would be more complex due to the distribution of M being marginalized over Z. No statistical tools for solving an EIC of that form currently exist. Solving the EIC for the the parameter Ψ(P)(a,g_M|a^*,W) for an unknown g_M|a^*,W is ongoing work.

 
  


§ TARGETED MINIMUM LOSS-BASED ESTIMATOR





Our proposed estimator 
 uses targeted minimum loss-based estimation (TMLE) <cit.>, 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 ĝ_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 g_M|a^*,W, since it does not target g_M|a^*,W (and is therefore not a full TMLE for the fixed parameters).  

The estimator integrates two previously developed TMLEs: one for stochastic interventions <cit.> and one for multiple time-point interventions <cit.>, which is built on the iterative/recursive g-computation approach <cit.>. 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 g_M|a^*,W as given. A similar EIC has been described previously <cit.>. The EIC for the parameter Ψ(P)(a, g_M|a^*,W) for a given g_M|a^*,W is given by:

    D^*(a, g_M|a^*,W)     = ∑_k=0^2 D_k^*(a, g_M|a^*,W),  where 
    
     
     
    D^*_0(a, g_M|a^*,W)     = Q̅^a_Z(W) - Ψ(P)(a, g_M|a^*,W)
    
        D^*_1(a, g_M|a^*,W)     = I(A=a)/P(A=a | W)(Q̅^g_M(Z,W) - Q̅^a_Z(W))
    
        D^*_2(a, g_M|a^*,W)     = I(A=a){I(M=1)g_M|a^*,W + I(M=0)(1-g_M|a^*,W) }/P(A=a | W){I(M=1)g_M|z,W + I(M=0)(1-g_M|z,W) }
       × (Y-Q̅_Y(M,Z,W)).

Substitution of g_M|a^*,W=ĝ_M|a^*,W yields the EIC used for the data-dependent parameter Ψ(P)(a, ĝ_M|a^*,W). The EIC for the parameter Ψ(P)(a, g_M|a^*,W) in which the stochastic intervention equals the unknown g_M|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 <cit.> 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 <cit.>. 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 g_M|a^*,W instead of ĝ_M|a^*,W.  

First, one estimates ĝ_M|a^*,W(W), which is the estimate of g_M|a^*,W(W), defined in Equation 1, using observed data. Consider a binary Z. We estimate g_Z|a^*,W(W)=P(Z=1 | A=a^*, W). We then estimate g_M|z,W(W)=P(M=1 | Z=z, W) for z ∈{0,1}. We use these quantities to calculate ĝ_M|a^*,W = ĝ_M|z=1,Wĝ_Z|a^*,W +  ĝ_M|z=0,W(1-ĝ_Z|a^*,W). We can obtain ĝ_Z|a^*,W(W) from a logistic regression of Z on A, W setting A=a^*, and ĝ_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.



  * 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. 
 
 



  * Estimate the weights to be used for the initial targeting step:
 h_1(a) = I(A=a){I(M=1)ĝ_M|a^*,W + I(M=0)(1-ĝ_M|a^*,W) }/P(A=a){I(M=1)g_M|Z,W + I(M=0)(1-g_M|Z,W) }, where estimates of g_M|Z,W are predicted probabilities from a logistic regression of M=m on Z and W. Let ĥ_1,n(a) denote the estimate of h_1(a). 

  * 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 ĥ_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 <cit.>.

  * Let Q̅^g_M,n(Z,W) be an estimate of Q̅^g_M(Z,W). To obtain Q̅^g_M,n(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 
 ĝ_M|a^*,W(W) (for a ∈{a, a^*}), and add them together (i.e., Q̅^ĝ_M,n(Z,W)= Q̂^*_Y(m=1, z, w)ĝ_M|a^*,W + Q̂^*_Y(m=0, z, w)(1-ĝ_M|a^*,W)).



  * 
 We now fit a regression of Q̅^ĝ,*_M,n(Z,W) on W among those with A=a. We call the predicted values from this regression Q̅^a_Z,n(W). The empirical mean of these predicted values is the TMLE estimate of Ψ(P)(a, ĝ_M|a^*,W).





  * 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,ĝ_M|1,W), 2) Ψ(P)(1,ĝ_M|0W), and 3) Ψ(P)(0,ĝ_M|0,W). 


  * The SDE can then be obtained by substituting estimates of parameters Ψ(P)(a,ĝ_M|a^*,W) - Ψ(P)(a^*,ĝ_M|a^*,W) and the SIE can be obtained by substituting estimates of parameters Ψ(P)(a,ĝ_M|a,W) - Ψ(P)(a,ĝ_M|a^*,W).



  * 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̅^a,*_Z,n(W)) divided by n. First, we estimate the EIC for each component of the data-dependent SDE/SIE, which we call EIC_Ψ(P)(a,ĝ_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 EIC_SDE = EIC_Ψ(P)(a,ĝ_M|a^*,W) - EIC_Ψ(P)(a^*,ĝ_M|a^*,W). The sample variance of this EIC divided by n is the influence curve (IC)-based variance of the data-dependent estimator. 




§ SIMULATION



 §.§ Data generating mechanism

We conduct a simulation study to examine finite sample performance of the TMLE estimators for the fixed SDE and SIE and data-dependent SDE and SIE from the data-generating mechanism (DGM) shown in Table <ref>. Under this DGM, the data-dependent SDE is: SDE=E(Y_1, ĝ_M|0,W) - E(Y_0, ĝ_M|0,W) and the SIE is: SIE=E(Y_1, ĝ_M|1,W) - E(Y_1, ĝ_M|0,W). The fixed versions are defined with respect to the unknown, true g_M|1,W and g_M|0,W. Table <ref> uses the same notation and SCM as in Section 2, with the addition of Δ, an indicator of selection into the sample (which corresponds to the MTO data used in the empirical illustration where one child from each family is selected to participate). 



We compare performance of the TMLE estimator to an inverse-probability weighted estimator (IPTW) and estimator that solves the EIC estimating equation (EE) but differs from TMLE in that it lacks the targeting steps and is not a plug-in estimator, so its estimates are not guaranteed to lie within the parameter space (which may be particularly relevant for small sample sizes). Variance for the fixed SDE and SIE parameters is calculated using 500 bootstrapped samples for each simulation iteration. Variance for the data-dependent SDE and SIE is calculated using the EIC. We show estimator performance in terms of absolute bias, percent bias, closeness to the efficiency bound (mean estimator standard error (SE) × the square root of the number of observations), 95% confidence interval (CI) coverage, and mean squared error (MSE) across 1,000 simulations for sample sizes of N=5,000, N=500, and N=100. In addition, we consider 1) correct specification of all models, and 2) misspecification of the Y model that included a term for Z only. 



 §.§ Performance

Table <ref> gives simulation results under correct model specification for fixed SDE and SIE using bootstrap-based variance. Table <ref> gives simulation results under correct model specification for data-dependent SDE and SIE using IC-based variance. Both Tables <ref> and <ref> show that the TMLE, IPTW, and EE estimators are consistent when all models are correctly specified, showing biases of around 1% or less under large sample size (N=5,000) and slightly larger biases with smaller sample sizes. The 95% CIs for the TMLE and EE estimators result in similar coverage that is close to 95%, except for estimation of the SIE with a sample size of N=100, which has coverage closer to 90%.
Confidence intervals for the IPTW estimator for the fixed parameter are close to 95% but are conservative and close to 100% for the data-dependent parameter. As expected, IPTW is less efficient than TMLE or EE; the TMLE and EE estimators perform similarly and close to the efficiency bound for all sample sizes.





Table <ref> gives simulation results under misspecification of the outcome model that only includes a term for Z for fixed SDE and SIE using bootstrap-based variance. Thus, comparing results in Table <ref> to Table <ref> demonstrates robustness to misspecification of the outcome model. As all three of the estimators evaluated are theoretically robust to misspecification of this model, we would expect similar results between the two Tables, and we see that is indeed the case.










§ EMPIRICAL ILLUSTRATION



 §.§ Overview and set-up

We now apply our proposed estimator to MTO: a longitudinal, randomized trial that is described above. Because we wish to use machine learning for this empirical illustration, we will estimate the data-dependent SDE of being randomized to receive a housing voucher (A) on marijuana use (Y) not mediated by change in school district (M) and the data-dependent SIE mediated by M among adolescent boys in the Boston site in the presence of an intermediate confounder (Z), moving with the voucher out of public housing. 

We restrict to adolescents less than 18 years old who were present at interim follow-up, as those participants had school data and were eligible to be asked about marijuana use. We restrict to boys in the Boston site as previous work has shown important quantitative and qualitative differences in MTO's effects by sex <cit.> and by city <cit.>. We choose to present results from a restricted analysis instead of a stratified analysis, as our goal is to illustrate the proposed method. A more thorough mediation analysis considering all sexes and sites is the subject of a future paper. 

Marijuana use was self-reported by adolescents at the interim follow-up, which occurred 4-7 years after baseline, and is defined as ever versus never use. Change in school district is defined as the school at follow-up and school at randomization being in the same district. Numerous baseline characteristics included individual and family sociodemographics, motivation for participating in the study, neighborhood perceptions, school-related characteristics of the adolescent, and predictive interactions. 

We use machine learning to flexibly and data-adaptively model the following relationships: instrument to intermediate confounder, intermediate confounder to mediator, and mediator to outcome. Specifically, we use least absolute shrinkage and selection operator (lasso) <cit.> and choose the model that improves 10-fold cross-validation prediction error, while always including age and race/ethnicity and relevant A, Z, and M variables. 



 §.§ Results

Figure <ref> 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). 





§ 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 <cit.>. 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 <cit.>) 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 <cit.>. 
 
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 ĝ_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 <cit.>.  
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 g_M|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 <cit.>. In addition, moving the targeting into the weights 
 improves computation time <cit.>.

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 <cit.>. 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 <cit.>. We also plan to focus future work on developing a full TMLE for the fixed SDE and SIE parameters. 

DeGruyter








§ GENERALIZATIONS TO OTHER STRUCTURAL CAUSAL MODELS



 §.§ Nonrandom Treatment

Let observed data: O=(W, A, Z, M, Y) with n i.i.d. copies O_1,...,O_n ∼ P_0, 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. We assume that M and Y are not affected by A except through its effect on Z. We assume that A is affected by {W}, Z is affected by {A, W}, 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: (U_W, U_A, U_Z, U_M, U_Y). Note that this SCM (including that U_Y and U_M 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). We can factorize the likelihood for this SCM as follows: P(O) = P(Y | M, Z, W)P(M|Z,W)P(Z | A, W)P(A | W)P(W).

The data-dependent, stochastic mediation estimand E(Y_a,ĝ_M|a^*W) can be identified via sequential regression, which provides the framework for our proposed estimator that follows. For intervention (A=a, M=ĝ_M|a^*,W), we have 
 Q̅_M^ĝ(Z,A,W)≡ E_ĝ_M|a^*,W(E(Y |M,Z,A,W)|Z,A,W), where we integrate out M under our stochastic intervention ĝ_M|a^*,W. This is accomplished by evaluating the inner expectation at each m and multiplying it by the probability that M=m under ĝ_M|a^*,W, summing over all m. We then integrate out Z and set A=a: Q̅_Z^a(W)≡ E_Z(Q̅_M^ĝ(Z,A,W)|A=a,W). Taking the empirical mean gives the statistical parameter: Ψ(P)(a,ĝ_M|a^*,W)=E_W(Q̅_Z^a(W)|W). 

The EIC for the parameter Ψ(P)(a, ĝ_M|a^*,W) is given by
    D^*(a, ĝ_M|a^*,W)     = ∑_k=0^2 D_k^*(a, ĝ_M|a^*,W),  where 
    
    D^*_0(a, ĝ_M|a^*,W)     = Q̅^a_Z(W) - Ψ(P)(a, ĝ_M|a^*,W)
    
        D^*_1(a, ĝ_M|a^*,W)     = I(A=a)/P(A=a | W)(Q̅^ĝ_M(Z,W) - Q̅^a_Z(W))
    
        D^*_2(a, ĝ_M|a^*,W)     = I(A=a){I(M=1)ĝ_M|a^*,W + I(M=0)(1-ĝ_M|a^*,W) }/P(A=a | W){I(M=1)g_M|Z,W + I(M=0)(1-g_M|Z,W) }(Y-Q̅_Y(M,Z,W)).


We now describe how to compute the TMLE. The estimation of ĝ_M|a^*,W(W) does not differ from that described in the main text. 



  * 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. 



  * Estimate the weights to be used for the initial targeting step:
 h_1(a) = I(A=a){I(M=1)ĝ_M|a^*,W + I(M=0)(1-ĝ_M|a^*,W) }/P(A=a){I(M=1)g_M|Z,W + I(M=0)(1-g_M|Z,W) }, where ĝ_M|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 h_1(a). 

  * 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.

  * Let Q̅^ĝ_M,n(Z,W) be an estimate of Q̅^ĝ_M(Z,W). To obtain  Q̅^ĝ_M,n(Z,W), we integrate out M 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,n(m=1,Z,W) and Q̅^*_Y,n(m=0,Z,W). Then, multiply these predicted values by their probabilities under ĝ_M|a^*,W(W) (for a ∈{a,a^*}), and add them together (i.e., Q̅^ĝ_M,n(Z,W) = Q̅^*_Y,n(m=1,Z,W)*ĝ_M|a^*,W + Q̅^*_Y,n(m=0,Z,W)*(1-ĝ_M|a^*,W)). 

  * 
 We now fit a regression of Q̅^ĝ,*_M,n(Z,W) on W among those with A=a. We call the predicted values from this regression Q̅^a_Z,n(W). 

  * Complete a second targeting step: logit(Q̅^a_Z,n(ϵ)(W)) = logit(Q̅^a_Z,n(W)) + ϵ h_2,n(a), where h_2,n(a) is an estimate of h_2(a) = I(A=a)/P(A=a | W) and P(A=a | W) can be estimated from a logistic regression model of A=a on W. Let ϵ_n again be the MLE fit of ϵ, which can be obtained by fitting an intercept-only weighted logistic regression model of Q̅^ĝ_M,n(Z,W) with logit(Q̅^a_Z,n(W)) as an offset and weights h_2,n(a). (Alternatively, we could fit an unweighted logistic regression model of Q̅^ĝ_M,n(Z,W) with logit(Q̅^a_Z,n(W)) as an offset and h_2,n(a) as a covariate, where ϵ_n is the fitted coefficient on h_2,n(a).) The update is given by Q̅^a,*_Z,n(W) = Q̅^a,a^*_Z,n(ϵ_n)(A,W).

  * The TMLE of Ψ(P)(a, ĝ_M|a^*,W) is the empirical mean of Q̅^a,*_Z,n(W).

  * 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,ĝ_M|1,W), 2) Ψ(P)(1,ĝ_M|0W), and 3) Ψ(P)(0,ĝ_M|0,W). 


  * The SDE can then be obtained by substituting estimates of parameters Ψ(P)(a,ĝ_M|a^*,W) - Ψ(P)(a^*,ĝ_M|a^*,W) and the SIE can be obtained by substituting estimates of parameters Ψ(P)(a,ĝ_M|a,W) - Ψ(P)(a,ĝ_M|a^*,W).



  * The variance of each estimate 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̅^a,*_Z,n(W)) divided by n. First, we estimate the EIC for each component of the SDE/SIE, which we call EIC_Ψ(P)(a,ĝ_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 EIC_SDE = EIC_Ψ(P)(a,ĝ_M|a^*,W) - EIC_Ψ(P)(a^*,ĝ_M|a^*,W). The sample variance of this EIC divided by n is the influence curve-based variance of the estimator. 




 §.§ There exist direct effects between A and M and between A and Y

Let observed data: O=(W, A, Z, M, Y) with n i.i.d. copies O_1,...,O_n ∼ P_0, 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. We assume that A is exogenous, Z is affected by {A, W}, M is affected by {A, Z, W}, and that Y is affected by {A, M, Z, W}. We assume exogenous random errors: (U_W, U_A, U_Z, U_M, U_Y). We can factorize the likelihood for this SCM as follows: P(O) = P(Y | M, Z, A,W)P(M|Z,A,W)P(Z | A, W)P(A)P(W).

Under this SCM, our proposed TMLE is efficient. The sequential regression used to identify the data-dependent, stochastic mediation estimands does not change from the that given in the main text for this SCM. 

The EIC for the parameter Ψ(P)(a, ĝ_M|a^*,W) is given by
    D^*(a, ĝ_M|a^*,W)     = ∑_k=0^2 D_k^*(a, ĝ_M|a^*,W),  where 
    
     
     
    D^*_0(a, ĝ_M|a^*,W)     = Q̅^a_Z(W) - Ψ(P)(a, ĝ_M|a^*,W)
    
        D^*_1(a, ĝ_M|a^*,W)     = I(A=a)/P(A=a)(Q̅^ĝ_M(Z,A,W) - Q̅^a_Z(W))
    
        D^*_2(a, ĝ_M|a^*,W)     = I(A=a){I(M=1)ĝ_M|a^*,W + I(M=0)(1-ĝ_M|a^*,W) }/P(A=a){I(M=1)g_M|Z,A,W + I(M=0)(1-g_M|Z,A,W) }(Y-Q̅_Y(M,Z,A,W)).


We now describe how to compute the TMLE. First, one estimates ĝ_M|a^*,W(W) = ∑_z=0^1 P(M=1|Z=z,A=a^*,W)P(Z=z | A=a^*, W). Consider a binary Z. We first estimate g_Z|a^*,W(W)=P(Z=1 | A=a^*, W). We then estimate g_M|z,a^*,W(W)=P(M=1 | Z=z, A=a^*, W) for z ∈{0,1}. We use these quantities to calculate ĝ_M|a^*,W = (ĝ_M|z=1,a^*,W×ĝ_Z|a^*,W) +  (ĝ_M|z=0,a^*,W× (1-ĝ_Z|a^*,W)). We can obtain ĝ_Z|a^*,W(W) from a logistic regression of Z on A, W setting A=a^*, and ĝ_M|z,a^*,W(W) from a logistic regression of M on Z, A, W, setting Z={0,1} and A=a^*. We will then use this data-dependent stochastic intervention in the TMLE, whose implementation is described as follows.



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

  * Estimate the weights to be used for the initial targeting step:
 h_1(a) = I(A=a){I(M=1)ĝ_M|a^*,W + I(M=0)(1-ĝ_M|a^*,W) }/P(A=a){I(M=1)g_M|Z,A,W + I(M=0)(1-g_M|Z,A,W) }, where ĝ_M|Z,A,W are predicted probabilities from a logistic regression of M=m on Z, A, and W. Let h_1,n(a) denote the estimate of h_1(a). 

  * Target the estimate of Q̅_Y,n(M,Z,A,W) by considering a univariate parametric submodel {Q̅_Y,n(M,Z,A,W)(ϵ):ϵ} defined as:
logit(Q̅_Y,n (ϵ)(M,Z,A,W)) = logit(Q̅_Y,n(M,Z,A,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,A,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,A,W) =  Q̅_Y,n(ϵ_n)(M,Z,A,W). Y can be bounded to the [0,1] scale.

  * Let Q̅^ĝ_M,n(Z,A,W) be an estimate of Q̅^ĝ_M(Z,A,W). To obtain  Q̅^ĝ_M,n(Z,A,W), we integrate out M from Q̅^*_Y,n(M,Z,A,W). First, we estimate Q̅^*_Y,n(M,Z,A,W) setting m=1 and m=0, giving Q̅^*_Y,n(m=1,Z,A,W) and Q̅^*_Y,n(m=0,Z,A,W). Then, multiply these predicted values by their probabilities under ĝ_M|a^*,W(W) (for a ∈{a,a^*}), and add them together (i.e., Q̅^ĝ_M,n(Z,A,W) = Q̅^*_Y,n(m=1,Z,A,W)*ĝ_M|a^*,W + Q̅^*_Y,n(m=0,Z,A,W)*(1-ĝ_M|a^*,W)). 

  * 
 We now fit a regression of Q̅^ĝ,*_M,n(Z,A,W) on W among those with A=a. We call the predicted values from this regression Q̅^a_Z,n(W). The empirical mean of these predicted values is the TMLE estimate of Ψ(P)(a, ĝ_M|a^*,W).





  * 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,ĝ_M|1,W), 2) Ψ(P)(1,ĝ_M|0W), and 3) Ψ(P)(0,ĝ_M|0,W). 


  * The SDE can then be obtained by substituting estimates of parameters Ψ(P)(a,ĝ_M|a^*,W) - Ψ(P)(a^*,ĝ_M|a^*,W) and the SIE can be obtained by substituting estimates of parameters Ψ(P)(a,ĝ_M|a,W) - Ψ(P)(a,ĝ_M|a^*,W).



  * The variance of each estimate 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̅^a,*_Z,n(W)) divided by n. First, we estimate the EIC for each component of the SDE/SIE, which we call EIC_Ψ(P)(a,ĝ_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 EIC_SDE = EIC_Ψ(P)(a,ĝ_M|a^*,W) - EIC_Ψ(P)(a^*,ĝ_M|a^*,W). The sample variance of this EIC divided by n is the influence curve-based variance of the estimator. 




§ FUNCTION CODE

functioncode.R
exampcode.R


