Traditionally mediation model can be represented in the following equations:
\[ Y = a X + b M + \epsilon_1 \]\[ M = c X + \epsilon_2 \]
That is, we’d like to study the effect of \(X\) on \(Y\), and we see the effect can be a direct effect, and an indirect effect, through \(M\).
Baron and Kenny’s (http://davidakenny.net/cm/mediate.htm) method is done in four steps. Modern approach tends to use SEM (structural equation modeling) to model these two equations directly.
R’s lavaan and Stata’s sem commands are powerful tools.
A simple example:
library(lavaan)set.seed(1234)n<-10000X<-rnorm(n)M<-0.5*X+rnorm(n)Y<-0.7*M+0.3*X+rnorm(n)Data<-data.frame(X =X, Y =Y, M =M)model<-' Y ~ a*X + b*M M ~ c*X # direct effect (a) # indirect effect (b*c) bc := b*c # total effect total := a + (b*c) 'fit<-sem(model, data =Data)summary(fit)
lavaan 0.6-21 ended normally after 1 iteration
Estimator ML
Optimization method NLMINB
Number of model parameters 5
Number of observations 10000
Model Test User Model:
Test statistic 0.000
Degrees of freedom 0
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Regressions:
Estimate Std.Err z-value P(>|z|)
Y ~
X (a) 0.286 0.011 25.233 0.000
M (b) 0.699 0.010 70.384 0.000
M ~
X (c) 0.511 0.010 50.161 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.Y 0.998 0.014 70.711 0.000
.M 1.012 0.014 70.711 0.000
Defined Parameters:
Estimate Std.Err z-value P(>|z|)
bc 0.357 0.009 40.849 0.000
total 0.643 0.012 51.956 0.000
Now let’s do it in stata.
setobs 10000gen x= rnormal()genm= .5*x + rnormal()geny= .7*m + .3*x +rnormal()sem (y <- x m) (m <- x)estat teffects
Number of observations (_N) was 0, now 10,000.
Endogenous variables
Observed: y m
Exogenous variables
Observed: x
Fitting target model:
Iteration 0: Log likelihood = -42567.918
Iteration 1: Log likelihood = -42567.918
Structural equation model Number of obs = 10,000
Estimation method: ml
Log likelihood = -42567.918
------------------------------------------------------------------------------
| OIM
| Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
Structural |
y |
m | .7136545 .0101594 70.25 0.000 .6937424 .7335666
x | .2971221 .0113774 26.12 0.000 .2748228 .3194215
_cons | -.0000391 .0101082 -0.00 0.997 -.0198508 .0197725
-----------+----------------------------------------------------------------
m |
x | .5038605 .0100014 50.38 0.000 .4842581 .5234628
_cons | .0187204 .0099478 1.88 0.060 -.0007769 .0382177
-------------+----------------------------------------------------------------
var(e.y)| 1.021391 .0144446 .9934684 1.050098
var(e.m)| .9895863 .0139949 .9625336 1.017399
------------------------------------------------------------------------------
LR test of model vs. saturated: chi2(0) = 0.00 Prob > chi2 = .
Direct effects
------------------------------------------------------------------------------
| OIM
| Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
Structural |
y |
m | .7136545 .0101594 70.25 0.000 .6937424 .7335666
x | .2971221 .0113774 26.12 0.000 .2748228 .3194215
-----------+----------------------------------------------------------------
m |
x | .5038605 .0100014 50.38 0.000 .4842581 .5234628
------------------------------------------------------------------------------
Indirect effects
------------------------------------------------------------------------------
| OIM
| Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
Structural |
y |
m | 0 (no path)
x | .3595823 .0087834 40.94 0.000 .3423672 .3767974
-----------+----------------------------------------------------------------
m |
x | 0 (no path)
------------------------------------------------------------------------------
Total effects
------------------------------------------------------------------------------
| OIM
| Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
Structural |
y |
m | .7136545 .0101594 70.25 0.000 .6937424 .7335666
x | .6567044 .0124172 52.89 0.000 .6323672 .6810417
-----------+----------------------------------------------------------------
m |
x | .5038605 .0100014 50.38 0.000 .4842581 .5234628
------------------------------------------------------------------------------
The above examples should have direct effect of .3 and indirect effect of .35, and total effect of .65.
28.2 Causal Mediation analysis
The traditional mediation analysis has been criticized for the lack of causal interpretation. Without manipulation of the mediator, it is hard to interpret the effects causally, because even if the treatment is from random experiments, the mediator is often not. Therefore there could be an unmeasured confounder that is causing both \(M\) and \(Y\).
R’s “mediation” package is for causal mediation analysis. It uses simulation to estimate the causal effects of treatment, under assumptions of sequential ignorability.
NoteIdentifying assumptions for natural (in)direct effects
Before reading the mediate() output causally, the following must hold (this is “sequential ignorability” plus the standard support and link conditions):
Treatment ignorability — no unmeasured treatment–outcome confounding given covariates \(W\).
Mediator ignorability — no unmeasured mediator–outcome confounding given \(W\) and treatment.
No treatment-induced confounding — no variable affected by the treatment confounds the mediator–outcome relationship (this is the strong, often-violated condition).
Positivity — both treatment and mediator have positive probability across \(W\).
Consistency — observed outcomes equal the corresponding potential outcomes.
Cross-world independence — the condition underlying natural effects, which no single experiment can test.
Mediator ignorability and the no-treatment-induced-confounding condition are the ones that usually break, because the mediator is rarely randomized.
It estimates the following quantities:
\[\tau_i = Y_i(1, M_i(1)) - Y_i(0, M_i(0))\]
This is the total treatment effect, which is to say, what’s the change in \(Y\) if we change each unit from control to treated, hypothetically?
Then this can be decomposed into the causal mediation effects:
\[\delta_i(t) = Y_i(t, M_i(1)) - Y_i(t, M_i(0))\]
for treatment status \(t=0,1\). This is to say, given treatment status, what’s the mediation effect?
\[\eta_i(t) = Y_i(1, M_i(t)) - Y_i(0, M_i(t)) \]
for treatment status \(t=0,1\). This is to say, given mediator status for each treatment status, what’s the direct effect?
Therefore there are four quantities estimated, direct and mediation effect for treated and control.
R’s “mediation” needs users to feed two models, outcome model and mediation model.
If we study the same data, we would expect it returns the same estimates as the traditional methods. However, the causal mediation models can be much more flexible in outcome and mediation models.
For example, in the case of binary outcome, the traditional approach will have difficulties. We can estimate the outcome model and mediator model jointly, but the total effects are not easy to decompose into direct and indirect effect (see Imai et al, page 320 https://imai.fas.harvard.edu/research/files/BaronKenny.pdf).
The causal mediation analysis framework is much more general.
library(mediation)data(framing)med.fit<-lm(emo~treat+age+educ+gender+income, data =framing)out.fit<-glm(cong_mesg~emo+treat+age+educ+gender+income, data =framing, family =binomial("probit"))set.seed(123)med.out<-mediate(med.fit, out.fit, treat ="treat", mediator ="emo", robustSE =TRUE, sims =100)summary(med.out)
“mediation” package has more functionalities, such as multilevel, interaction of treatment and mediator, etc.
Stata’s sem and gsem commands can model different situations, but the direct effect and indirect effects are not easy to compute, especially when you have binary outcome, or other non-continuous outcome situations. They are not designed for causal mediation analysis.
28.4 Going further: formal causal mediation
The mediation package above implements sequential ignorability — a specific set of identifying assumptions. The causal mediation chapter covers the assumptions more formally, introduces the controlled direct effect (CDE) and natural direct/indirect effects (NDE/NIE) using potential-outcome notation, and discusses identification with unmeasured mediator–outcome confounding. If your setting involves a non-binary outcome, treatment–mediator interaction, or you need a doubly-robust estimator, start there.
---title: "Mediation analysis in R and Stata"date: "2019-08-06"---## Mediation analysisTraditionally mediation model can be represented in the following equations:$$ Y = a X + b M + \epsilon_1 $$$$ M = c X + \epsilon_2 $$That is, we'd like to study the effect of $X$ on $Y$, and we see the effect can be a direct effect, and an indirect effect, through $M$.Baron and Kenny's (http://davidakenny.net/cm/mediate.htm) method is done in four steps. Modern approach tends to use SEM (structural equation modeling) to model these two equations directly.R's lavaan and Stata's sem commands are powerful tools.A simple example:```{r}#| echo: true#| message: falselibrary(lavaan)set.seed(1234)n <-10000X <-rnorm(n)M <-0.5*X +rnorm(n)Y <-0.7*M +0.3*X +rnorm(n)Data <-data.frame(X = X, Y = Y, M = M)model <-' Y ~ a*X + b*M M ~ c*X # direct effect (a) # indirect effect (b*c) bc := b*c # total effect total := a + (b*c) 'fit <-sem(model, data = Data)summary(fit)```Now let's do it in stata.```{r}#| echo: false#| message: falselibrary(Statamarkdown)``````{stata}*| cache: trueset obs 10000gen x= rnormal()gen m= .5*x + rnormal()gen y= .7*m + .3*x +rnormal()sem (y <- x m) (m <- x)estat teffects```The above examples should have direct effect of .3 and indirect effect of .35, and total effect of .65.## Causal Mediation analysisThe traditional mediation analysis has been criticized for the lack of causal interpretation. Without manipulation of the mediator, it is hard to interpret the effects causally, because even if the treatment is from random experiments, the mediator is often not. Therefore there could be an unmeasured confounder that is causing both $M$ and $Y$.R's "mediation" package is for causal mediation analysis. It uses simulation to estimate the causal effects of treatment, under assumptions of sequential ignorability.::: {.callout-note}## Identifying assumptions for natural (in)direct effectsBefore reading the `mediate()` output causally, the following must hold (this is "sequential ignorability" plus the standard support and link conditions):1. **Treatment ignorability** — no unmeasured treatment–outcome confounding given covariates $W$.2. **Mediator ignorability** — no unmeasured mediator–outcome confounding given $W$ and treatment.3. **No treatment-induced confounding** — no variable affected by the treatment confounds the mediator–outcome relationship (this is the strong, often-violated condition).4. **Positivity** — both treatment and mediator have positive probability across $W$.5. **Consistency** — observed outcomes equal the corresponding potential outcomes.6. **Cross-world independence** — the condition underlying natural effects, which no single experiment can test.Mediator ignorability and the no-treatment-induced-confounding condition are the ones that usually break, because the mediator is rarely randomized.:::It estimates the following quantities:$$\tau_i = Y_i(1, M_i(1)) - Y_i(0, M_i(0))$$This is the total treatment effect, which is to say, what's the change in $Y$ if we change each unit from control to treated, hypothetically?Then this can be decomposed into the causal mediation effects:$$\delta_i(t) = Y_i(t, M_i(1)) - Y_i(t, M_i(0))$$for treatment status $t=0,1$. This is to say, given treatment status, what's the mediation effect?$$\eta_i(t) = Y_i(1, M_i(t)) - Y_i(0, M_i(t)) $$for treatment status $t=0,1$. This is to say, given mediator status for each treatment status, what's the direct effect?Therefore there are four quantities estimated, direct and mediation effect for treated and control.R's "mediation" needs users to feed two models, outcome model and mediation model.If we study the same data, we would expect it returns the same estimates as the traditional methods. However, the causal mediation models can be much more flexible in outcome and mediation models.```{r}#| echo: true#| message: falselibrary(mediation)med.fit <-lm(M ~ X, data=Data)summary(med.fit)out.fit <-lm(Y ~ X + M,data=Data)summary(out.fit)set.seed(123)med.out <-mediate(med.fit, out.fit, treat ="X", mediator ="M", sims =100)summary(med.out)```## Binary outcomeFor example, in the case of binary outcome, the traditional approach will have difficulties. We can estimate the outcome model and mediator model jointly, but the total effects are not easy to decompose into direct and indirect effect (see Imai et al, page 320 https://imai.fas.harvard.edu/research/files/BaronKenny.pdf).The causal mediation analysis framework is much more general.```{r}#| echo: true#| message: falselibrary(mediation)data(framing)med.fit <-lm(emo ~ treat + age + educ + gender + income, data = framing)out.fit <-glm(cong_mesg ~ emo + treat + age + educ + gender + income,data = framing, family =binomial("probit"))set.seed(123)med.out <-mediate(med.fit, out.fit, treat ="treat", mediator ="emo",robustSE =TRUE, sims =100)summary(med.out)```"mediation" package has more functionalities, such as multilevel, interaction of treatment and mediator, etc.Stata's sem and gsem commands can model different situations, but the direct effect and indirect effects are not easy to compute, especially when you have binary outcome, or other non-continuous outcome situations. They are not designed for causal mediation analysis.## Going further: formal causal mediationThe `mediation` package above implements sequential ignorability — a specificset of identifying assumptions. The [causal mediation chapter](causal-mediation.qmd)covers the assumptions more formally, introduces the controlled direct effect(CDE) and natural direct/indirect effects (NDE/NIE) using potential-outcomenotation, and discusses identification with unmeasured mediator–outcomeconfounding. If your setting involves a non-binary outcome, treatment–mediatorinteraction, or you need a doubly-robust estimator, start there.---<!-- see-also-footer -->*Systematic treatment: [R](https://xiangao.github.io/causal_econometrics_guide/mediation.html) · [Julia](https://xiangao.github.io/causal_econometrics_julia/mediation.html).*