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.3) ended normally after 67 iterations
#>   Estimator                                 MCPLSc-MLM
#>   Link                                          LINEAR
#>                                                       
#>   Number of observations                          5000
#>   Number of iterations                              67
#>   Number of latent variables                         4
#>   Number of observed variables                      18
#> 
#> Fit Measures:
#>   Chi-Square                                   130.171
#>   Degrees of Freedom                                49
#>   SRMR                                           0.011
#>   RMSEA                                          0.018
#> 
#> R-squared (indicators):
#>   x1                                             0.859
#>   x2                                             0.685
#>   x3                                             0.772
#>   z1                                             0.838
#>   z2                                             0.698
#>   z3                                             0.760
#>   y1                                             0.838
#>   y2                                             0.725
#>   y3                                             0.751
#>   w1                                             0.840
#>   w2                                             0.697
#>   w3                                             0.769
#> 
#> R-squared (latents):
#>   Y                                              0.326
#>   W                                              0.250
#> 
#> Latent Variables:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   X =~          
#>     x1              0.927      0.011   87.527    0.000
#>     x2              0.828      0.011   75.744    0.000
#>     x3              0.879      0.007  124.785    0.000
#>   Z =~          
#>     z1              0.916      0.011   81.649    0.000
#>     z2              0.836      0.008   98.562    0.000
#>     z3              0.872      0.009   99.158    0.000
#>   Y =~          
#>     y1              0.915      0.008  113.219    0.000
#>     y2              0.852      0.013   67.255    0.000
#>     y3              0.867      0.015   59.517    0.000
#>   W =~          
#>     w1              0.917      0.010   87.825    0.000
#>     w2              0.835      0.014   61.636    0.000
#>     w3              0.877      0.015   57.172    0.000
#> 
#> Regressions:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   Y ~           
#>     X               0.292      0.019   15.688    0.000
#>     Z               0.444      0.041   10.823    0.000
#>   W ~           
#>     X               0.394      0.045    8.821    0.000
#>     Z               0.248      0.059    4.176    0.000
#> 
#> Covariances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   X ~~          
#>     Z               0.176      0.013   13.387    0.000
#>   Y~X ~~        
#>     Y~1            -0.003      0.005   -0.611    0.541
#>   Y~Z ~~        
#>     Y~1            -0.025      0.014   -1.793    0.073
#>     Y~X             0.011      0.009    1.205    0.228
#>   W~X ~~        
#>     W~1             0.002      0.010    0.234    0.815
#>   W~Z ~~        
#>     W~1             0.009      0.011    0.824    0.410
#>     W~X             0.010      0.014    0.694    0.487
#> 
#> Variances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>     X               1.000                             
#>     Z               1.000                             
#>    .Y               0.674      0.040   16.967    0.000
#>    .W               0.750      0.055   13.765    0.000
#>    .x1              0.141      0.020    7.197    0.000
#>    .x2              0.315      0.018   17.433    0.000
#>    .x3              0.228      0.012   18.395    0.000
#>    .z1              0.162      0.021    7.882    0.000
#>    .z2              0.302      0.014   21.308    0.000
#>    .z3              0.240      0.015   15.693    0.000
#>    .y1              0.162      0.015   10.935    0.000
#>    .y2              0.275      0.022   12.748    0.000
#>    .y3              0.249      0.025    9.865    0.000
#>    .w1              0.160      0.019    8.341    0.000
#>    .w2              0.303      0.023   13.374    0.000
#>    .w3              0.231      0.027    8.575    0.000
#>     Y~1             0.084      0.021    4.045    0.000
#>     Y~X             0.018      0.006    3.171    0.002
#>     Y~Z             0.102      0.019    5.287    0.000
#>     W~1             0.057      0.010    5.693    0.000
#>     W~X             0.094      0.022    4.350    0.000
#>     W~Z             0.150      0.026    5.865    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.3) ended normally after 65 iterations
#>   Estimator                              MCOrdPLSc-MLM
#>   Link                                          PROBIT
#>                                                       
#>   Number of observations                          5000
#>   Number of iterations                              65
#>   Number of latent variables                         4
#>   Number of observed variables                      18
#> 
#> Fit Measures:
#>   Chi-Square                                   148.880
#>   Degrees of Freedom                                49
#>   SRMR                                           0.012
#>   RMSEA                                          0.020
#> 
#> R-squared (indicators):
#>   x1                                             0.867
#>   x2                                             0.672
#>   x3                                             0.789
#>   z1                                             0.840
#>   z2                                             0.716
#>   z3                                             0.757
#>   y1                                             0.848
#>   y2                                             0.711
#>   y3                                             0.754
#>   w1                                             0.842
#>   w2                                             0.678
#>   w3                                             0.785
#> 
#> R-squared (latents):
#>   Y                                              0.331
#>   W                                              0.234
#> 
#> Latent Variables:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   X =~          
#>     x1              0.931      0.015   62.928    0.000
#>     x2              0.820      0.010   83.920    0.000
#>     x3              0.888      0.012   73.611    0.000
#>   Z =~          
#>     z1              0.917      0.014   66.876    0.000
#>     z2              0.846      0.014   58.468    0.000
#>     z3              0.870      0.015   58.937    0.000
#>   Y =~          
#>     y1              0.921      0.016   56.143    0.000
#>     y2              0.843      0.024   35.529    0.000
#>     y3              0.868      0.016   52.969    0.000
#>   W =~          
#>     w1              0.918      0.016   55.713    0.000
#>     w2              0.824      0.021   38.545    0.000
#>     w3              0.886      0.024   37.354    0.000
#> 
#> Regressions:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   Y ~           
#>     X               0.291      0.028   10.468    0.000
#>     Z               0.450      0.044   10.134    0.000
#>   W ~           
#>     X               0.375      0.049    7.730    0.000
#>     Z               0.248      0.056    4.440    0.000
#> 
#> Covariances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   X ~~          
#>     Z               0.168      0.013   13.193    0.000
#>   Y~X ~~        
#>     Y~1            -0.006      0.007   -0.834    0.404
#>   Y~Z ~~        
#>     Y~1            -0.022      0.016   -1.411    0.158
#>     Y~X             0.013      0.011    1.180    0.238
#>   W~X ~~        
#>     W~1             0.004      0.013    0.328    0.743
#>   W~Z ~~        
#>     W~1             0.012      0.020    0.615    0.538
#>     W~X             0.013      0.013    0.982    0.326
#> 
#> Thresholds:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>     x1|t1          -2.853      0.032  -87.991    0.000
#>     x1|t2          -2.005      0.039  -51.742    0.000
#>     x1|t3          -0.945      0.019  -49.812    0.000
#>     x1|t4           0.140      0.019    7.438    0.000
#>     x1|t5           1.090      0.022   49.254    0.000
#>     x1|t6           2.179      0.042   52.338    0.000
#>     x2|t1          -1.945      0.036  -53.466    0.000
#>     x2|t2          -0.981      0.018  -53.919    0.000
#>     x2|t3           0.288      0.015   18.616    0.000
#>     x2|t4           0.966      0.018   54.492    0.000
#>     x2|t5           2.258      0.050   45.207    0.000
#>     x2|t6           3.095      0.033   95.070    0.000
#>     x3|t1          -2.014      0.040  -50.970    0.000
#>     x3|t2          -1.145      0.022  -52.150    0.000
#>     x3|t3          -0.176      0.019   -9.379    0.000
#>     x3|t4           0.670      0.019   34.500    0.000
#>     x3|t5           1.864      0.036   51.996    0.000
#>     x3|t6           2.708      0.032   85.391    0.000
#>     z1|t1          -2.130      0.045  -47.777    0.000
#>     z1|t2          -1.330      0.023  -58.791    0.000
#>     z1|t3          -0.213      0.017  -12.364    0.000
#>     z1|t4           0.841      0.023   36.979    0.000
#>     z1|t5           1.890      0.046   41.231    0.000
#>     z1|t6           2.607      0.035   74.689    0.000
#>     z2|t1          -2.476      0.034  -73.870    0.000
#>     z2|t2          -1.442      0.024  -58.951    0.000
#>     z2|t3          -0.284      0.022  -13.008    0.000
#>     z2|t4           0.813      0.016   52.368    0.000
#>     z2|t5           1.528      0.026   58.037    0.000
#>     z2|t6           2.569      0.042   61.004    0.000
#>     z3|t1          -2.659      0.046  -57.352    0.000
#>     z3|t2          -2.083      0.034  -61.064    0.000
#>     z3|t3          -1.161      0.024  -48.605    0.000
#>     z3|t4           0.095      0.020    4.645    0.000
#>     z3|t5           0.977      0.026   37.979    0.000
#>     z3|t6           1.964      0.046   42.312    0.000
#>     y1|t1          -2.670      0.097  -27.629    0.000
#>     y1|t2          -1.769      0.085  -20.773    0.000
#>     y1|t3          -0.652      0.050  -13.020    0.000
#>     y1|t4           0.395      0.045    8.723    0.000
#>     y1|t5           1.434      0.075   19.176    0.000
#>     y1|t6           2.397      0.080   29.996    0.000
#>     y2|t1          -2.549      0.052  -49.315    0.000
#>     y2|t2          -1.704      0.081  -20.925    0.000
#>     y2|t3          -0.555      0.047  -11.832    0.000
#>     y2|t4           0.222      0.034    6.454    0.000
#>     y2|t5           1.400      0.050   28.144    0.000
#>     y2|t6           2.360      0.075   31.630    0.000
#>     y3|t1          -1.906      0.100  -18.965    0.000
#>     y3|t2          -0.938      0.047  -20.131    0.000
#>     y3|t3          -0.091      0.048   -1.912    0.056
#>     y3|t4           1.104      0.045   24.685    0.000
#>     y3|t5           1.985      0.065   30.470    0.000
#>     y3|t6           3.094      0.140   22.077    0.000
#>     w1|t1          -2.787      0.090  -31.140    0.000
#>     w1|t2          -1.653      0.059  -27.812    0.000
#>     w1|t3          -0.773      0.047  -16.624    0.000
#>     w1|t4           0.259      0.035    7.350    0.000
#>     w1|t5           1.073      0.056   19.082    0.000
#>     w1|t6           2.373      0.089   26.738    0.000
#>     w2|t1          -2.486      0.037  -67.589    0.000
#>     w2|t2          -1.645      0.063  -26.270    0.000
#>     w2|t3          -0.596      0.035  -17.091    0.000
#>     w2|t4           0.561      0.029   19.133    0.000
#>     w2|t5           1.408      0.046   30.338    0.000
#>     w2|t6           2.455      0.072   34.246    0.000
#>     w3|t1          -2.376      0.094  -25.401    0.000
#>     w3|t2          -1.411      0.071  -19.840    0.000
#>     w3|t3          -0.715      0.040  -17.726    0.000
#>     w3|t4           0.566      0.038   15.026    0.000
#>     w3|t5           1.528      0.070   21.944    0.000
#>     w3|t6           2.265      0.075   30.185    0.000
#> 
#> Variances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>     X               1.000                             
#>     Z               1.000                             
#>    .Y               0.669      0.050   13.379    0.000
#>    .W               0.766      0.050   15.271    0.000
#>    .x1              0.133      0.028    4.819    0.000
#>    .x2              0.328      0.016   20.455    0.000
#>    .x3              0.211      0.021    9.836    0.000
#>    .z1              0.160      0.025    6.349    0.000
#>    .z2              0.284      0.024   11.616    0.000
#>    .z3              0.243      0.026    9.459    0.000
#>    .y1              0.152      0.030    5.049    0.000
#>    .y2              0.289      0.040    7.227    0.000
#>    .y3              0.246      0.028    8.633    0.000
#>    .w1              0.158      0.030    5.227    0.000
#>    .w2              0.322      0.035    9.140    0.000
#>    .w3              0.215      0.042    5.103    0.000
#>     Y~1             0.083      0.027    3.095    0.002
#>     Y~X             0.018      0.006    3.247    0.001
#>     Y~Z             0.102      0.027    3.797    0.000
#>     W~1             0.060      0.017    3.594    0.000
#>     W~X             0.102      0.017    6.096    0.000
#>     W~Z             0.154      0.040    3.896    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.3) ended normally after 58 iterations
#>   Estimator                                  MCPLSc-MLM
#>   Link                                           LINEAR
#>                                                        
#>   Number of observations                          10000
#>   Number of iterations                               58
#>   Number of latent variables                          1
#>   Number of observed variables                        9
#> 
#> Fit Measures:
#>   Chi-Square                                     23.841
#>   Degrees of Freedom                                 10
#>   SRMR                                            0.003
#>   RMSEA                                           0.012
#> 
#> R-squared (indicators):
#>   y1                                              0.890
#>   y2                                              0.787
#>   y3                                              0.814
#> 
#> R-squared (latents):
#>   f                                               0.122
#> 
#> Latent Variables:
#>                  Estimate  Std.Error   z.value  P(>|z|)
#>   f =~          
#>     y1              0.943      0.008   115.270    0.000
#>     y2              0.887      0.010    87.170    0.000
#>     y3              0.902      0.010    93.935    0.000
#> 
#> Regressions:
#>                  Estimate  Std.Error   z.value  P(>|z|)
#>   f ~           
#>     x1              0.236      0.005    49.927    0.000
#>     x2              0.162      0.006    29.183    0.000
#>     x3              0.077      0.007    11.349    0.000
#>     w1              0.127      0.038     3.317    0.001
#>     w2              0.091      0.038     2.425    0.015
#> 
#> Covariances:
#>                  Estimate  Std.Error   z.value  P(>|z|)
#>   x1 ~~         
#>     x2              0.103      0.010    10.495    0.000
#>     x3              0.005      0.011     0.470    0.639
#>     w1              0.001      0.000  1674.785    0.000
#>     w2              0.001      0.000     9.337    0.000
#>   x2 ~~         
#>     x3              0.095      0.011     8.629    0.000
#>     w1             -0.002      0.000  -104.395    0.000
#>     w2             -0.002      0.000   -12.718    0.000
#>   x3 ~~         
#>     w1              0.000      0.000    15.925    0.000
#>     w2              0.000      0.000    -0.586    0.558
#>   w1 ~~         
#>     w2             -0.040      0.050    -0.800    0.424
#> 
#> Variances:
#>                  Estimate  Std.Error   z.value  P(>|z|)
#>    .f               0.878      0.012    72.391    0.000
#>     x1              1.000                              
#>     x2              1.000                              
#>     x3              1.000                              
#>     w1              1.000                              
#>     w2              1.000                              
#>    .y1              0.110      0.015     7.153    0.000
#>    .y2              0.213      0.018    11.789    0.000
#>    .y3              0.186      0.017    10.698    0.000
#>     f~1             0.615      0.037    16.817    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.3) ended normally after 76 iterations
#>   Estimator                               MCOrdPLSc-MLM
#>   Link                                           PROBIT
#>                                                        
#>   Number of observations                          10000
#>   Number of iterations                               76
#>   Number of latent variables                          1
#>   Number of observed variables                        9
#> 
#> Fit Measures:
#>   Chi-Square                                      9.496
#>   Degrees of Freedom                                 10
#>   SRMR                                            0.003
#>   RMSEA                                           0.000
#> 
#> R-squared (indicators):
#>   y1                                              0.881
#>   y2                                              0.786
#>   y3                                              0.816
#> 
#> R-squared (latents):
#>   f                                               0.121
#> 
#> Latent Variables:
#>                  Estimate  Std.Error   z.value  P(>|z|)
#>   f =~          
#>     y1              0.939      0.013    71.575    0.000
#>     y2              0.887      0.016    55.955    0.000
#>     y3              0.903      0.017    52.208    0.000
#> 
#> Regressions:
#>                  Estimate  Std.Error   z.value  P(>|z|)
#>   f ~           
#>     x1              0.242      0.008    30.435    0.000
#>     x2              0.158      0.008    19.805    0.000
#>     x3              0.080      0.008    10.085    0.000
#>     w1              0.123      0.035     3.527    0.000
#>     w2              0.078      0.044     1.767    0.077
#> 
#> Covariances:
#>                  Estimate  Std.Error   z.value  P(>|z|)
#>   x1 ~~         
#>     x2              0.108      0.011    10.188    0.000
#>     x3              0.013      0.011     1.125    0.261
#>     w1              0.000      0.006     0.030    0.976
#>     w2              0.001      0.006     0.096    0.924
#>   x2 ~~         
#>     x3              0.099      0.012     8.341    0.000
#>     w1             -0.004      0.005    -0.848    0.397
#>     w2             -0.001      0.005    -0.154    0.877
#>   x3 ~~         
#>     w1             -0.002      0.004    -0.546    0.585
#>     w2              0.005      0.004     1.190    0.234
#>   w1 ~~         
#>     w2             -0.025      0.056    -0.443    0.657
#> 
#> Thresholds:
#>                  Estimate  Std.Error   z.value  P(>|z|)
#>     y1|t1          -2.531      0.078   -32.527    0.000
#>     y1|t2          -1.743      0.070   -24.735    0.000
#>     y1|t3          -0.427      0.034   -12.565    0.000
#>     y1|t4           0.316      0.038     8.424    0.000
#>     y1|t5           1.343      0.038    35.266    0.000
#>     y1|t6           2.410      0.056    43.211    0.000
#>     y2|t1          -2.811      0.086   -32.839    0.000
#>     y2|t2          -1.873      0.067   -27.941    0.000
#>     y2|t3          -0.984      0.041   -24.181    0.000
#>     y2|t4           0.171      0.034     5.092    0.000
#>     y2|t5           1.034      0.040    25.635    0.000
#>     y2|t6           2.278      0.105    21.693    0.000
#>     y3|t1          -2.204      0.115   -19.098    0.000
#>     y3|t2          -1.245      0.053   -23.617    0.000
#>     y3|t3          -0.014      0.038    -0.374    0.709
#>     y3|t4           0.682      0.031    21.686    0.000
#>     y3|t5           1.816      0.054    33.591    0.000
#>     y3|t6           2.917      0.060    48.565    0.000
#>     x1|t1          -2.394      0.030   -80.441    0.000
#>     x1|t2          -1.347      0.015   -91.330    0.000
#>     x1|t3          -0.664      0.010   -67.424    0.000
#>     x1|t4           0.503      0.008    64.502    0.000
#>     x1|t5           1.514      0.018    81.960    0.000
#>     x1|t6           2.609      0.028    93.574    0.000
#>     x2|t1          -2.832      0.027  -106.084    0.000
#>     x2|t2          -2.025      0.028   -73.604    0.000
#>     x2|t3          -1.012      0.012   -84.876    0.000
#>     x2|t4          -0.115      0.007   -15.790    0.000
#>     x2|t5           0.810      0.011    72.103    0.000
#>     x2|t6           1.906      0.021    90.760    0.000
#>     x3|t1          -2.932      0.027  -106.999    0.000
#>     x3|t2          -1.955      0.021   -91.852    0.000
#>     x3|t3          -0.847      0.009   -90.678    0.000
#>     x3|t4          -0.149      0.007   -20.522    0.000
#>     x3|t5           1.008      0.012    86.268    0.000
#>     x3|t6           2.225      0.050    44.800    0.000
#>     w1|t1          -2.841      0.155   -18.391    0.000
#>     w1|t2          -2.008      0.136   -14.748    0.000
#>     w1|t3          -0.959      0.072   -13.388    0.000
#>     w1|t4           0.109      0.067     1.622    0.105
#>     w1|t5           1.249      0.092    13.637    0.000
#>     w1|t6           2.602      0.158    16.458    0.000
#>     w2|t1          -1.992      0.121   -16.433    0.000
#>     w2|t2          -0.935      0.079   -11.794    0.000
#>     w2|t3          -0.119      0.084    -1.426    0.154
#>     w2|t4           0.734      0.070    10.490    0.000
#>     w2|t5           1.810      0.143    12.625    0.000
#>     w2|t6           2.531      0.186    13.603    0.000
#> 
#> Variances:
#>                  Estimate  Std.Error   z.value  P(>|z|)
#>    .f               0.879      0.011    81.668    0.000
#>     x1              1.000                              
#>     x2              1.000                              
#>     x3              1.000                              
#>     w1              1.000                              
#>     w2              1.000                              
#>    .y1              0.119      0.025     4.819    0.000
#>    .y2              0.214      0.028     7.602    0.000
#>    .y3              0.184      0.031     5.895    0.000
#>     f~1             0.614      0.037    16.430    0.000