Direct, indirect, and interaction effects based on principal stratification with a binary mediator

1 Introduction

Given a binary treatment D , a binary mediator M , and an outcome/response variable Y , the causal chain of interest in mediation analysis is

The total effect of D on Y consists of the direct effect of D on Y and the indirect effect through M . This is an important issue in many disciplines of science, as reviewed by MacKinnon et al. [1], Pearl [2], Imai et al. [3], TenHave and Joffe [4], Preacher [5], VanderWeele [6,7], and Nguyen et al. [8], among others.

The total effect of D on Y can be found in various ways, such as matching, regression adjustment and inverse probability weighting (see, e.g., Lee and Lee [9] and Choi and Lee [10] for reviews on treatment effect estimators). Traditionally, the total effect has been decomposed into sub-effects with linear structural-form (SF) models for M as a function of D and Y as a function of ( D , M ) , e.g., with some α and β parameters,

where ( ε , U ) are the error terms with E ( ε ∣ D ) = 0 and E ( U ∣ D , M ) = 0 , α d is the effect of D on M , β d is the direct effect of D on Y , β m is the direct effect of M on Y , and β d m is the interaction effect of D M on Y . Here, β m α d is the “indirect effect” of D on Y through M , which is the “effect of M on Y ” times the “effect of D on M .” The standard Baron-Kenny approach (Baron and Kenny [11]) does not include the interaction part β d m D M , which turns out to be a major complicating term for total effect decomposition.

Linear SFs such as (1.1) are intuitive and easy to understand, but they are subject to specification errors and not grounded in the modern counterfactual causal framework. Once we adopt nonparametric approaches to avoid SF misspecifications and use potential versions of ( M , Y ) , total effect decomposition is no longer straightforward as with (1.1). Note that when exogenous covariates X need to be controlled, an easy way to do it generalizing (1.1) is

where α ( X ) s and β ( X ) s are the functions of X . Specifying them as linear, ordinary least-squares estimator (OLS) can be applied to (1.2). This approach turns out to be the main estimator proposed in this article, although (1.2) will be derived formally using the potential versions of M and Y .

Consider two potential versions M d , d = 0 , 1 , of M corresponding to D = 0 , 1 , and the four potential responses Y d m of Y for D = 0 , 1 and M = 0 , 1 . Also, define the potential responses Y d , d = 0 , 1 , corresponding to D = 0 , 1 “when M is allowed to take its natural course given D = d ”:

Then, the mean/average total effect is

In the literature, two-way decompositions of the total effect appeared first, as can be seen in Pearl [12] and Robins [13], among others:

The terms in { ⋅ } are called the “natural direct effects” of D on Y , and the terms in (⋅) are the “natural indirect effects.” In contrast, the “controlled direct effect” with M = m is E ( Y 1 , m − Y 0 , m ) , but there is no definition of “controlled indirect effect” that is generally agreed on. One glaring problem with (1.4) is that the decomposition is not unique: which one to take between (a) and (b)? Another problem is that the presence of the interaction effect is not clear in (1.4).

VanderWeele [14] proposed a three-way decomposition separating the aforementioned effects, and VanderWeele [15] proposed a four-way decomposition, which includes the existing two- and three-way decompositions in the literature as special cases by merging different terms in the four-way decomposition.

In this article, we define four subject types with ( M 0 , M 1 ) :

When D is endogenous with a binary instrumental variable δ for D , Imbens and Angrist [16] and Angrist et al. [17] classified the subjects using the potential treatments ( D 0 , D 1 ) corresponding to δ = 0 , 1 , and the classification is analogous to (1.5), e.g., compliers are those with ( D 0 = 0 , D 1 = 1 ) . Frangakis and Rubin [18] called the classification based on ( D 0 , D 1 ) “principal stratification,” and doing analogously, we call (1.5) “mediative principal stratification” based on ( M 0 , M 1 ) . We address only binary M in this article; allowing non-binary M is not straightforward and thus left for a future research.

With many decompositions of the total effect in the literature, the question is which one to use. There seems no single best answer to this question, but once the mediator types are taken into account, the following three-way decomposition seems well suited:

“M3M” stands for “Main 3-way Mediator-based decomposition.” M3M appeared in VanderWeele [15] with different notation, who also referred to VanderWeele and Tchetgen Tchetgen [19]. The following explains the three terms/effects in M3M.

First, the direct effect Y 10 − Y 00 occurs when D changes while the presence of M is nullified by m = 0 in Y d m . Second, the indirect effect through M is the product of the “direct effect Y 01 − Y 00 of M on Y (with D nullified by d = 0 in Y d m )” times the “effect M 1 − M 0 of D on M ,” which occurs only to CP and DF (with the opposite signs) as only they have M 1 ≠ M 0 . Third, the interaction effect of D M is the “net effect of D M ,” which is the “gross effect of D M ” minus the “direct/partial” effects of D and M (Choi and Lee [20]):

This effect occurs only to CP and AT because only they have D M = D M 1 = 1 , which is why M 1 ( = 1 only for CP and AT) is attached to Δ Y ± in M3M. Note that the interaction effect can be viewed either as the direct effect of D moderated by the level of M or as the direct effect of M moderated by the level of D .

Our approach based on “mediative principal stratification” looks at which effects are associated with which type, so that each type’s contributions to the total effect suggest the appropriate decomposition of the total effect, as was seen just above. Our approach will further show that, despite the apparent “symmetry” in (1.4)(a) and (b), (1.4)(b) is better than (1.4)(a) because (1.4)(a) does not identify any effect in M3M, whereas (1.4)(b) identifies the indirect effect and the sum of the direct and interaction effects.

In the remainder of this article, Section 2 shows the details of M3M. Section 3 addresses identifying all sub-effects in M3M, and maps out our estimation strategy. Section 4 presents the effect estimators. Sections 5 and 6 provide simulation and empirical studies. Finally, Section 6 concludes this article. Appendix contains most proofs.

2 Main three-way mediator-based decomposition

In this section, first, we formally put forth M3M as a preferred decomposition and explain why. Second, we compare M3M to other decompositions in the literature. Third, we illustrate various decompositions, using simple SFs as in (1.1).

3 Identification and estimation strategy

Our identification conditions for M3M with X are (“ ∐ ” stands for independence):

Conditions C(a) and C(b) are the “ignorability” of confounders in the treatment-mediator, treatment-outcome, and mediator-outcome relationships. C(a) and C(b) operate in two stages: as D precedes M , which, in turn, precedes Y , the first stage is D being independent of all potential future variables given X , and the second stage is ( M 0 , M 1 ) being independent of all potential future versions of Y given ( D , X ) . C(c) is a support-overlap condition for D ∣ X and M ∣ ( D , X ) ; for example, C(c) implies

Slightly different conditions from C(a) and C(b) appeared in Imai et al. [3]:

where the joint distributions of ( M 0 , M 1 ) and ( Y 00 , Y 01 , Y 10 , Y 11 ) do not appear, differently from C(a) and C(b). Also, Petersen et al. [21] assumed

Using marginal independence instead of joint independence, we can relax C(a) and C(b), but we continue to assume C(a) and C(b) for simplicity. In the following, we present “causal reduced forms (CRF’s)” for M and Y , which form the basis for identification.

Regarding (3.1), it holds for any form of M as long as M 1 − M 0 makes sense. However, since M3M is based on binary M , we continue to assume binary M . There are two “cells” D = 0 , 1 , and ψ 1 ( X ) and ψ d ( X ) are nonparametrically identified using

This point can be understood by considering a nonparametric estimation of M on X for the D = 0 , 1 groups, separately.

Analogously to (3.1), (3.2) holds also for any form of Y , as long as Y 10 − Y 00 , Y 01 − Y 00 , and Δ Y ± make sense; For example, Y can be binary, continuous, etc. There are four cells formed by D = 0 , 1 and M = 0 , 1 , and analogously to the identification of { ψ 1 ( X ) , ψ d ( X ) } , { μ 1 ( X ) , μ d ( X ) , μ m ( X ) , μ d m ( X ) } are nonparametrically identified using the X -conditional means on the cells.

A model is a SF, if it is a data-generating process with parameters of interest governing the behaviors of subjects. A model is a RF, if it is derived from some SF’s; e.g., a SF for Y has D and M on the right-hand side, and substituting out the SF for M yields the RF for Y with only D on the right-hand side. The parameters in a RF are not of direct interest, as they are functions of the underlying SF parameters. CRF falls between SF and RF, as it is a derived form or RF but with causal parameters of interest, such as μ d ( X ) , μ m ( X ) , and μ d m ( X ) in (3.2). CRF’s similar to (3.1) and (3.2) appeared in Lee [22,23], Mao and Li [24], Choi et al. [25], and Lee et al. [26].

To understand the effect decomposition better, substitute (3.1) into (3.2):

The slope of D is the “ X -conditioned total effect” consisting of direct ( μ d ( X ) ), indirect ( μ m ( X ) ψ d ( X ) ) and interaction ( μ d m ( X ) { ψ 1 ( X ) + ψ d ( X ) } ) effects. Compare these to the constant-effect versions in (2.9): β d + β m α d + β d m { α 1 + E ( X ′ ) α x + α d } .

If desired, { ψ 1 ( X ) , ψ d ( X ) , μ 1 ( X ) , μ d ( X ) , μ m ( X ) , μ d m ( X ) } in (3.1) and (3.2) can be estimated nonparametrically. However, more practical would be specifying those functions as, e.g., linear functions of X , and then, apply OLS to the resulting M and Y models. The details of this approach will be seen in the next section.

The OLS just mentioned is reminiscent of the traditional approach with (1.1), which raises the question: how much is the aforementioned OLS different from the traditional OLS to (1.1)? The answer is that there are three critical differences.

First, whereas the SF’s in (1.1), i.e., the data-generating processes, hold only for continuous M and Y in general, the aforementioned CRF’s for M and Y hold more generally for any forms of M and Y .

Second, (1.1) can be generalized to account for effect heterogeneity as in (1.2). However, proceeding with the SF’s in (1.2) makes it unclear what kind of direct, indirect, and interaction effects are actually estimated: e.g., are they (2.3) or M3M?

Third, if one starts off with (1.2), then it is not clear whether α d ( X ) , β d ( X ) , β m ( X ) , and β d m ( X ) in (1.2) should be restricted or not. In contrast, the definitions of the unknown functions of X in (3.1) and (3.2) show how they might be restricted. For example, since ψ 1 ( X ) ≡ E ( M 0 ∣ X ) , for binary M as in our setup, ψ 1 ( X ) cannot be specified just as X ′ α for a parameter α ; rather, ψ 1 ( X ) = Φ ( X ′ α ) is more suitable for a distribution function Φ ( · ) . Another example is μ d ( X ) ≡ E ( Y 10 − Y 00 ∣ X ) : if Y is binary, E ( Y 10 − Y 00 ∣ X ) should be bounded by [ − 1 , 1 ] , which can be accommodated by a smooth function with the range on [ − 1 , 1 ] , such as the arctan function or Φ ( X ′ α 10 ) − Φ ( X ′ α 00 ) , where Φ ( X ′ α 10 ) is for E ( Y 10 ∣ X ) and Φ ( X ′ α 00 ) is for E ( Y 00 ∣ X ) . If these nonlinear functions are used, then the aforementioned OLS should be replaced by a nonlinear least-squares estimator, which will be further discussed in the next section.

4 Effect estimators using OLS and generalized method of moment (GMM)

To estimate the effects in M3M, we linearly approximate all ψ ( X ) and μ ( X ) terms in the CRF’s (3.1) and (3.2), and then apply OLS to (3.1) and (3.2). This is summarized in Proposition 3, which is followed by discussions on nonlinear estimation when some ψ ( X ) and μ ( X ) terms are nonlinearly specified.

Let X α , X 1 , X d , X m , X d m consist of elements of X and their functions, with X j of dimension k j × 1 , j = α , 1 , d , m , d m . Linearly approximate the functions of X in (3.1):

α 1 ′ X α is for ψ 1 ( X ) in (3.1) and α d ′ X α is for ψ d ( X ) . For simplicity, we use the same X α in α 1 ′ X α and α d ′ X α , but we can certainly use different covariates, if desired.

Doing analogously for (3.2), we have

where β 1 ′ X 1 , β d ′ X d , β m ′ X m , and β d m ′ X d m are for μ 1 ( X ) , μ d ( X ) , μ m ( X ) , and μ d m ( X ) .

One might use different covariates X α , X 1 , X d , X m , and X d m as in (4.1) and (4.2), but this approach would be following the conventional SF view as in (1.2). Rather, since our CRF’s are not SF’s and the X -conditional causal parameters are of RF type (e.g., ψ d ( X ) ≡ E ( M 1 − M 0 ∣ X ) is a RF function), it is better to control for all X to ensure C(a) to C(c), i.e., setting X = X α = X 1 = X d = X m = X d m in (4.1) and (4.2) would be fine, as will be done in our simulation and empirical studies.