How to automate mediation analysis in R using structural equation modelling structure?

Question

Suppose I have a data set df containing data on 2 media variables Fb and TV and Sales as the dependent variable. Now in order to investigate the indirect effect of Facebook on Sales via TV which acts as a mediator variable the following R code is used.

specmod <- "

#Path c (direct effect)

Sales ~ c*Facebook

#Path a

TV ~ a*Facebook

#Path b

Sales ~ b*TV

#Indirect effect (a*b)

ab := a*b

"

set.seed(1214)

fitmod <- sem(specmod, data = df, se = "bootstrap", bootstrap = 100)

parameterEstimates(fitmod, ci = TRUE, level = 0.95, boot.ci.type = "perc")

But if I have 9 media variables and want to check the mediation effect of all the variables then how should I proceed? The predictor and the mediator variable for each model iteration will be different combinations of the media variables. How can I automate the above code with different model iterations.

Answer 1

This is not an elegant solution but it is simple in terms of programming:

set.seed(581745)
dat <- MASS::mvrnorm(100, rep(0, 7), diag(7))
colnames(dat) <- c("y", paste0("m", 1:3), paste0("x", 1:3))
head(dat)
#>                y          m1         m2         m3          x1          x2
#> [1,]  0.19461669 -0.64146102 -0.7205590  0.7927056 -0.28878053  2.43030745
#> [2,] -1.04997220  1.48515283 -1.2392574  0.3022706 -0.12727050 -0.69668095
#> [3,]  0.31000671 -0.11000903  0.3791038 -1.1362097 -0.39751624 -0.32798805
#> [4,]  0.03023118  0.02059596  0.6318862 -0.7585338 -0.04914222  1.52437138
#> [5,]  0.20720111 -0.28528315 -0.5581268  0.3280841  0.68032932  0.48519911
#> [6,] -1.81659430 -0.11340674  0.3485680  0.3389938  1.03419571  0.09771754
#>               x3
#> [1,] -0.83107742
#> [2,] -1.64153599
#> [3,] -1.17252108
#> [4,]  1.11641882
#> [5,]  0.33163842
#> [6,]  0.08189787

x <- c("x1", "x2", "x3")
m <- c("m1", "m2", "m3")

mods <- paste0(
"# Path c (direct effect)
y ~ c*", x, "\n",
"# Path a", "\n",
m, " ~ a*", x, "\n",
"# Path b
y ~ b*", m, "\n",
"#Indirect effect (a*b)
ab := a*b")

cat(mods[1])
#> # Path c (direct effect)
#> y ~ c*x1
#> # Path a
#> m1 ~ a*x1
#> # Path b
#> y ~ b*m1
#> #Indirect effect (a*b)
#> ab := a*b
cat(mods[2])
#> # Path c (direct effect)
#> y ~ c*x2
#> # Path a
#> m2 ~ a*x2
#> # Path b
#> y ~ b*m2
#> #Indirect effect (a*b)
#> ab := a*b
cat(mods[3])
#> # Path c (direct effect)
#> y ~ c*x3
#> # Path a
#> m3 ~ a*x3
#> # Path b
#> y ~ b*m3
#> #Indirect effect (a*b)
#> ab := a*b

library(lavaan)
#> This is lavaan 0.6-15
#> lavaan is FREE software! Please report any bugs.
fits <- lapply(mods, sem, data = dat, se = "boot", iseed = 1234, bootstrap = 50)
lapply(fits, parameterEstimates)
#> [[1]]
#>   lhs op rhs label    est    se      z pvalue ci.lower ci.upper
#> 1   y  ~  x1     c  0.170 0.081  2.111  0.035   -0.012    0.336
#> 2  m1  ~  x1     a  0.072 0.086  0.839  0.402   -0.133    0.226
#> 3   y  ~  m1     b -0.121 0.087 -1.387  0.165   -0.310    0.067
#> 4   y ~~   y        0.974 0.119  8.187  0.000    0.703    1.185
#> 5  m1 ~~  m1        0.862 0.114  7.555  0.000    0.609    1.099
#> 6  x1 ~~  x1        0.980 0.000     NA     NA    0.980    0.980
#> 7  ab := a*b    ab -0.009 0.017 -0.524  0.600   -0.056    0.027
#> 
#> [[2]]
#>   lhs op rhs label   est    se     z pvalue ci.lower ci.upper
#> 1   y  ~  x2     c 0.029 0.090 0.319  0.750   -0.131    0.296
#> 2  m2  ~  x2     a 0.101 0.108 0.932  0.352   -0.084    0.385
#> 3   y  ~  m2     b 0.033 0.084 0.389  0.697   -0.117    0.244
#> 4   y ~~   y       1.010 0.117 8.656  0.000    0.728    1.228
#> 5  m2 ~~  m2       1.198 0.214 5.592  0.000    0.840    1.732
#> 6  x2 ~~  x2       1.103 0.000    NA     NA    1.103    1.103
#> 7  ab := a*b    ab 0.003 0.019 0.172  0.864   -0.019    0.084
#> 
#> [[3]]
#>   lhs op rhs label    est    se      z pvalue ci.lower ci.upper
#> 1   y  ~  x3     c -0.065 0.097 -0.675  0.499   -0.240    0.119
#> 2  m3  ~  x3     a -0.054 0.097 -0.558  0.577   -0.279    0.182
#> 3   y  ~  m3     b -0.177 0.099 -1.779  0.075   -0.378    0.062
#> 4   y ~~   y        0.986 0.121  8.150  0.000    0.720    1.245
#> 5  m3 ~~  m3        0.718 0.112  6.420  0.000    0.492    0.929
#> 6  x3 ~~  x3        1.118 0.000     NA     NA    1.118    1.118
#> 7  ab := a*b    ab  0.010 0.018  0.534  0.594   -0.041    0.052

First, form the vectors of x-variables and m-variables. Second, use paste() to create a vector of the model syntax strings. The call to cat() is optional, for verifying the generated model syntax strings. Third, simply use lapply() to fit each of the generated models. If bootstrap CIs are needed, need to set iseed to ensure that the same set of bootstrap samples are used in all models.

lavaan has semList() for fitting a model to a list of datasets. However, I am not aware of a similar function to fit a list of models to the same dataset. If yes, that function can be used in place of lapply().