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,
  boot.R    = 50
)
summary(fit_slopes_cont)
#> plssem (0.1.1) ended normally after 3 iterations
#> 
#>   Estimator                                   PLSc-MLM
#>   Link                                          PROBIT
#>                                                       
#>   Number of observations                          5000
#>   Number of iterations                               3
#>   Number of latent variables                         4
#>   Number of observed variables                      18
#> 
#> Fit Measures:
#>   Chi-Square                                   127.210
#>   Degrees of Freedom                                49
#>   SRMR                                           0.008
#>   RMSEA                                          0.018
#> 
#> 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.542
#>   W                                              0.553
#> 
#> Latent Variables:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   X =~          
#>     x1              0.927      0.010   96.515    0.000
#>     x2              0.826      0.010   79.841    0.000
#>     x3              0.879      0.007  118.521    0.000
#>   Z =~          
#>     z1              0.916      0.009  106.571    0.000
#>     z2              0.833      0.011   74.956    0.000
#>     z3              0.871      0.009   97.399    0.000
#>   Y =~          
#>     y1              0.915      0.009  105.080    0.000
#>     y2              0.852      0.013   65.547    0.000
#>     y3              0.866      0.013   67.926    0.000
#>   W =~          
#>     w1              0.916      0.011   84.423    0.000
#>     w2              0.834      0.013   62.526    0.000
#>     w3              0.878      0.011   83.457    0.000
#> 
#> Regressions:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   Y ~           
#>     X               0.293      0.022   13.346    0.000
#>     Z               0.444      0.041   10.747    0.000
#>   W ~           
#>     X               0.393      0.030   13.077    0.000
#>     Z               0.248      0.054    4.591    0.000
#> 
#> Intercepts:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>    .Y               0.010      0.005    2.149    0.032
#>    .W               0.010      0.008    1.312    0.190
#> 
#> Covariances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   X ~~          
#>     Z               0.175      0.013   13.799    0.000
#>   Y~X ~~        
#>     Y~1            -0.004      0.005   -0.707    0.480
#>   Y~Z ~~        
#>     Y~1            -0.025      0.011   -2.337    0.019
#>     Y~X             0.012      0.009    1.413    0.158
#>   W~X ~~        
#>     W~1             0.003      0.012    0.263    0.793
#>   W~Z ~~        
#>     W~1             0.009      0.014    0.645    0.519
#>     W~X             0.013      0.012    1.120    0.263
#> 
#> Variances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>     X               1.000                             
#>     Z               1.000                             
#>    .Y               0.458      0.028   16.444    0.000
#>    .W               0.447      0.031   14.645    0.000
#>    .x1              0.140      0.018    7.836    0.000
#>    .x2              0.317      0.017   18.641    0.000
#>    .x3              0.228      0.013   17.494    0.000
#>    .z1              0.161      0.016   10.218    0.000
#>    .z2              0.306      0.019   16.521    0.000
#>    .z3              0.241      0.016   15.538    0.000
#>    .y1              0.162      0.016   10.185    0.000
#>    .y2              0.274      0.022   12.375    0.000
#>    .y3              0.250      0.022   11.442    0.000
#>    .w1              0.160      0.020    8.065    0.000
#>    .w2              0.305      0.022   13.799    0.000
#>    .w3              0.229      0.018   12.426    0.000
#>     Y~1             0.086      0.021    4.087    0.000
#>     Y~X             0.018      0.007    2.747    0.006
#>     Y~Z             0.105      0.015    6.773    0.000
#>     W~1             0.059      0.012    4.977    0.000
#>     W~X             0.095      0.014    6.676    0.000
#>     W~Z             0.144      0.027    5.270    0.000

Ordered Indicators

fit_slopes_ord <- pls(
  slopes_model,
  data      = randomSlopesOrdered,
  bootstrap = TRUE,
  boot.R    = 50,
  ordered   = colnames(randomSlopesOrdered) # explicitly specify variables as ordered
)
summary(fit_slopes_ord)
#> plssem (0.1.1) 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
#> 
#> Fit Measures:
#>   Chi-Square                                   263.372
#>   Degrees of Freedom                                49
#>   SRMR                                           0.010
#>   RMSEA                                          0.030
#> 
#> 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.539
#>   W                                              0.550
#> 
#> Latent Variables:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   X =~          
#>     x1              0.933      0.013   74.052    0.000
#>     x2              0.818      0.011   74.996    0.000
#>     x3              0.888      0.011   83.460    0.000
#>   Z =~          
#>     z1              0.917      0.013   72.433    0.000
#>     z2              0.845      0.012   72.353    0.000
#>     z3              0.871      0.013   67.845    0.000
#>   Y =~          
#>     y1              0.918      0.011   80.060    0.000
#>     y2              0.843      0.013   65.426    0.000
#>     y3              0.869      0.014   63.020    0.000
#>   W =~          
#>     w1              0.914      0.015   62.284    0.000
#>     w2              0.825      0.019   43.446    0.000
#>     w3              0.886      0.014   63.524    0.000
#> 
#> Regressions:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   Y ~           
#>     X               0.290      0.020   14.553    0.000
#>     Z               0.451      0.032   13.963    0.000
#>   W ~           
#>     X               0.391      0.042    9.281    0.000
#>     Z               0.245      0.061    4.008    0.000
#> 
#> Intercepts:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>    .Y               0.011      0.005    2.500    0.012
#>    .W               0.009      0.007    1.370    0.171
#> 
#> Covariances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   X ~~          
#>     Z               0.168      0.012   13.759    0.000
#>   Y~X ~~        
#>     Y~1            -0.007      0.004   -1.526    0.127
#>   Y~Z ~~        
#>     Y~1            -0.025      0.013   -1.902    0.057
#>     Y~X             0.012      0.010    1.277    0.202
#>   W~X ~~        
#>     W~1             0.004      0.010    0.392    0.695
#>   W~Z ~~        
#>     W~1             0.010      0.014    0.697    0.486
#>     W~X             0.016      0.009    1.803    0.071
#> 
#> Variances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>     X               1.000                             
#>     Z               1.000                             
#>    .Y               0.461      0.025   18.742    0.000
#>    .W               0.450      0.040   11.276    0.000
#>    .x1              0.130      0.023    5.590    0.000
#>    .x2              0.331      0.018   18.508    0.000
#>    .x3              0.211      0.019   11.156    0.000
#>    .z1              0.159      0.023    6.800    0.000
#>    .z2              0.286      0.020   14.542    0.000
#>    .z3              0.242      0.022   10.839    0.000
#>    .y1              0.156      0.021    7.436    0.000
#>    .y2              0.289      0.022   13.316    0.000
#>    .y3              0.246      0.024   10.263    0.000
#>    .w1              0.165      0.027    6.208    0.000
#>    .w2              0.319      0.031   10.240    0.000
#>    .w3              0.215      0.025    8.708    0.000
#>     Y~1             0.084      0.023    3.693    0.000
#>     Y~X             0.018      0.005    3.806    0.000
#>     Y~Z             0.101      0.015    6.807    0.000
#>     W~1             0.058      0.012    4.732    0.000
#>     W~X             0.099      0.016    6.375    0.000
#>     W~Z             0.142      0.026    5.494    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,
  boot.R    = 50
)
summary(fit_intercepts_cont)
#> plssem (0.1.1) ended normally after 2 iterations
#> 
#>   Estimator                                   PLSc-MLM
#>   Link                                          PROBIT
#>                                                       
#>   Number of observations                         10000
#>   Number of iterations                               2
#>   Number of latent variables                         1
#>   Number of observed variables                       9
#> 
#> Fit Measures:
#>   Chi-Square                                    24.372
#>   Degrees of Freedom                                10
#>   SRMR                                           0.003
#>   RMSEA                                          0.012
#> 
#> R-squared (indicators):
#>   y1                                             0.891
#>   y2                                             0.785
#>   y3                                             0.814
#> 
#> R-squared (latents):
#>   f                                              0.744
#> 
#> Latent Variables:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   f =~          
#>     y1              0.944      0.009   99.515    0.000
#>     y2              0.886      0.009   98.843    0.000
#>     y3              0.902      0.009   95.623    0.000
#> 
#> Regressions:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   f ~           
#>     x1              0.238      0.008   28.898    0.000
#>     x2              0.162      0.007   23.392    0.000
#>     x3              0.077      0.006   13.149    0.000
#>     w1              0.128      0.038    3.378    0.001
#>     w2              0.091      0.033    2.746    0.006
#> 
#> Intercepts:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>    .f              -0.013      0.008   -1.636    0.102
#> 
#> Covariances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   x1 ~~         
#>     x2              0.104      0.011    9.068    0.000
#>     x3              0.004      0.010    0.426    0.670
#>     w1              0.000                             
#>     w2              0.000                             
#>   x2 ~~         
#>     x3              0.097      0.010    9.333    0.000
#>     w1              0.000                             
#>     w2              0.000                             
#>   x3 ~~         
#>     w1              0.000                             
#>     w2              0.000                             
#>   w1 ~~         
#>     w2             -0.041      0.050   -0.819    0.413
#> 
#> Variances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>    .f               0.256      0.014   17.816    0.000
#>     x1              1.000                             
#>     x2              1.000                             
#>     x3              1.000                             
#>     w1              1.000                             
#>     w2              1.000                             
#>    .y1              0.109      0.018    6.056    0.000
#>    .y2              0.215      0.016   13.486    0.000
#>    .y3              0.186      0.017   10.985    0.000
#>     f~1             0.621      0.021   30.153    0.000

Ordered Indicators

fit_intercepts_ord <- pls(
  intercepts_model,
  data      = randomInterceptsOrdered,
  bootstrap = TRUE,
  boot.R    = 50,
  ordered   = colnames(randomInterceptsOrdered) # explicitly specify variables as ordered
)
summary(fit_intercepts_ord)
#> plssem (0.1.1) 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
#> 
#> Fit Measures:
#>   Chi-Square                                    25.508
#>   Degrees of Freedom                                10
#>   SRMR                                           0.003
#>   RMSEA                                          0.012
#> 
#> R-squared (indicators):
#>   y1                                             0.885
#>   y2                                             0.788
#>   y3                                             0.809
#> 
#> R-squared (latents):
#>   f                                              0.723
#> 
#> Latent Variables:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   f =~          
#>     y1              0.941      0.012   78.652    0.000
#>     y2              0.888      0.013   69.164    0.000
#>     y3              0.899      0.014   63.554    0.000
#> 
#> Regressions:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   f ~           
#>     x1              0.241      0.009   27.424    0.000
#>     x2              0.158      0.008   20.447    0.000
#>     x3              0.080      0.007   11.671    0.000
#>     w1              0.121      0.044    2.726    0.006
#>     w2              0.081      0.043    1.897    0.058
#> 
#> Intercepts:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>    .f              -0.012      0.008   -1.422    0.155
#> 
#> Covariances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   x1 ~~         
#>     x2              0.110      0.011   10.357    0.000
#>     x3              0.012      0.013    0.899    0.369
#>     w1              0.001      0.004    0.214    0.831
#>     w2              0.001      0.003    0.361    0.718
#>   x2 ~~         
#>     x3              0.099      0.011    8.741    0.000
#>     w1             -0.003      0.003   -1.012    0.311
#>     w2             -0.001      0.003   -0.455    0.649
#>   x3 ~~         
#>     w1             -0.003      0.003   -1.034    0.301
#>     w2              0.003      0.003    1.071    0.284
#>   w1 ~~         
#>     w2             -0.025      0.054   -0.454    0.650
#> 
#> Variances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>    .f               0.277      0.015   18.173    0.000
#>     x1              1.000                             
#>     x2              1.000                             
#>     x3              1.000                             
#>     w1              1.000                             
#>     w2              1.000                             
#>    .y1              0.115      0.023    5.067    0.000
#>    .y2              0.212      0.023    9.284    0.000
#>    .y3              0.191      0.025    7.543    0.000
#>     f~1             0.601      0.021   28.497    0.000