How to compare mediation models with one or two mediators?

Question

I'm looking to compare the indirect effect sizes of two models:

IV1 -> IV2 -> DV

and

IV1 -> IV2 -> IV3 -> DV

Is there a way to do this in R?

I know how to run the first model using the mediation package. But I'm not sure how to have multiple mediators.

Answer 1

You can easily do that with the package lavaan. Here are the steps :

1. you create the second mediation model;
2. you fix the mediator 2's parameters to 0;
3. you run both models;
4. you can compare them with the anova() function.

Here is an example with an hypothetical data set.

# A simple data set
set.seed(42)
jd <- MASS::mvrnorm(120,
                    mu = rep(0, 5),
                    Sigma = matrix(c(  1, .5,-.2,-.1,
                                      .5,  1,-.4,-.2,
                                     -.2,-.4,  1, .5,
                                     -.1,-.2, .5,  1,), 
                                   4, 4),
                    empirical = TRUE)
colnames(jd) <- c("X","M1","M2","Y")
jd <- as.data.frame(jd)

# The two models
# 1) The second model
model2 <- ' # direct effect
             Y ~  c * X
            # mediator1
             M1 ~ a * X
             Y  ~ b * M1
            # mediator2
             M2 ~ d * X 
             M2 ~ e * M1
             Y  ~ f * M2
            # indirect effects
             ab  := a * b
             df  := d * f
             aef := a * e * f
            # total effect
             total := c + (a*b + d*f + a*e*f)
         '
# 2) Fix parameters to 0
model1 <- ' # direct effect
             Y ~  c * X
            # mediator1
             M1 ~ a * X
             Y  ~ b * M1
            # mediator2
             M2 ~ 0 * X 
             M2 ~ 0 * M1
             Y  ~ 0 * M2
            # indirect effects
             ab  := a * b
            # total effect
             total := c + (a*b)
         '
# 3) Carry both models
library(lavaan)
mod1 <- sem(model = model1, data = jd)
mod2 <- sem(model = model2, data = jd)
# summary(mod1)
# summary(mod2)

# 4) Compare
anova(mod2, mod1)

The output :

Chi-Squared Difference Test

     Df    AIC     BIC  Chisq Chisq diff   RMSEA Df diff Pr(>Chisq)    
mod2  0 946.66  971.74  0.000                                          
mod1  3 991.20 1007.93 50.546     50.546 0.36342       3  6.114e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The Pr(>Chisq) column tells you that the model are significantly different and the one with the lower AIC is the best.

Make sure you have the same variables and a different number of free parameters, i.e., why we added the second mediator in the simple mediation model.