Traditionally mediation model can be represented in the following equestions:
\[ 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
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 = -43667.84
Iteration 1: Log likelihood = -43667.84
Structural equation model Number of obs = 10,000
Estimation method: ml
Log likelihood = -43667.84
------------------------------------------------------------------------------
| OIM
| Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
Structural |
y |
m | .7186314 .0113407 63.37 0.000 .6964041 .7408587
x | .3054332 .0127003 24.05 0.000 .2805411 .3303253
_cons | .199257 .0112834 17.66 0.000 .1771418 .2213721
-----------+----------------------------------------------------------------
m |
x | .5038605 .0100014 50.38 0.000 .4842581 .5234628
_cons | .0187204 .0099478 1.88 0.060 -.0007769 .0382177
-------------+----------------------------------------------------------------
var(e.y)| 1.272711 .0179989 1.237919 1.308482
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 | .7186314 .0113407 63.37 0.000 .6964041 .7408587
x | .3054332 .0127003 24.05 0.000 .2805411 .3303253
-----------+----------------------------------------------------------------
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 | .36209 .009182 39.43 0.000 .3440937 .3800863
-----------+----------------------------------------------------------------
m |
x | 0 (no path)
------------------------------------------------------------------------------
Total effects
------------------------------------------------------------------------------
| OIM
| Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
Structural |
y |
m | .7186314 .0113407 63.37 0.000 .6964041 .7408587
x | .6675231 .0134277 49.71 0.000 .6412053 .693841
-----------+----------------------------------------------------------------
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. It estimates the following quantities:
\[\tau_i = Y_i(1, M_i(1)) - Y_i(0, Mi(0))\]
This is the total treatment effefct, 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 stauts 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 tranditional 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.
“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.
Source Code
---title: "Mediation analysis in R and Stata"date: "2019-08-06"---## Mediation analysisTraditionally mediation model can be represented in the following equestions:$$ 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 + .4*m*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. It estimates the following quantities:$$\tau_i = Y_i(1, M_i(1)) - Y_i(0, Mi(0))$$This is the total treatment effefct, 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 stauts 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 tranditional 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.