45 Causal mediation analysis

Published

July 21, 2022

45.1 Traditional mediation analysis

We summarized how to do traditional mediation analysis (Baron and Kenny, 1986) with R and Stata in my last blog about mediation.

Traditionally mediation model can be represented in the following equations:

\[ Y = a A + b M + \epsilon_1 \tag{45.1}\]

\[ M = c A + \epsilon_2 \tag{45.2}\]

There are three main variables, A (treatment, or exposure), M (mediator) and Y (outcome). There are two pathways from A to Y, a direct effect, and an indirect effect through A -> M -> Y.

The traditional approach is to model these two equations directly, and get the estimates of direct and indirect effects.

45.2 Causal mediation analysis

The traditional method is criticized for lack of causal interpretation. The causal mediation analysis makes assumptions up front. If those assumptions are met, then you have a causal interpretation.

In a random experiment, usually A is randomly assigned, but M is not. If there is a confounder on the M to Y pathway, then we don’t have a causal effect, or say we cannot estimate causal effect without confounder bias.

Causal mediation analysis needs that too; it’s not like if you use a causal mediation package then you get the causal effect.

The other shortcoming of the traditional mediation method is that it assumes a homogeneous treatment indirect effect; that is, there is no $A$ - $M$ interaction. With causal mediation, this interaction can be allowed.

45.2.1 Controlled direct effect

Suppose we can set $A$ and $M$ at will to any $(a, m)$, then we have the potential outcome $Y(a,m)$.

Controlled direct effect (CDE) is defined as \[ CDE(m) = E(Y(1,m)-Y(0,m)) \tag{45.3}\] That is, setting $M$ to $m$, what is the effect of $A$ on $Y$?

45.2.1.1 Assumptions

What are the assumptions for CDE to be identified (estimable)?

Conditional treatment randomization. Suppose we observe confounders W, which can be a joint set of confounders for the A to Y pathway, or the M to Y pathway.

\[ Y(a,m) \perp A|W \tag{45.4}\]

This is the usual conditional independence (unconfoundedness, ignorability, etc.) assumption. This is saying the assignment of treatment, given covariates W, has nothing to do with potential outcomes.
Conditional mediator randomization.

\[ Y(a,m) \perp M|W, A=a \tag{45.5}\]

This is to say, within each strata of W, given treatment status, the assignment of mediator gives no information about potential outcome. This randomization is usually not implemented during many experiments (random trials). This is the assumption that makes a lot of mediation hard to make a causal claim.
Positivity.

There are two positivity assumptions:

\[ P(A=a | W=w) > 0 \tag{45.6}\] for all w.

This is the usual positivity assumption.

\[ P(M=m | W=w, A=a) > 0 \tag{45.7}\]

for all w. This is the mediator positivity.

If these assumptions can be met, then we can do causal mediation analysis. In reality, the ignorability assumptions are untestable. Positivity is not testable either, though it is partly diagnosable: estimated propensity distributions and observed cell counts reveal practical violations, and a fitted propensity close to zero or one in some region is a warning. What cannot be established from a finite sample is the population condition itself. Consistency is the third requirement, alongside exchangeability and positivity.

45.2.1.2 Estimation

CDE can be estimated using G-computation method, or IPW, or the doubly-robust methods, one of them is AIPW (augmented IPW).

G-computation is to model the outcome equation:

\[ CDE_G(m) = E_W \bigl[ Q(1, m, W) - Q(0, m, W) \bigr], \qquad Q(a, m, w) = E[Y \mid A=a, M=m, W=w] \tag{45.8}\] \[ CDE_G(m) = \sum_w \bigl[ Q(1, m, w) - Q(0, m, w) \bigr] P(W=w) \tag{45.9}\]

Writing the inner conditional mean as $Q$ keeps the two operations apart: $Q$ is a regression on observed variables, and the outer expectation is over the covariate distribution. Both equations are the same statement, the second spelling out the outer expectation for discrete $W$.

IPW is to model the treatment assignment equation and the mediator assignment equation:

\[ E[Y(a,m)] = E[{\frac{I(A=a, M=m)}{P(A=a, M=m | W)}} Y] \tag{45.10}\]

Therefore,

\[ CDE_{ipw}(m) = E[{\frac{I(A=1,M=m)}{g_m(m|1,W)g_A(1|W)} - \frac{I(A=0,M=m)}{g_m(m|0,W)g_A(0|W)}} Y] \tag{45.11}\]

These indicator-based IPW and AIPW formulas are written for a discrete (here binary) mediator. For a continuous mediator the indicator $I(M=m)$ and the probability $g_m$ must be replaced by a conditional density (or a kernel / stochastic-intervention formulation); the formulas do not apply as written.

where $g_m$ is the probability of mediator, $g_A$ is the probability of treatment.

For the AIPW estimator:

\[ CDE_{AIPW} = CDE_G(m) + B(\bar Q, g_M, g_A) \tag{45.12}\] where $\bar Q$ is the mean outcome function.

\[ B(\bar Q, g_M, g_A) = \frac{1}{n} \sum_n {\frac{I(M=m, A=1)}{g_m(m|1, W) g_A(1|W)}[Y-\bar Q(m,1,W)]} - \frac{1}{n} \sum_n {\frac{I(M=m, A=0)}{g_m(m|0, W) g_A(0|W)}[Y-\bar Q(m,0,W)]} \tag{45.13}\]

45.2.1.3 simulation example

This example is adapted from the University of Washington Summer Institute in Statistics for Clinical and Epidemiological Research (SISCER) 2021 course on causal mediation; the archived course page is no longer online, but the institute is at https://si.biostat.washington.edu/.

We simulate $n = 5000$ observations. There are two independent standard normal covariates, $W_1$ confounding the treatment-outcome relationship and $W_2$ confounding the mediator-outcome relationship. Treatment is $A \sim \text{Bernoulli}(\Lambda(-1 + W_1/2))$, the mediator is $M \sim \text{Bernoulli}(\Lambda(-2 + A/2 + W_2/3))$, and the outcome is $Y \sim \text{Bernoulli}(\Lambda(-1 + A - M/2 + W_1/3 + W_2/3))$. All three are binary.

The target is $CDE(0)$, the effect of treatment holding the mediator at 0 for everyone. Because the DGP is fully specified we can integrate it out rather than simulate: $CDE(0) = E_W[\Lambda(W_1/3 + W_2/3) - \Lambda(-1 + W_1/3 + W_2/3)] = 0.222$.

First, G-computation:

Code

set.seed(1234)
n <- 5000
# confounder of A/Y
W1 <- rnorm(n)
# confounder of M/Y
W2 <- rnorm(n)
# treatment
A <- rbinom(n, 1, plogis(-1 + W1 / 2))
# binary mediator
M <- rbinom(n, 1, plogis(-2 + A / 2 + W2 / 3))
# binary outcome
Y <- rbinom(n, 1, plogis(-1 + A - M / 2 + W1 / 3 + W2 / 3))
full_data <- data.frame(W1 = W1, W2 = W2, A = A, M = M, Y = Y)

# fit outcome regression
or_fit <- glm(Y ~ A + M + W1 + W2, family = binomial(), data = full_data)
# new data setting A and M
data_A1_M0 <- data_A0_M0 <- full_data
data_A1_M0$A <- 1; data_A1_M0$M <- 0
data_A0_M0$A <- 0; data_A0_M0$M <- 0
# predict on new data
Qbar_A1_M0 <- predict(or_fit, newdata = data_A1_M0, type = "response")
Qbar_A0_M0 <- predict(or_fit, newdata = data_A0_M0, type = "response")
# gcomp estimate of CDE(0)
mean(Qbar_A1_M0 - Qbar_A0_M0)

[1] 0.2515527

G-computation gives 0.2516.

Second, IPW:

Code

# model for P(A = 1 | W)
ps_fit1 <- glm(A ~ W1 + W2, family = binomial(), data = full_data)
P_A1_W <- predict(ps_fit1, type = "response")
P_A0_W <- 1 - P_A1_W
# model for P(M = 0 | A, W)
ps_fit2 <- glm(M ~ A + W1 + W2, family = binomial(), data = full_data)
# P(M = 0 | A = 1, W)
data_A1 <- full_data; data_A1$A <- 1
P_M0_A1_W <- 1 - predict(ps_fit2, newdata = data_A1, type = "response")
# P(M = 0 | A = 0, W)
data_A0 <- full_data; data_A0$A <- 0
P_M0_A0_W <- 1 - predict(ps_fit2, newdata = data_A0, type = "response")
# ipw estimate of CDE(0)
mean( (A == 1) / P_A1_W * (M == 0) / P_M0_A1_W * Y ) -
mean( (A == 0) / P_A0_W * (M == 0) / P_M0_A0_W * Y )

[1] 0.2553091

IPW gives 0.2553.

AIPW:

Code

# aipw estimate of E[Y(1,0)]
aiptw_EY_A1_M0 <- mean(Qbar_A1_M0) +
mean( (A == 1) / P_A1_W * (M == 0) / P_M0_A1_W * (Y - Qbar_A1_M0) )
# aipw estimate of E[Y(0,0)]
aiptw_EY_A0_M0 <- mean(Qbar_A0_M0) +
mean( (A == 0) / P_A0_W * (M == 0) / P_M0_A0_W * (Y - Qbar_A0_M0) )
# aipw estimate of CDE(0)
aiptw_EY_A1_M0 - aiptw_EY_A0_M0

[1] 0.2554265

AIPW gives 0.2554.

The three estimates cluster at 0.25 against a true 0.222, all about 0.03 high. That is one unlucky draw rather than bias: every model here is correctly specified, and over 400 replications of this design the g-computation estimator averages 0.2206 with a sampling standard deviation of 0.0159, so this seed sits about 1.9 standard deviations above the truth. Note that the three agree with each other far more closely than any of them agrees with the truth — they share the same data and the same fitted nuisance models, so their agreement checks the implementation, not the sampling error.

45.2.2 Natural direct and indirect effect

CDE is to study the effect of treatment, given the level of mediator. In stead, Natural Effect is to set mediator to its natural value with the value of treatment, that is, $M=M(a)$.

\[ ATE = NIE + NDE = (E[ Y(1,M(1))] - E[Y(1,M(0))]) + (E[Y(1, M(0))] - E[Y(0,M(0))]) \tag{45.14}\]

There are two different ways to decompose the $ATE$ in terms of $NIE$ and $NDE$. Above is one of them.

The advantage of NDE and NIE comparing to CDE is that it’s more “natural”; that is, you don’t set the level of mediator deterministically. And it can decompose the ATE into direct and indirect effects.

However, there is an additional assumption for NDE and NIE identified.

Assumption 4:

\[ Y(a, m) \perp M(a^*) | W \tag{45.15}\]

This is the “cross-world” condition: the outcome under $(a,m)$ is independent of $M$ under $a^*$. These two situations cannot happen in the same world; you cannot set $A$ to both $a$ and $a^*$. No experiment can implement it.

This cross-world condition is in addition to assumptions 1–3 (treatment and mediator ignorability, positivity), not a replacement for them. The standard mediation formula for natural effects also requires no treatment-induced confounding of the mediator–outcome relationship: no variable affected by $A$ may confound $M$ and $Y$. When such a confounder exists, natural effects are not identified by this formula and interventional direct/indirect effects are used instead.

Simulation example from the same website materials:

Code

# fit outcome regression (include interaction because we can)
or_fit <- glm(Y ~ A + M + W1 + W2 + A*M + M*W1,
family = binomial(), data = full_data)
# need E(Y | A = 0/1, M = 0/1, W1 = W1i, W2 = W2i)
get_EY_a_m_Wi <- function(full_data, or_fit, a, m){
data_Aa_Mm_Wi <- full_data
data_Aa_Mm_Wi$A <- a; data_Aa_Mm_Wi$M <- m
predict(or_fit, newdata = data_Aa_Mm_Wi, type = "response")
}
EY_A0_M0_Wi <- get_EY_a_m_Wi(full_data, or_fit, a = 0, m = 0)
EY_A0_M1_Wi <- get_EY_a_m_Wi(full_data, or_fit, a = 0, m = 1)
EY_A1_M0_Wi <- get_EY_a_m_Wi(full_data, or_fit, a = 1, m = 0)
EY_A1_M1_Wi <- get_EY_a_m_Wi(full_data, or_fit, a = 1, m = 1)

# include interactions -- why not?
med_fit <- glm(M ~ A*W1 + W1*W2, family = binomial(), data = full_data)
# estimates of P(M = m | A = a, W = W_i)
get_Pm_a_Wi <- function(full_data, med_fit, a, m){
data_Aa_Wi <- full_data; data_Aa_Wi$A <- a
p <- predict(med_fit, newdata = data_Aa_Wi, type = "response")
if(m == 1){
p
}else{
1 - p
}
}
PM0_A0_Wi <- get_Pm_a_Wi(full_data, med_fit, a = 0, m = 0)
PM1_A0_Wi <- get_Pm_a_Wi(full_data, med_fit, a = 0, m = 1)
PM0_A1_Wi <- get_Pm_a_Wi(full_data, med_fit, a = 1, m = 0)
PM1_A1_Wi <- get_Pm_a_Wi(full_data, med_fit, a = 1, m = 1)

# E(E(Y | A = 1, M, W) | A = 1, W)
EY1M1_Wi <- EY_A1_M1_Wi * PM1_A1_Wi + EY_A1_M0_Wi * PM0_A1_Wi
# E(E(Y | A = 0, M, W) | A = 1, W)
EY0M1_Wi <- EY_A0_M1_Wi * PM1_A1_Wi + EY_A0_M0_Wi * PM0_A1_Wi
# E(E(Y | A = 1, M, W) | A = 0, W)
EY1M0_Wi <- EY_A1_M1_Wi * PM1_A0_Wi + EY_A1_M0_Wi * PM0_A0_Wi
# E(E(Y | A = 0, M, W) | A = 0, W)
EY0M0_Wi <- EY_A0_M1_Wi * PM1_A0_Wi + EY_A0_M0_Wi * PM0_A0_Wi

# estimate of E[Y(1, M(1))]
E_Y1M1 <- mean(EY1M1_Wi)
# estimate of E[Y(0, M(1))]
E_Y0M1 <- mean(EY0M1_Wi)
# estimate of E[Y(1, M(0))]
E_Y1M0 <- mean(EY1M0_Wi)
# estimate of E[Y(0, M(0))]
E_Y0M0 <- mean(EY0M0_Wi)

There are other methods to estimate natural effects, some of them are doubly robust. The “mediation” package can be used to do it.

45.2.3 When there is an exposure-induced confounder

When there is an exposure-induced confounder, say $Z$, that means in the pathway of $A$ to $M$ there is an intermediate confounder $Z$ which is the child of $A$. When that happens and we observe $Z$, $CDE$ can still be estimated, but $NIE$ and $NDE$ cannot. The reason is that the presence of $Z$ implies that the cross-world assumption cannot hold.

45.2.4 Interventional direct and indirect effects

People are not happy with the cross-world assumption in general. Interventional direct and indirect effects are introduced to avoid this assumption and still be able to decompose the $ATE$.

\[ IIE + IDE = (E[ Y(1,G_1)] - E[Y(1,G_0)]) + (E[Y(1,G_0)]-E[Y(0,G_0)]) = E[Y(1,G_1)] - E[Y(0,G_0)] \tag{45.16}\]

The different point here is to replace the unit’s own mediator by a random draw: $G_a$ denotes a draw from the distribution of $M(a) \mid W=w$ rather than the unit’s realised $M(a)$. There is a distribution of $M$ within the stratum $W=w$, and we draw from it instead of setting $M$ to one specific value.

Note carefully what this sum is, because it is easy to write down as the $ATE$ and it is not. $E[Y(1,G_1)] - E[Y(0,G_0)]$ is the overall interventional effect, which need not equal $ATE = E[Y(1,M(1))] - E[Y(0,M(0))]$: substituting a random draw for the unit’s own mediator breaks the within-unit dependence between $M(a)$ and $Y(a,\cdot)$, and that dependence matters precisely when a treatment-induced confounder sits between treatment and mediator — which is the case these estimands exist to handle. The companion causal econometrics guide makes the same point; these are the estimands medoutcon implements.

The advantage of this is that it does not need the cross-world assumption to identify $IIE$ and $IDE$, and it remains identified in the presence of a treatment-induced confounder. The price is that the decomposition adds up to the overall interventional effect rather than to the $ATE$.

45.2.4.1 an example

Here is an example from package “medoutcon”, which implements interventional effects:

Code

library(data.table)
library(tidyverse)
library(medoutcon)
set.seed(1584)

# produces a simple data set based on ca causal model with mediation
make_example_data <- function(n_obs = 1000) {
  ## baseline covariates
  w_1 <- rbinom(n_obs, 1, prob = 0.6)
  w_2 <- rbinom(n_obs, 1, prob = 0.3)
  w_3 <- rbinom(n_obs, 1, prob = pmin(0.2 + (w_1 + w_2) / 3, 1))
  w <- cbind(w_1, w_2, w_3)
  w_names <- paste("W", seq_len(ncol(w)), sep = "_")

  ## exposure
  a <- as.numeric(rbinom(n_obs, 1, plogis(rowSums(w) - 2)))

  ## mediator-outcome confounder affected by treatment
  z <- rbinom(n_obs, 1, plogis(rowMeans(-log(2) + w - a) + 0.2))

  ## mediator -- could be multivariate
  m <- rbinom(n_obs, 1, plogis(rowSums(log(3) * w[, -3] + a - z)))
  m_names <- "M"

  ## outcome
  y <- rbinom(n_obs, 1, plogis(1 / (rowSums(w) - z + a + m)))

  ## construct output
  dat <- as.data.table(cbind(w = w, a = a, z = z, m = m, y = y))
  setnames(dat, c(w_names, "A", "Z", m_names, "Y"))
  return(dat)
}

# set seed and simulate example data
example_data <- make_example_data()
w_names <- str_subset(colnames(example_data), "W")
m_names <- str_subset(colnames(example_data), "M")

# quick look at the data
head(example_data)

     W_1   W_2   W_3     A     Z     M     Y
   <num> <num> <num> <num> <num> <num> <num>
1:     1     0     1     0     0     0     1
2:     0     1     0     0     0     1     0
3:     1     1     1     1     0     1     1
4:     0     1     1     0     0     1     0
5:     0     0     0     0     0     1     1
6:     1     0     1     1     0     1     0

Code

# compute one-step estimate of the interventional direct effect
os_de <- medoutcon(W = example_data[, ..w_names],
                   A = example_data$A,
                   Z = example_data$Z,
                   M = example_data[, ..m_names],
                   Y = example_data$Y,
                   effect = "direct",
                   estimator = "onestep")
os_de

# compute targeted minimum loss estimate of the interventional direct effect
tmle_de <- medoutcon(W = example_data[, ..w_names],
                     A = example_data$A,
                     Z = example_data$Z,
                     M = example_data[, ..m_names],
                     Y = example_data$Y,
                     effect = "direct",
                     estimator = "tmle")
tmle_de

45.2.4.2 potential problem

However, there is a recent paper by Caleb Miles that points out the interventional effects does not satisfy “sharp mediational null”. “Intuitively, an indirect effect measure should be null whenever there is no individual-level effect for anyone in the population.” He found a counter example that the interventional effects cannot pass this “sharp mediational null” criterion.

More recently, Ivan Diaz and coauthors have proposed alternative mediation estimands (discussed under the name “causal influence” in some of that line of work) aimed at satisfying the sharp mediational null where interventional effects fail to. We have not worked through the details of that proposal here, and won’t attempt to summarize its identification conditions or estimator without first working through the derivation carefully – readers interested in this specific gap should go directly to Miles’s counterexample paper and the more recent Diaz et al. work rather than relying on a secondhand summary.

45.3 Summary

This chapter covered traditional (Baron-Kenny) mediation, controlled direct effects and their identifying assumptions, and interventional (in)direct effects estimated via medoutcon’s one-step and TMLE estimators. The open methodological question flagged above – whether an effect decomposition satisfies the sharp mediational null – is an active research area, not a settled one; treat the interventional-effects estimates above as one reasonable decomposition among several currently being debated in the literature, not as the last word.

Systematic treatment: R · Julia.

--- title: "Causal mediation analysis" date: "2022-07-21" --- ## Traditional mediation analysis We summarized how to do traditional mediation analysis (Baron and Kenny, 1986) with R and Stata in my last blog about mediation. Traditionally mediation model can be represented in the following equations: $$ Y = a A + b M + \epsilon_1 $$ {#eq-causal-mediation-1} $$ M = c A + \epsilon_2 $$ {#eq-causal-mediation-2} There are three main variables, A (treatment, or exposure), M (mediator) and Y (outcome). There are two pathways from A to Y, a direct effect, and an indirect effect through A -\> M -\> Y. The traditional approach is to model these two equations directly, and get the estimates of direct and indirect effects. ## Causal mediation analysis The traditional method is criticized for lack of causal interpretation. The causal mediation analysis makes assumptions up front. If those assumptions are met, then you have a causal interpretation. In a random experiment, usually A is randomly assigned, but M is not. If there is a confounder on the M to Y pathway, then we don't have a causal effect, or say we cannot estimate causal effect without confounder bias. Causal mediation analysis needs that too; it's not like if you use a causal mediation package then you get the causal effect. The other shortcoming of the traditional mediation method is that it assumes a homogeneous treatment indirect effect; that is, there is no $A$ - $M$ interaction. With causal mediation, this interaction can be allowed. ### Controlled direct effect Suppose we can set $A$ and $M$ at will to any $(a, m)$, then we have the potential outcome $Y(a,m)$. Controlled direct effect (CDE) is defined as $$ CDE(m) = E(Y(1,m)-Y(0,m)) $$ {#eq-causal-mediation-3} That is, setting $M$ to $m$, what is the effect of $A$ on $Y$? #### Assumptions What are the assumptions for CDE to be identified (estimable)? 1. Conditional treatment randomization. Suppose we observe confounders W, which can be a joint set of confounders for the A to Y pathway, or the M to Y pathway. $$ Y(a,m) \perp A|W $$ {#eq-causal-mediation-4} This is the usual conditional independence (unconfoundedness, ignorability, etc.) assumption. This is saying the assignment of treatment, given covariates W, has nothing to do with potential outcomes. 2. Conditional mediator randomization. $$ Y(a,m) \perp M|W, A=a $$ {#eq-causal-mediation-5} This is to say, within each strata of W, given treatment status, the assignment of mediator gives no information about potential outcome. This randomization is usually not implemented during many experiments (random trials). This is the assumption that makes a lot of mediation hard to make a causal claim. 3. Positivity. There are two positivity assumptions: $$ P(A=a | W=w) > 0 $$ {#eq-causal-mediation-6} for all w. This is the usual positivity assumption. $$ P(M=m | W=w, A=a) > 0 $$ {#eq-causal-mediation-7} for all w. This is the mediator positivity. If these assumptions can be met, then we can do causal mediation analysis. In reality, the ignorability assumptions are untestable. Positivity is not testable either, though it is partly diagnosable: estimated propensity distributions and observed cell counts reveal practical violations, and a fitted propensity close to zero or one in some region is a warning. What cannot be established from a finite sample is the population condition itself. Consistency is the third requirement, alongside exchangeability and positivity. #### Estimation CDE can be estimated using G-computation method, or IPW, or the doubly-robust methods, one of them is AIPW (augmented IPW). G-computation is to model the outcome equation: $$ CDE_G(m) = E_W \bigl[ Q(1, m, W) - Q(0, m, W) \bigr], \qquad Q(a, m, w) = E[Y \mid A=a, M=m, W=w] $$ {#eq-causal-mediation-8} $$ CDE_G(m) = \sum_w \bigl[ Q(1, m, w) - Q(0, m, w) \bigr] P(W=w) $$ {#eq-causal-mediation-9} Writing the inner conditional mean as $Q$ keeps the two operations apart: $Q$ is a regression on observed variables, and the outer expectation is over the covariate distribution. Both equations are the same statement, the second spelling out the outer expectation for discrete $W$. IPW is to model the treatment assignment equation and the mediator assignment equation: $$ E[Y(a,m)] = E[{\frac{I(A=a, M=m)}{P(A=a, M=m | W)}} Y] $$ {#eq-causal-mediation-10} Therefore, $$ CDE_{ipw}(m) = E[{\frac{I(A=1,M=m)}{g_m(m|1,W)g_A(1|W)} - \frac{I(A=0,M=m)}{g_m(m|0,W)g_A(0|W)}} Y] $$ {#eq-causal-mediation-11} These indicator-based IPW and AIPW formulas are written for a **discrete (here binary) mediator**. For a continuous mediator the indicator $I(M=m)$ and the probability $g_m$ must be replaced by a conditional density (or a kernel / stochastic-intervention formulation); the formulas do not apply as written. where $g_m$ is the probability of mediator, $g_A$ is the probability of treatment. For the AIPW estimator: $$ CDE_{AIPW} = CDE_G(m) + B(\bar Q, g_M, g_A) $$ {#eq-causal-mediation-12} where $\bar Q$ is the mean outcome function. $$ B(\bar Q, g_M, g_A) = \frac{1}{n} \sum_n {\frac{I(M=m, A=1)}{g_m(m|1, W) g_A(1|W)}[Y-\bar Q(m,1,W)]} - \frac{1}{n} \sum_n {\frac{I(M=m, A=0)}{g_m(m|0, W) g_A(0|W)}[Y-\bar Q(m,0,W)]} $$ {#eq-causal-mediation-13} #### simulation example This example is adapted from the University of Washington Summer Institute in Statistics for Clinical and Epidemiological Research (SISCER) 2021 course on causal mediation; the archived course page is no longer online, but the institute is at <https://si.biostat.washington.edu/>. We simulate $n = 5000$ observations. There are two independent standard normal covariates, $W_1$ confounding the treatment-outcome relationship and $W_2$ confounding the mediator-outcome relationship. Treatment is $A \sim \text{Bernoulli}(\Lambda(-1 + W_1/2))$, the mediator is $M \sim \text{Bernoulli}(\Lambda(-2 + A/2 + W_2/3))$, and the outcome is $Y \sim \text{Bernoulli}(\Lambda(-1 + A - M/2 + W_1/3 + W_2/3))$. All three are binary. The target is $CDE(0)$, the effect of treatment holding the mediator at 0 for everyone. Because the DGP is fully specified we can integrate it out rather than simulate: $CDE(0) = E_W[\Lambda(W_1/3 + W_2/3) - \Lambda(-1 + W_1/3 + W_2/3)] = 0.222$. First, G-computation: ```{r} #| echo: true #| message: false set.seed(1234) n <- 5000 # confounder of A/Y W1 <- rnorm(n) # confounder of M/Y W2 <- rnorm(n) # treatment A <- rbinom(n, 1, plogis(-1 + W1 / 2)) # binary mediator M <- rbinom(n, 1, plogis(-2 + A / 2 + W2 / 3)) # binary outcome Y <- rbinom(n, 1, plogis(-1 + A - M / 2 + W1 / 3 + W2 / 3)) full_data <- data.frame(W1 = W1, W2 = W2, A = A, M = M, Y = Y) # fit outcome regression or_fit <- glm(Y ~ A + M + W1 + W2, family = binomial(), data = full_data) # new data setting A and M data_A1_M0 <- data_A0_M0 <- full_data data_A1_M0$A <- 1; data_A1_M0$M <- 0 data_A0_M0$A <- 0; data_A0_M0$M <- 0 # predict on new data Qbar_A1_M0 <- predict(or_fit, newdata = data_A1_M0, type = "response") Qbar_A0_M0 <- predict(or_fit, newdata = data_A0_M0, type = "response") # gcomp estimate of CDE(0) mean(Qbar_A1_M0 - Qbar_A0_M0) ``` G-computation gives 0.2516. Second, IPW: ```{r} #| echo: true #| message: false # model for P(A = 1 | W) ps_fit1 <- glm(A ~ W1 + W2, family = binomial(), data = full_data) P_A1_W <- predict(ps_fit1, type = "response") P_A0_W <- 1 - P_A1_W # model for P(M = 0 | A, W) ps_fit2 <- glm(M ~ A + W1 + W2, family = binomial(), data = full_data) # P(M = 0 | A = 1, W) data_A1 <- full_data; data_A1$A <- 1 P_M0_A1_W <- 1 - predict(ps_fit2, newdata = data_A1, type = "response") # P(M = 0 | A = 0, W) data_A0 <- full_data; data_A0$A <- 0 P_M0_A0_W <- 1 - predict(ps_fit2, newdata = data_A0, type = "response") # ipw estimate of CDE(0) mean( (A == 1) / P_A1_W * (M == 0) / P_M0_A1_W * Y ) - mean( (A == 0) / P_A0_W * (M == 0) / P_M0_A0_W * Y ) ``` IPW gives 0.2553. AIPW: ```{r} #| echo: true #| message: false # aipw estimate of E[Y(1,0)] aiptw_EY_A1_M0 <- mean(Qbar_A1_M0) + mean( (A == 1) / P_A1_W * (M == 0) / P_M0_A1_W * (Y - Qbar_A1_M0) ) # aipw estimate of E[Y(0,0)] aiptw_EY_A0_M0 <- mean(Qbar_A0_M0) + mean( (A == 0) / P_A0_W * (M == 0) / P_M0_A0_W * (Y - Qbar_A0_M0) ) # aipw estimate of CDE(0) aiptw_EY_A1_M0 - aiptw_EY_A0_M0 ``` AIPW gives 0.2554. The three estimates cluster at 0.25 against a true 0.222, all about 0.03 high. That is one unlucky draw rather than bias: every model here is correctly specified, and over 400 replications of this design the g-computation estimator averages 0.2206 with a sampling standard deviation of 0.0159, so this seed sits about 1.9 standard deviations above the truth. Note that the three agree with each other far more closely than any of them agrees with the truth — they share the same data and the same fitted nuisance models, so their agreement checks the implementation, not the sampling error. ### Natural direct and indirect effect CDE is to study the effect of treatment, given the level of mediator. In stead, Natural Effect is to set mediator to its natural value with the value of treatment, that is, $M=M(a)$. $$ ATE = NIE + NDE = (E[ Y(1,M(1))] - E[Y(1,M(0))]) + (E[Y(1, M(0))] - E[Y(0,M(0))]) $$ {#eq-causal-mediation-14} There are two different ways to decompose the $ATE$ in terms of $NIE$ and $NDE$. Above is one of them. The advantage of NDE and NIE comparing to CDE is that it's more "natural"; that is, you don't set the level of mediator deterministically. And it can decompose the ATE into direct and indirect effects. However, there is an additional assumption for NDE and NIE identified. Assumption 4: $$ Y(a, m) \perp M(a^*) | W$$ {#eq-causal-mediation-15} This is the "cross-world" condition: the outcome under $(a,m)$ is independent of $M$ under $a^*$. These two situations cannot happen in the same world; you cannot set $A$ to both $a$ and $a^*$. No experiment can implement it. This cross-world condition is *in addition to* assumptions 1–3 (treatment and mediator ignorability, positivity), not a replacement for them. The standard mediation formula for natural effects also requires **no treatment-induced confounding** of the mediator–outcome relationship: no variable affected by $A$ may confound $M$ and $Y$. When such a confounder exists, natural effects are not identified by this formula and interventional direct/indirect effects are used instead. Simulation example from the same website materials: ```{r} #| echo: true #| message: false # fit outcome regression (include interaction because we can) or_fit <- glm(Y ~ A + M + W1 + W2 + A*M + M*W1, family = binomial(), data = full_data) # need E(Y | A = 0/1, M = 0/1, W1 = W1i, W2 = W2i) get_EY_a_m_Wi <- function(full_data, or_fit, a, m){ data_Aa_Mm_Wi <- full_data data_Aa_Mm_Wi$A <- a; data_Aa_Mm_Wi$M <- m predict(or_fit, newdata = data_Aa_Mm_Wi, type = "response") } EY_A0_M0_Wi <- get_EY_a_m_Wi(full_data, or_fit, a = 0, m = 0) EY_A0_M1_Wi <- get_EY_a_m_Wi(full_data, or_fit, a = 0, m = 1) EY_A1_M0_Wi <- get_EY_a_m_Wi(full_data, or_fit, a = 1, m = 0) EY_A1_M1_Wi <- get_EY_a_m_Wi(full_data, or_fit, a = 1, m = 1) # include interactions -- why not? med_fit <- glm(M ~ A*W1 + W1*W2, family = binomial(), data = full_data) # estimates of P(M = m | A = a, W = W_i) get_Pm_a_Wi <- function(full_data, med_fit, a, m){ data_Aa_Wi <- full_data; data_Aa_Wi$A <- a p <- predict(med_fit, newdata = data_Aa_Wi, type = "response") if(m == 1){ p }else{ 1 - p } } PM0_A0_Wi <- get_Pm_a_Wi(full_data, med_fit, a = 0, m = 0) PM1_A0_Wi <- get_Pm_a_Wi(full_data, med_fit, a = 0, m = 1) PM0_A1_Wi <- get_Pm_a_Wi(full_data, med_fit, a = 1, m = 0) PM1_A1_Wi <- get_Pm_a_Wi(full_data, med_fit, a = 1, m = 1) # E(E(Y | A = 1, M, W) | A = 1, W) EY1M1_Wi <- EY_A1_M1_Wi * PM1_A1_Wi + EY_A1_M0_Wi * PM0_A1_Wi # E(E(Y | A = 0, M, W) | A = 1, W) EY0M1_Wi <- EY_A0_M1_Wi * PM1_A1_Wi + EY_A0_M0_Wi * PM0_A1_Wi # E(E(Y | A = 1, M, W) | A = 0, W) EY1M0_Wi <- EY_A1_M1_Wi * PM1_A0_Wi + EY_A1_M0_Wi * PM0_A0_Wi # E(E(Y | A = 0, M, W) | A = 0, W) EY0M0_Wi <- EY_A0_M1_Wi * PM1_A0_Wi + EY_A0_M0_Wi * PM0_A0_Wi # estimate of E[Y(1, M(1))] E_Y1M1 <- mean(EY1M1_Wi) # estimate of E[Y(0, M(1))] E_Y0M1 <- mean(EY0M1_Wi) # estimate of E[Y(1, M(0))] E_Y1M0 <- mean(EY1M0_Wi) # estimate of E[Y(0, M(0))] E_Y0M0 <- mean(EY0M0_Wi) ``` There are other methods to estimate natural effects, some of them are doubly robust. The "mediation" package can be used to do it. ### When there is an exposure-induced confounder When there is an exposure-induced confounder, say $Z$, that means in the pathway of $A$ to $M$ there is an intermediate confounder $Z$ which is the child of $A$. When that happens and we observe $Z$, $CDE$ can still be estimated, but $NIE$ and $NDE$ cannot. The reason is that the presence of $Z$ implies that the cross-world assumption cannot hold. ### Interventional direct and indirect effects People are not happy with the cross-world assumption in general. Interventional direct and indirect effects are introduced to avoid this assumption and still be able to decompose the $ATE$. $$ IIE + IDE = (E[ Y(1,G_1)] - E[Y(1,G_0)]) + (E[Y(1,G_0)]-E[Y(0,G_0)]) = E[Y(1,G_1)] - E[Y(0,G_0)] $$ {#eq-causal-mediation-16} The different point here is to replace the unit's own mediator by a random draw: $G_a$ denotes a draw from the distribution of $M(a) \mid W=w$ rather than the unit's realised $M(a)$. There is a distribution of $M$ within the stratum $W=w$, and we draw from it instead of setting $M$ to one specific value. Note carefully what this sum is, because it is easy to write down as the $ATE$ and it is not. $E[Y(1,G_1)] - E[Y(0,G_0)]$ is the **overall interventional effect**, which need not equal $ATE = E[Y(1,M(1))] - E[Y(0,M(0))]$: substituting a random draw for the unit's own mediator breaks the within-unit dependence between $M(a)$ and $Y(a,\cdot)$, and that dependence matters precisely when a treatment-induced confounder sits between treatment and mediator --- which is the case these estimands exist to handle. The companion [causal econometrics guide](https://xiangao.github.io/causal_econometrics_guide/mediation.html) makes the same point; these are the estimands `medoutcon` implements. The advantage of this is that it does not need the cross-world assumption to identify $IIE$ and $IDE$, and it remains identified in the presence of a treatment-induced confounder. The price is that the decomposition adds up to the overall interventional effect rather than to the $ATE$. #### an example Here is an example from package "medoutcon", which implements interventional effects: ```{r} #| echo: true #| message: false library(data.table) library(tidyverse) library(medoutcon) set.seed(1584) # produces a simple data set based on ca causal model with mediation make_example_data <- function(n_obs = 1000) { ## baseline covariates w_1 <- rbinom(n_obs, 1, prob = 0.6) w_2 <- rbinom(n_obs, 1, prob = 0.3) w_3 <- rbinom(n_obs, 1, prob = pmin(0.2 + (w_1 + w_2) / 3, 1)) w <- cbind(w_1, w_2, w_3) w_names <- paste("W", seq_len(ncol(w)), sep = "_") ## exposure a <- as.numeric(rbinom(n_obs, 1, plogis(rowSums(w) - 2))) ## mediator-outcome confounder affected by treatment z <- rbinom(n_obs, 1, plogis(rowMeans(-log(2) + w - a) + 0.2)) ## mediator -- could be multivariate m <- rbinom(n_obs, 1, plogis(rowSums(log(3) * w[, -3] + a - z))) m_names <- "M" ## outcome y <- rbinom(n_obs, 1, plogis(1 / (rowSums(w) - z + a + m))) ## construct output dat <- as.data.table(cbind(w = w, a = a, z = z, m = m, y = y)) setnames(dat, c(w_names, "A", "Z", m_names, "Y")) return(dat) } # set seed and simulate example data example_data <- make_example_data() w_names <- str_subset(colnames(example_data), "W") m_names <- str_subset(colnames(example_data), "M") # quick look at the data head(example_data) # compute one-step estimate of the interventional direct effect os_de <- medoutcon(W = example_data[, ..w_names], A = example_data$A, Z = example_data$Z, M = example_data[, ..m_names], Y = example_data$Y, effect = "direct", estimator = "onestep") os_de # compute targeted minimum loss estimate of the interventional direct effect tmle_de <- medoutcon(W = example_data[, ..w_names], A = example_data$A, Z = example_data$Z, M = example_data[, ..m_names], Y = example_data$Y, effect = "direct", estimator = "tmle") tmle_de ``` #### potential problem However, there is a recent paper by Caleb Miles that points out the interventional effects does not satisfy "sharp mediational null". "Intuitively, an indirect effect measure should be null whenever there is no individual-level effect for anyone in the population." He found a counter example that the interventional effects cannot pass this "sharp mediational null" criterion. More recently, Ivan Diaz and coauthors have proposed alternative mediation estimands (discussed under the name "causal influence" in some of that line of work) aimed at satisfying the sharp mediational null where interventional effects fail to. We have not worked through the details of that proposal here, and won't attempt to summarize its identification conditions or estimator without first working through the derivation carefully -- readers interested in this specific gap should go directly to Miles's counterexample paper and the more recent Diaz et al. work rather than relying on a secondhand summary. ## Summary This chapter covered traditional (Baron-Kenny) mediation, controlled direct effects and their identifying assumptions, and interventional (in)direct effects estimated via `medoutcon`'s one-step and TMLE estimators. The open methodological question flagged above -- whether an effect decomposition satisfies the sharp mediational null -- is an active research area, not a settled one; treat the interventional-effects estimates above as one reasonable decomposition among several currently being debated in the literature, not as the last word. ---  *Systematic treatment: [R](https://xiangao.github.io/causal_econometrics_guide/mediation.html) · [Julia](https://xiangao.github.io/causal_econometrics_julia/mediation.html).*