Multiple-mediator analysis with lavaan

I wrote this brief introductory post for my friend Simon. I want to show how easy the transition from SPSS to R can be. In the specific case of mediation analysis the transition to R can be very smooth because, thanks to lavaan, the R knowledge required to use the package is minimal. Analysis of mediator effects in lavaan requires only the specification of the model, all the other processes are automated by the package. So, after reading in the data, running the test is trivial.

This time, to keep the focus on the mediation analysis I will skip reading-in the data and generate a synthetic dataset instead. This is because otherwise I would have to spend the next paragraph explaining the dataset and the variables it contains and I really want to only focus on the analysis.

Create toy data:

# make up some data with one mediator
set.seed(06052017)
X <- rnorm(100)
M <- 0.5*X + rnorm(100)
Y <- 0.7*M + rnorm(100)
Data <- data.frame(X = X, Y = Y, M = M)

# add another mediator
M2 <- -0.35*X + rnorm(100)
Y <- 0.7*M2 + 0.48*-M + rnorm(100)
Data <- data.frame(X = X, Y = Y, M1 = M, M2 = M2)

As shown in the lavaan website performing a mediation analysis is as simple as typing in the code below:

simpleMediation <- '
	Y ~ b * M + c * X
	M ~ a * X
	indirect := a * b
	total    := c + (a * b)
'
require(lavaan)
fit <- sem(model = simpleMediation, data = Data)
summary(fit)

For multiple mediators one simply need to extend the model recycling the code of the first mediator variable:

multipleMediation <- '
Y ~ b1 * M1 + b2 * M2 + c * X
M1 ~ a1 * X
M2 ~ a2 * X
indirect1 := a1 * b1
indirect2 := a2 * b2
total := c + (a1 * b1) + (a2 * b2)
M1 ~~ M2
'
fit <- sem(model = multipleMediation, data = Data)
summary(fit)

Note that with multiple mediators we must add the covariance of the two mediators to the model. Covariances are added using the notation below:

M1 ~~ M2

There are two ways to test the null hypothesis that the indirect effect are equal to each other. The first is to specify a contrast for the two indirect effects. In the definition of the contrast the two indirect effects are subtracted. If it is significant the two indirect effects differ.

contrast := indirect1 – indirect2

# including a contrast in the model
contrastsMediation <- '
Y ~ b1 * M1 + b2 * M2 + c * X
M1 ~ a1 * X
M2 ~ a2 * X
indirect1 := a1 * b1
indirect2 := a2 * b2
contrast := indirect1 – indirect2
total := c + (a1 * b1) + (a2 * b2)
M1 ~~ M2
'
fit <- sem(model = contrastsMediation, data = Data)
summary(fit)

The second option to determine whether the indirect effects differ is to set a constrain in the model specifying the two indirect effect to be equal.

indirect1 == indirect2

Then, with the anova function one can compare the models and determine which one is better. Including the constrain and comparing the models is simple:

constrainedMediation <- '
Y ~ b1 * M1 + b2 * M2 + c * X
M1 ~ a1 * X
M2 ~ a2 * X
indirect1 := a1 * b1
indirect2 := a2 * b2
total := c + (a1 * b1) + (a2 * b2)
# covariances
M1 ~~ M2
# constrain
indirect1 == indirect2
'
noConstrFit <- sem(model = multipleMediation, data = Data)
constrFit <- sem(model = constrainedMediation, data = Data)
anova(noConstrFit, constrFit)

In my case the test is not significant so there is no evidence that the indirect effects are different.

For these toy models there is no further need of customizing the calls to sem. However, when performing a proper analysis one might prefer to have bootstrapped confidence intervals. The bootstrap standard errors and the quantity of bootstrap samples to perform are controlled in the call to the sem function with the additional arguments se' andbootstrap’, as in:

fit <- sem(
	model = contrastsMediation,
	data  = Data,
 	se = "bootstrap",
 	bootstrap = 5000 # 1000 is the default
)

(Bootstrap) confidence interval can be extracted with the function calls 1) summary, 2) parameterEstimates, or 3) bootstrapLavaan

summary(fit, fit.measures=TRUE, standardize=TRUE, rsquare=TRUE,
estimates = TRUE, ci = TRUE)

parameterEstimates(fit, boot.ci.type="bca.simple")

bootstrapLavaan(fit, R=5000, type="bollen.stine",
FUN=fitMeasures, fit.measures="chisq")

NOTE that bootstrapLavaan will re-compute the bootstrap samples requiring to wait as long as it took the sem function to run if called with the bootstrap option.

Since this post is longer than I wanted it to be, I will leave as a brief introduction to mediation with lavaan. In this follow up post I describe multiple mediation with lavaan using an actual dataset. And if you might be interested in how to automate the plotting of lavaan’s output into a pretty graph read this post. On github is the whole code in one .R file.

—- UPDATE Feb 25th 2020 —-
I wrote a few posts that follow up on this from a structural equation modelling perspective which I believe can clarify some of the issues raised in the comments below. Here is the first post.

For customization of the plots which can be created using semPlot see this post too.