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.

Advertisements
Multiple-mediator analysis with lavaan

14 thoughts on “Multiple-mediator analysis with lavaan

  1. Patricia says:

    Your blog is super helpful, thank you for the contents. Could you explain a little bit more why a one must add the covariance between the two moderators (M1 ~~ M2)? Could you suggest a paper that explores this? Thank you.

    Like

    1. Hi Patricia, I do not know of a paper exploring this topic in specific. Preacher and Hayes (2008) talk about this briefly, maybe you can find more references in their paper. That is also the paper I used to translate the mPlus formula to lavaan, see the `with’ clause in the mplus statements in appendix A.

      Like

    1. Yes it is. You will have to extend the model to accommodate for a third mediator. Starting from modeling the mediator itself, adding indirect and total effects plus modeling the covariances of the 3 mediators which would be m1~~m2, m1~~m3 and m2~~m3. I suggest you then check with semPlot whether you created all the paths… Goodluck!

      Like

  2. LEMAITRE says:

    Thanks for you help. I did what you recommended to do and it worked.
    Now, I’m little bit confused. For practicing, I changed the order of the two mediators in your example and I noticed that the “total effect” is not the same. Is there a reason for that?

    Defined Parameters:
    Estimate Std.Err z-value P(>|z|)
    indirect1 -0.237 0.087 -2.720 0.007
    indirect2 -0.287 0.088 -3.257 0.001
    total -0.328 0.206 -1.596 0.110

    Defined Parameters:
    Estimate Std.Err z-value P(>|z|)
    indirect1 -0.287 0.088 -3.257 0.001
    indirect2 -0.237 0.087 -2.720 0.007
    total -0.429 0.200 -2.139 0.032

    Liked by 1 person

    1. I could not reproduce your data. How did you invert the mediators? Using this formula:
      Y ~ b2 * M2 + b1 * M1 + c * X
      I get the exact same indirect and total effects for both models. If you like you could mail me your model specification [varati at hotmail dot com] and I can see whether I figure out what is different between our approaches. By the way. Your comment point out to me that I made a mistake in the specification of the total effect which I did not spot before.

      total := c + (a1 * b1) + (a1 * b1)
      should be:
      total := c + (a1 * b1) + (a2 * b2)

      Thank you very much for sending me into the right direction!

      Like

    1. That explains the difference. The reason you got different total effect is the mistake in the formula of the total effect

      total := c + (a1 * b1) + (a1 * b1)

      should be

      total := c + (a1 * b1) + (a2 * b2)

      If you try the proper formula the total effects are identical. The reason why it happened is that, the wrong formula was using the first mediator for both estimates. Therefore, inverting the order of the mediators yielded different estimates.

      Like

  3. Bogdan says:

    When I write this syntax for the model, I receive this error:

    lavaan ERROR: unknown label(s) in variable definition(s): a1 b1 a2 a3 a4 b2 a5 a6 b6 c1 c2 c3

    I’m not using R for a long time, so this could be a beginner’s question. Thank you for posting this! It is really helpful.

    Like

    1. Hi Bogdan, could you please provide more context? In my example I use those as coefficients which then estimated by lavaan. Because you do not provide enough context I cannot know whether those are you variables’ names or coefficients which should be estimated…

      Like

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s