SEM model with multiple mediators and multiple independent variables in lavaan

Question

I have cross-sectional data and I am trying to specify a model with multiple mediations.

My independent variable (IV) is measured by a tool with 24 items, which make up 5 subscales (latent variables), which in turn load onto a total "higher-order" factor. I try to estimate this model in two different ways: (1) A five-factor model (without a higher order factor) in which all 5 subscales are allowed to correlate and (2) a higher-order model with a TOTAL latent variable made up of those 5 suscales.The first model has five correlated latents (FNR, FOB...FAA) factors, with variance fixed to 1.

My first model has an IV (FTOTAL), 4 mediators (ER 1, 2, 3 and 4), and one DV (PH). The relationship between ER 4 and the IV is mediated by one of the other mediators (ER3), making it a mediated mediation. I was able to specify the SEM for the first model with the higher order TOTAL factor without a problem, and the code is shown below:

"#Measurements model
  FNR =~ FNR1 + FNR2 + FNR3 +FNR4 +FNR5
  
  FOB =~ FOB1 + FOB2 +FOB3 +FOB4
  
  FDS =~ FDS1 +FDS2 +FDS3 + FDS4 + FDS5
  
  FNJ =~ FNJ1 + FNJ2 + FNJ3 +FNJ4 + FNJ5
  
  FAA =~ FAA1 + FAA2 +FAA3 + FAA4 +FAA5
  
  FTOTAL =~ FNR + FOB + FDS + FNJ+ FAA 

#Regressions

ER3~ a*FTOTAL
ER4~ b*RSTOTAL +FTOTAL
ER1 ~ u1*FTOTAL
ER2 ~ u2*FTOTAL

PHQTOTAL ~ z1*ER1 + z2*ER2 + d*FTOTAL + c*ER4 + ER3

indirect1 := u1 * z1
indirect2 := u2 * z2
indirect3 := a*b*c
total := d + (u1 * z1) + (u2 * z2) + a*b*c


#Residual correlations
CRTOTAL~~SUPTOTAL+SLEEPTOTAL+RSTOTAL
SUPTOTAL~~SLEEPTOTAL+RSTOTAL

"

fitPHtotal <- sem(model = multipleMediationPH, data = SEMDATA, std.lv=TRUE)
summary(fitPHtotal)

However, I cannot figure out how to specify the model with all the 5 subscales as independent variables in the same model. I tried following the same logic and but including the 5 IVs in the model and it did not work out. Could anyone suggest a solution? or is the only way to just run 5 different models with one subscale as an IV at a time?

Thank you in advance for the help.

score 2 · Accepted Answer · answered Mar 30 '21 at 16:51

Since you did not provide data. I show you how to do it using the HolzingerSwineford1939 which comes with the lavaan package.

First, a mediation using a second-order latent factor (3 first-order factors):

library(lavaan)
#> This is lavaan 0.6-8
#> lavaan is FREE software! Please report any bugs.
model_2L <- "
visual  =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed   =~ x7 + x8 + x9
higher =~ visual + textual + speed

#grade will be your Y
#higher order latent factor will be your X
#agemo will be your M
grade ~ c*higher + b*agemo
agemo ~ a*higher

# indirect effect (a*b)
ab := a*b
# total effect
total := c + (a*b)
"

fit_2L <- sem(model = model_2L, data = HolzingerSwineford1939)
summary(object = fit_2L, std=T)
#> lavaan 0.6-8 ended normally after 48 iterations
#> 
#>   Estimator                                         ML
#>   Optimization method                           NLMINB
#>   Number of model parameters                        26
#>                                                       
#>                                                   Used       Total
#>   Number of observations                           300         301
#>                                                                   
#> Model Test User Model:
#>                                                       
#>   Test statistic                               116.110
#>   Degrees of freedom                                40
#>   P-value (Chi-square)                           0.000
#> 
#> Parameter Estimates:
#> 
#>   Standard errors                             Standard
#>   Information                                 Expected
#>   Information saturated (h1) model          Structured
#> 
#> Latent Variables:
#>                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
#>   visual =~                                                             
#>     x1                1.000                               0.849    0.727
#>     x2                0.621    0.109    5.680    0.000    0.527    0.448
#>     x3                0.824    0.124    6.641    0.000    0.699    0.619
#>   textual =~                                                            
#>     x4                1.000                               0.990    0.851
#>     x5                1.117    0.066   16.998    0.000    1.106    0.859
#>     x6                0.922    0.056   16.563    0.000    0.913    0.834
#>   speed =~                                                              
#>     x7                1.000                               0.648    0.595
#>     x8                1.130    0.148    7.612    0.000    0.732    0.726
#>     x9                1.010    0.135    7.465    0.000    0.655    0.649
#>   higher =~                                                             
#>     visual            1.000                               0.673    0.673
#>     textual           0.849    0.185    4.586    0.000    0.490    0.490
#>     speed             0.810    0.179    4.519    0.000    0.714    0.714
#> 
#> Regressions:
#>                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
#>   grade ~                                                               
#>     higher     (c)    0.421    0.089    4.730    0.000    0.241    0.482
#>     agemo      (b)   -0.004    0.008   -0.519    0.604   -0.004   -0.029
#>   agemo ~                                                               
#>     higher     (a)    0.322    0.469    0.687    0.492    0.184    0.053
#> 
#> Variances:
#>                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
#>    .x1                0.641    0.110    5.822    0.000    0.641    0.471
#>    .x2                1.108    0.102   10.848    0.000    1.108    0.799
#>    .x3                0.786    0.094    8.398    0.000    0.786    0.616
#>    .x4                0.373    0.048    7.750    0.000    0.373    0.276
#>    .x5                0.436    0.058    7.453    0.000    0.436    0.263
#>    .x6                0.364    0.044    8.369    0.000    0.364    0.304
#>    .x7                0.767    0.080    9.629    0.000    0.767    0.646
#>    .x8                0.482    0.070    6.924    0.000    0.482    0.474
#>    .x9                0.589    0.068    8.686    0.000    0.589    0.579
#>    .grade             0.192    0.020    9.767    0.000    0.192    0.768
#>    .agemo            11.881    0.972   12.220    0.000   11.881    0.997
#>    .visual            0.394    0.111    3.535    0.000    0.547    0.547
#>    .textual           0.745    0.101    7.397    0.000    0.760    0.760
#>    .speed             0.206    0.062    3.312    0.001    0.490    0.490
#>     higher            0.327    0.097    3.375    0.001    1.000    1.000
#> 
#> Defined Parameters:
#>                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
#>     ab               -0.001    0.004   -0.366    0.715   -0.001   -0.002
#>     total             0.420    0.089    4.728    0.000    0.240    0.481

Second, a mediation using a three first-order factors. Three indirect effects and three total effects are estimated:

library(lavaan)
model_1L <- "
visual  =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed   =~ x7 + x8 + x9

#grade will be your Y
#higher order latent factor will be your X
#agemo will be your M
grade ~ c1*visual + c2*textual + c3*speed + b*agemo
agemo ~ a1*visual + a2*textual + a3*speed

# indirect effect (a*b)
a1b := a1*b
a2b := a2*b
a3b := a3*b

# total effect
total1 := c1 + (a1*b)
total2 := c2 + (a2*b)
total3 := c3 + (a3*b)
"

fit_1L <- sem(model = model_1L, data = HolzingerSwineford1939)
summary(object = fit_1L, std=T)
#> lavaan 0.6-8 ended normally after 55 iterations
#> 
#>   Estimator                                         ML
#>   Optimization method                           NLMINB
#>   Number of model parameters                        30
#>                                                       
#>                                                   Used       Total
#>   Number of observations                           300         301
#>                                                                   
#> Model Test User Model:
#>                                                       
#>   Test statistic                               101.925
#>   Degrees of freedom                                36
#>   P-value (Chi-square)                           0.000
#> 
#> Parameter Estimates:
#> 
#>   Standard errors                             Standard
#>   Information                                 Expected
#>   Information saturated (h1) model          Structured
#> 
#> Latent Variables:
#>                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
#>   visual =~                                                             
#>     x1                1.000                               0.904    0.775
#>     x2                0.555    0.100    5.564    0.000    0.501    0.426
#>     x3                0.724    0.109    6.657    0.000    0.655    0.580
#>   textual =~                                                            
#>     x4                1.000                               0.993    0.853
#>     x5                1.108    0.065   17.017    0.000    1.101    0.855
#>     x6                0.921    0.055   16.667    0.000    0.915    0.836
#>   speed =~                                                              
#>     x7                1.000                               0.668    0.613
#>     x8                1.115    0.142    7.840    0.000    0.744    0.737
#>     x9                0.945    0.125    7.540    0.000    0.631    0.625
#> 
#> Regressions:
#>                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
#>   grade ~                                                               
#>     visual    (c1)    0.012    0.048    0.246    0.806    0.011    0.021
#>     textual   (c2)    0.048    0.035    1.376    0.169    0.047    0.095
#>     speed     (c3)    0.295    0.063    4.689    0.000    0.197    0.394
#>     agemo      (b)   -0.003    0.008   -0.361    0.718   -0.003   -0.020
#>   agemo ~                                                               
#>     visual    (a1)    0.354    0.355    0.996    0.319    0.320    0.093
#>     textual   (a2)   -0.233    0.256   -0.912    0.362   -0.231   -0.067
#>     speed     (a3)    0.098    0.421    0.232    0.817    0.065    0.019
#> 
#> Covariances:
#>                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
#>   visual ~~                                                             
#>     textual           0.412    0.074    5.565    0.000    0.459    0.459
#>     speed             0.265    0.058    4.554    0.000    0.438    0.438
#>   textual ~~                                                            
#>     speed             0.180    0.052    3.448    0.001    0.271    0.271
#> 
#> Variances:
#>                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
#>    .x1                0.545    0.115    4.747    0.000    0.545    0.400
#>    .x2                1.135    0.102   11.115    0.000    1.135    0.819
#>    .x3                0.846    0.091    9.322    0.000    0.846    0.664
#>    .x4                0.368    0.048    7.698    0.000    0.368    0.272
#>    .x5                0.447    0.058    7.657    0.000    0.447    0.270
#>    .x6                0.361    0.043    8.343    0.000    0.361    0.301
#>    .x7                0.741    0.079    9.422    0.000    0.741    0.624
#>    .x8                0.465    0.069    6.724    0.000    0.465    0.456
#>    .x9                0.620    0.067    9.217    0.000    0.620    0.609
#>    .grade             0.201    0.018   11.307    0.000    0.201    0.806
#>    .agemo            11.813    0.969   12.191    0.000   11.813    0.991
#>     visual            0.817    0.147    5.564    0.000    1.000    1.000
#>     textual           0.986    0.113    8.752    0.000    1.000    1.000
#>     speed             0.446    0.091    4.906    0.000    1.000    1.000
#> 
#> Defined Parameters:
#>                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
#>     a1b              -0.001    0.003   -0.344    0.731   -0.001   -0.002
#>     a2b               0.001    0.002    0.335    0.738    0.001    0.001
#>     a3b              -0.000    0.002   -0.183    0.855   -0.000   -0.000
#>     total1            0.011    0.048    0.226    0.821    0.010    0.020
#>     total2            0.048    0.035    1.399    0.162    0.048    0.096
#>     total3            0.295    0.063    4.685    0.000    0.197    0.394

^{Created on 2021-03-30 by the reprex package (v1.0.0)}

SEM model with multiple mediators and multiple independent variables in lavaan

1 Answers1