Multilevel Models with plssem

This vignette shows examples of multilevel random slopes and intercept models, with both continuous and ordinal data.

Random Slopes Model

slopes_model <- "
  X =~ x1 + x2 + x3
  Z =~ z1 + z2 + z3
  Y =~ y1 + y2 + y3
  W =~ w1 + w2 + w3
  Y ~ X + Z + (1 + X + Z | cluster)
  W ~ X + Z + (1 + X + Z | cluster)
"

Continuous Indicators

fit_slopes_cont <- pls(
  slopes_model,
  data      = randomSlopes,
  bootstrap = TRUE,
  sample    = 50
)
summary(fit_slopes_cont)
#> plssem (0.1.0) ended normally after 3 iterations
#> 
#>   Estimator                                   PLSc-MLM
#>   Link                                          LINEAR
#>                                                       
#>   Number of observations                          5000
#>   Number of iterations                               3
#>   Number of latent variables                         4
#>   Number of observed variables                      18
#> 
#> R-squared (indicators):
#>   x1                                             0.860
#>   x2                                             0.683
#>   x3                                             0.772
#>   z1                                             0.839
#>   z2                                             0.694
#>   z3                                             0.759
#>   y1                                             0.838
#>   y2                                             0.726
#>   y3                                             0.750
#>   w1                                             0.840
#>   w2                                             0.695
#>   w3                                             0.771
#> 
#> R-squared (latents):
#>   Y                                              0.533
#>   W                                              0.547
#> 
#> Latent Variables:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   X =~          
#>     x1              0.927      0.009   98.707    0.000
#>     x2              0.826      0.010   79.620    0.000
#>     x3              0.879      0.008  114.894    0.000
#>   Z =~          
#>     z1              0.916      0.010   87.450    0.000
#>     z2              0.833      0.012   71.560    0.000
#>     z3              0.871      0.009   94.854    0.000
#>   Y =~          
#>     y1              0.915      0.011   85.376    0.000
#>     y2              0.852      0.011   76.649    0.000
#>     y3              0.866      0.014   61.430    0.000
#>   W =~          
#>     w1              0.916      0.011   86.168    0.000
#>     w2              0.834      0.016   50.799    0.000
#>     w3              0.878      0.014   64.442    0.000
#> 
#> Regressions:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   Y ~           
#>     X               0.293      0.021   13.836    0.000
#>     Z               0.444      0.038   11.540    0.000
#>   W ~           
#>     X               0.393      0.033   11.946    0.000
#>     Z               0.248      0.042    5.868    0.000
#> 
#> Intercepts:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>    .Y               0.009      0.005    1.790    0.074
#>    .W               0.009      0.008    1.249    0.212
#> 
#> Covariances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   X ~~          
#>     Z               0.175      0.012   14.107    0.000
#>   Y~X ~~        
#>     Y~1            -0.003      0.006   -0.579    0.562
#>   Y~Z ~~        
#>     Y~1            -0.024      0.014   -1.769    0.077
#>     Y~X             0.012      0.009    1.422    0.155
#>   W~X ~~        
#>     W~1             0.003      0.012    0.256    0.798
#>   W~Z ~~        
#>     W~1             0.009      0.011    0.815    0.415
#>     W~X             0.013      0.012    1.112    0.266
#> 
#> Variances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>     X               1.000                             
#>     Z               1.000                             
#>    .Y               0.467      0.028   16.699    0.000
#>    .W               0.453      0.031   14.606    0.000
#>    .x1              0.140      0.017    8.039    0.000
#>    .x2              0.317      0.017   18.434    0.000
#>    .x3              0.228      0.013   16.940    0.000
#>    .z1              0.161      0.019    8.381    0.000
#>    .z2              0.306      0.019   15.729    0.000
#>    .z3              0.241      0.016   15.072    0.000
#>    .y1              0.162      0.020    8.280    0.000
#>    .y2              0.274      0.019   14.492    0.000
#>    .y3              0.250      0.024   10.343    0.000
#>    .w1              0.160      0.020    8.201    0.000
#>    .w2              0.305      0.027   11.250    0.000
#>    .w3              0.229      0.024    9.599    0.000
#>     Y~1             0.079      0.019    4.236    0.000
#>     Y~X             0.018      0.004    4.099    0.000
#>     Y~Z             0.105      0.017    6.253    0.000
#>     W~1             0.054      0.012    4.705    0.000
#>     W~X             0.095      0.015    6.326    0.000
#>     W~Z             0.144      0.026    5.483    0.000

Ordered Indicators

fit_slopes_ord <- pls(
  slopes_model,
  data      = randomSlopesOrdered,
  bootstrap = TRUE,
  sample    = 50,
  ordered   = colnames(randomSlopesOrdered) # explicitly specify variables as ordered
)
summary(fit_slopes_ord)
#> plssem (0.1.0) ended normally after 3 iterations
#> 
#>   Estimator                                OrdPLSc-MLM
#>   Link                                          PROBIT
#>                                                       
#>   Number of observations                          5000
#>   Number of iterations                               3
#>   Number of latent variables                         4
#>   Number of observed variables                      18
#> 
#> R-squared (indicators):
#>   x1                                             0.870
#>   x2                                             0.669
#>   x3                                             0.789
#>   z1                                             0.841
#>   z2                                             0.714
#>   z3                                             0.758
#>   y1                                             0.844
#>   y2                                             0.711
#>   y3                                             0.754
#>   w1                                             0.835
#>   w2                                             0.681
#>   w3                                             0.785
#> 
#> R-squared (latents):
#>   Y                                              0.531
#>   W                                              0.544
#> 
#> Latent Variables:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   X =~          
#>     x1              0.933      0.011   87.644    0.000
#>     x2              0.818      0.011   72.002    0.000
#>     x3              0.888      0.010   86.496    0.000
#>   Z =~          
#>     z1              0.917      0.009  100.749    0.000
#>     z2              0.845      0.013   66.353    0.000
#>     z3              0.871      0.010   83.421    0.000
#>   Y =~          
#>     y1              0.918      0.009   97.607    0.000
#>     y2              0.843      0.012   71.023    0.000
#>     y3              0.869      0.014   62.273    0.000
#>   W =~          
#>     w1              0.914      0.013   69.631    0.000
#>     w2              0.825      0.014   57.186    0.000
#>     w3              0.886      0.016   56.953    0.000
#> 
#> Regressions:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   Y ~           
#>     X               0.290      0.022   13.000    0.000
#>     Z               0.451      0.033   13.676    0.000
#>   W ~           
#>     X               0.391      0.033   11.720    0.000
#>     Z               0.245      0.052    4.729    0.000
#> 
#> Intercepts:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>    .Y               0.011      0.005    2.204    0.028
#>    .W               0.009      0.008    1.045    0.296
#> 
#> Covariances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   X ~~          
#>     Z               0.168      0.014   11.947    0.000
#>   Y~X ~~        
#>     Y~1            -0.006      0.005   -1.401    0.161
#>   Y~Z ~~        
#>     Y~1            -0.024      0.013   -1.825    0.068
#>     Y~X             0.012      0.009    1.400    0.161
#>   W~X ~~        
#>     W~1             0.004      0.011    0.335    0.738
#>   W~Z ~~        
#>     W~1             0.009      0.011    0.874    0.382
#>     W~X             0.016      0.013    1.305    0.192
#> 
#> Variances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>     X               1.000                             
#>     Z               1.000                             
#>    .Y               0.469      0.028   16.941    0.000
#>    .W               0.456      0.035   12.917    0.000
#>    .x1              0.130      0.020    6.552    0.000
#>    .x2              0.331      0.019   17.840    0.000
#>    .x3              0.211      0.018   11.629    0.000
#>    .z1              0.159      0.017    9.487    0.000
#>    .z2              0.286      0.021   13.325    0.000
#>    .z3              0.242      0.018   13.359    0.000
#>    .y1              0.156      0.017    9.075    0.000
#>    .y2              0.289      0.020   14.466    0.000
#>    .y3              0.246      0.024   10.183    0.000
#>    .w1              0.165      0.024    6.885    0.000
#>    .w2              0.319      0.024   13.424    0.000
#>    .w3              0.215      0.027    7.858    0.000
#>     Y~1             0.076      0.018    4.329    0.000
#>     Y~X             0.018      0.005    3.726    0.000
#>     Y~Z             0.101      0.014    7.370    0.000
#>     W~1             0.052      0.012    4.550    0.000
#>     W~X             0.099      0.016    6.141    0.000
#>     W~Z             0.142      0.027    5.336    0.000

Random Intercepts Model

intercepts_model <- '
  f =~ y1 + y2 + y3
  f ~ x1 + x2 + x3 + w1 + w2 + (1 | cluster)
'

Continuous Indicators

fit_intercepts_cont <- pls(
  intercepts_model,
  data      = randomIntercepts,
  bootstrap = TRUE,
  sample    = 50
)
summary(fit_intercepts_cont)
#> plssem (0.1.0) ended normally after 2 iterations
#> 
#>   Estimator                                   PLSc-MLM
#>   Link                                          LINEAR
#>                                                       
#>   Number of observations                         10000
#>   Number of iterations                               2
#>   Number of latent variables                         1
#>   Number of observed variables                       9
#> 
#> R-squared (indicators):
#>   y1                                             0.891
#>   y2                                             0.785
#>   y3                                             0.814
#> 
#> R-squared (latents):
#>   f                                              0.705
#> 
#> Latent Variables:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   f =~          
#>     y1              0.944      0.008  120.524    0.000
#>     y2              0.886      0.009   96.581    0.000
#>     y3              0.902      0.010   87.518    0.000
#> 
#> Regressions:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   f ~           
#>     x1              0.238      0.007   32.191    0.000
#>     x2              0.162      0.007   23.900    0.000
#>     x3              0.077      0.005   15.071    0.000
#>     w1              0.128      0.032    4.024    0.000
#>     w2              0.091      0.037    2.486    0.013
#> 
#> Intercepts:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>    .f              -0.012      0.008   -1.617    0.106
#> 
#> Covariances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   x1 ~~         
#>     x2              0.104      0.009   11.591    0.000
#>     x3              0.004      0.012    0.352    0.725
#>     w1              0.000                             
#>     w2              0.000                             
#>   x2 ~~         
#>     x3              0.097      0.013    7.568    0.000
#>     w1              0.000                             
#>     w2              0.000                             
#>   x3 ~~         
#>     w1              0.000                             
#>     w2              0.000                             
#>   w1 ~~         
#>     w2             -0.041      0.044   -0.923    0.356
#> 
#> Variances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>    .f               0.295      0.015   19.227    0.000
#>     x1              1.000                             
#>     x2              1.000                             
#>     x3              1.000                             
#>     w1              1.000                             
#>     w2              1.000                             
#>    .y1              0.109      0.015    7.343    0.000
#>    .y2              0.215      0.016   13.196    0.000
#>    .y3              0.186      0.019   10.018    0.000
#>     f~1             0.582      0.021   28.051    0.000

Ordered Indicators

fit_intercepts_ord <- pls(
  intercepts_model,
  data      = randomInterceptsOrdered,
  bootstrap = TRUE,
  sample    = 50,
  ordered   = colnames(randomInterceptsOrdered) # explicitly specify variables as ordered
)
summary(fit_intercepts_ord)
#> plssem (0.1.0) ended normally after 2 iterations
#> 
#>   Estimator                                OrdPLSc-MLM
#>   Link                                          PROBIT
#>                                                       
#>   Number of observations                         10000
#>   Number of iterations                               2
#>   Number of latent variables                         1
#>   Number of observed variables                       9
#> 
#> R-squared (indicators):
#>   y1                                             0.885
#>   y2                                             0.788
#>   y3                                             0.809
#> 
#> R-squared (latents):
#>   f                                              0.684
#> 
#> Latent Variables:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   f =~          
#>     y1              0.941      0.011   83.225    0.000
#>     y2              0.888      0.014   64.791    0.000
#>     y3              0.899      0.012   73.257    0.000
#> 
#> Regressions:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   f ~           
#>     x1              0.241      0.008   30.016    0.000
#>     x2              0.158      0.009   18.199    0.000
#>     x3              0.080      0.007   12.154    0.000
#>     w1              0.121      0.049    2.480    0.013
#>     w2              0.081      0.044    1.861    0.063
#> 
#> Intercepts:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>    .f              -0.011      0.005   -2.184    0.029
#> 
#> Covariances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   x1 ~~         
#>     x2              0.110      0.011    9.693    0.000
#>     x3              0.012      0.011    1.124    0.261
#>     w1              0.001      0.004    0.186    0.853
#>     w2              0.001      0.004    0.279    0.781
#>   x2 ~~         
#>     x3              0.099      0.012    7.955    0.000
#>     w1             -0.003      0.003   -1.097    0.273
#>     w2             -0.001      0.003   -0.464    0.643
#>   x3 ~~         
#>     w1             -0.003      0.003   -0.873    0.382
#>     w2              0.003      0.003    0.961    0.337
#>   w1 ~~         
#>     w2             -0.025      0.058   -0.430    0.667
#> 
#> Variances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>    .f               0.316      0.014   22.558    0.000
#>     x1              1.000                             
#>     x2              1.000                             
#>     x3              1.000                             
#>     w1              1.000                             
#>     w2              1.000                             
#>    .y1              0.115      0.021    5.405    0.000
#>    .y2              0.212      0.024    8.651    0.000
#>    .y3              0.191      0.022    8.696    0.000
#>     f~1             0.562      0.021   27.261    0.000