Chapter 22 Lavaan Lab 19: Multilevel SEM

In this lab, we will:

build a multilevel CFA model
add covariates at both the between and the within level

Load up the lavaan library:

library(lavaan)

Let’s read in a Mplus example dataset from an online location

Data <- read.table("http://statmodel.com/usersguide/chap9/ex9.6.dat")
names(Data) <- c("y1", "y2", "y3", "y4", "x1", "x2", "w", "clus")

Take a look at the matrix:

head(Data)

##          y1        y2        y3        y4        x1        x2         w clus
## 1  2.203250  1.858861  1.738477  2.244863  1.142800 -0.796987 -0.149501    1
## 2  1.934917  2.127876  0.083120  2.509436  1.949033 -0.122764 -0.149501    1
## 3  0.321955  0.977231 -0.835405  0.558367 -0.716481 -0.767064 -0.149501    1
## 4  0.073154 -1.743092 -2.310271 -1.514332 -2.649131  0.637570 -0.149501    1
## 5 -1.214906  0.452618  0.372610 -1.790372 -0.262916  0.302564 -0.149501    1
## 6  0.298330 -1.820272  0.561335 -2.090582 -0.944963  1.363045  0.319335    2

dim(Data)

## [1] 1000    8

length(unique(Data$clus))

## [1] 110

there are 1000 individual observations in 110 clusters
cluster sizes: 5, 10, 15
4 measures at the within level y1, y2, y3, y4
2 covariates at the within level: x1, x2
1 covariate at the between level: w

22.1 PART I: Multilevel CFA 1: within-only construct

This model specifies the latent variable only at the within level:

model1 <- '
  level: 1
    fw =~ y1 + y2 + y3 + y4

  level: 2
    y1 ~~ y1 + y2 + y3 + y4
    y2 ~~ y2 + y2 + y3
    y3 ~~ y3 + y4
    y4 ~~ y4

    # all variances and covariances are freely estimated
'

Fit the model:

model1fit <- lavaan::sem(model1, data = Data, 
                         cluster = "clus",
                         estimator = 'MLR',
                         fixed.x = FALSE)

lavaan::summary(model1fit, fit.measures = T, standardized = T)

## lavaan 0.6-12 ended normally after 57 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        21
## 
##   Number of observations                          1000
##   Number of clusters [clus]                        110
## 
## Model Test User Model:
##                                               Standard      Robust
##   Test Statistic                                52.360      64.738
##   Degrees of freedom                                 3           3
##   P-value (Chi-square)                           0.000       0.000
##   Scaling correction factor                                  0.809
##     Yuan-Bentler correction (Mplus variant)                       
## 
## Model Test Baseline Model:
## 
##   Test statistic                              2619.397    2815.244
##   Degrees of freedom                                12          12
##   P-value                                        0.000       0.000
##   Scaling correction factor                                  0.930
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.981       0.978
##   Tucker-Lewis Index (TLI)                       0.924       0.912
##                                                                   
##   Robust Comparative Fit Index (CFI)                         0.981
##   Robust Tucker-Lewis Index (TLI)                            0.923
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -6900.518   -6900.518
##   Scaling correction factor                                  1.050
##       for the MLR correction                                      
##   Loglikelihood unrestricted model (H1)      -6874.338   -6874.338
##   Scaling correction factor                                  1.020
##       for the MLR correction                                      
##                                                                   
##   Akaike (AIC)                               13843.036   13843.036
##   Bayesian (BIC)                             13946.099   13946.099
##   Sample-size adjusted Bayesian (BIC)        13879.402   13879.402
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.128       0.143
##   90 Percent confidence interval - lower         0.099       0.111
##   90 Percent confidence interval - upper         0.160       0.178
##   P-value RMSEA <= 0.05                          0.000       0.000
##                                                                   
##   Robust RMSEA                                               0.129
##   90 Percent confidence interval - lower                     0.103
##   90 Percent confidence interval - upper                     0.157
## 
## Standardized Root Mean Square Residual (corr metric):
## 
##   SRMR (within covariance matrix)                0.022       0.022
##   SRMR (between covariance matrix)               0.334       0.334
## 
## Parameter Estimates:
## 
##   Standard errors                             Sandwich
##   Information bread                           Observed
##   Observed information based on                Hessian
## 
## 
## Level 1 [within]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fw =~                                                                 
##     y1                1.000                               1.568    0.843
##     y2                1.019    0.035   28.881    0.000    1.598    0.857
##     y3                0.995    0.041   24.267    0.000    1.560    0.832
##     y4                1.034    0.035   29.613    0.000    1.622    0.851
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.000                               0.000    0.000
##    .y2                0.000                               0.000    0.000
##    .y3                0.000                               0.000    0.000
##    .y4                0.000                               0.000    0.000
##     fw                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.997    0.070   14.220    0.000    0.997    0.289
##    .y2                0.919    0.061   14.990    0.000    0.919    0.265
##    .y3                1.082    0.057   18.977    0.000    1.082    0.308
##    .y4                0.997    0.069   14.422    0.000    0.997    0.275
##     fw                2.458    0.186   13.209    0.000    1.000    1.000
## 
## 
## Level 2 [clus]:
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .y1 ~~                                                                 
##    .y2                0.069    0.027    2.558    0.011    0.069    0.935
##    .y3                0.080    0.027    2.916    0.004    0.080    1.141
##    .y4                0.052    0.024    2.157    0.031    0.052    0.764
##  .y2 ~~                                                                 
##    .y3                0.055    0.023    2.399    0.016    0.055    0.935
##  .y3 ~~                                                                 
##    .y4                0.040    0.026    1.537    0.124    0.040    0.727
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1               -0.063    0.089   -0.716    0.474   -0.063   -0.215
##    .y2               -0.060    0.088   -0.674    0.500   -0.060   -0.239
##    .y3               -0.026    0.086   -0.301    0.763   -0.026   -0.109
##    .y4               -0.005    0.089   -0.053    0.958   -0.005   -0.020
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.088    0.033    2.679    0.007    0.088    1.000
##    .y2                0.062    0.037    1.687    0.092    0.062    1.000
##    .y3                0.056    0.035    1.608    0.108    0.056    1.000
##    .y4                0.054    0.038    1.423    0.155    0.054    1.000

22.2 PART II: Multilevel CFA 2: Between-only construct

Example: construct reflects self-reported ‘school climate’ measured by a questionnaire filled in by the school principles
We will only have one response for each school
We may collect other variables from students/teachers in the schools though

Note that the following model syntax:

model2.wrong <- '
  level: 1
  # perhaps other level-1 variables

  level: 2
    fb =~ y1 + y2 + y3 + y4
'

model2fit <- lavaan::sem(model2.wrong, data = Data, 
                         cluster = "clus",
                         estimator = 'MLR',
                         fixed.x = FALSE)

Error in lav_partable_vnames(tmp.lav, type = "ov", level = tmp.level.values[l]) : lavaan ERROR: level column does not contain value `1'

won’t work because there is nothing at level 1. Instead, specify this model just like a regular CFA model:

model2 <- '
    fb =~ y1 + y2 + y3 + y4
'

Fit the model:

model2fit <- lavaan::sem(model2, data = Data, 
                         cluster = "clus",
                         estimator = 'MLR',
                         fixed.x = FALSE)

lavaan::summary(model2fit, fit.measures = T, standardized = T)

## lavaan 0.6-12 ended normally after 20 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        12
## 
##   Number of observations                          1000
##   Number of clusters [clus]                        110
## 
## Model Test User Model:
##                                               Standard      Robust
##   Test Statistic                                 0.108       0.107
##   Degrees of freedom                                 2           2
##   P-value (Chi-square)                           0.947       0.948
##   Scaling correction factor                                  1.013
##     Yuan-Bentler correction (Mplus variant)                       
## 
## Model Test Baseline Model:
## 
##   Test statistic                              2729.902    5889.388
##   Degrees of freedom                                 6           6
##   P-value                                        0.000       0.000
##   Scaling correction factor                                  0.464
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    1.000       1.000
##   Tucker-Lewis Index (TLI)                       1.002       1.001
##                                                                   
##   Robust Comparative Fit Index (CFI)                         1.000
##   Robust Tucker-Lewis Index (TLI)                            1.002
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -6907.183   -6907.183
##   Scaling correction factor                                  1.205
##       for the MLR correction                                      
##   Loglikelihood unrestricted model (H1)      -6907.129   -6907.129
##   Scaling correction factor                                  1.178
##       for the MLR correction                                      
##                                                                   
##   Akaike (AIC)                               13838.367   13838.367
##   Bayesian (BIC)                             13897.260   13897.260
##   Sample-size adjusted Bayesian (BIC)        13859.147   13859.147
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.000       0.000
##   90 Percent confidence interval - lower         0.000       0.000
##   90 Percent confidence interval - upper         0.007       0.003
##   P-value RMSEA <= 0.05                          0.995       0.996
##                                                                   
##   Robust RMSEA                                               0.000
##   90 Percent confidence interval - lower                     0.000
##   90 Percent confidence interval - upper                     0.006
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.001       0.001
## 
## Parameter Estimates:
## 
##   Standard errors                        Robust.cluster
##   Information                                  Observed
##   Observed information based on                 Hessian
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fb =~                                                                 
##     y1                1.000                               1.638    0.856
##     y2                0.995    0.031   31.694    0.000    1.629    0.860
##     y3                0.984    0.037   26.944    0.000    1.612    0.842
##     y4                1.003    0.029   34.582    0.000    1.643    0.851
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1               -0.064    0.090   -0.707    0.479   -0.064   -0.033
##    .y2               -0.058    0.089   -0.652    0.515   -0.058   -0.031
##    .y3               -0.026    0.087   -0.303    0.762   -0.026   -0.014
##    .y4               -0.011    0.091   -0.119    0.905   -0.011   -0.006
##     fb                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.983    0.065   15.151    0.000    0.983    0.268
##    .y2                0.937    0.060   15.657    0.000    0.937    0.261
##    .y3                1.065    0.058   18.477    0.000    1.065    0.291
##    .y4                1.030    0.064   16.008    0.000    1.030    0.276
##     fb                2.684    0.178   15.101    0.000    1.000    1.000

22.3 PART III: Multilevel CFA 3: Shared cross-level construct

This model specifies the latent variable both at the within and the between level;
However, the latent variable only makes sense at the between level so SEM model is only built at the between level;
The indicators are correlated at the within level;

model3 <- '
  level: 1
    y1 ~~ y1 + y2 + y3 + y4
    y2 ~~ y2 + y3 + y4
    y3 ~~ y3 + y4
    y4 ~~ y4

  level: 2
    fs =~ y1 + y2 + y3 + y4

 # Fix Significant Heywood Cases
  y2 ~~ v2*y2
  y4 ~~ v4*y4
  v2 > 0
  v4 > 0 
'

Fit the model:

model3fit <- lavaan::sem(model3, data = Data, 
                         cluster = "clus",
                         estimator = 'MLR',
                         fixed.x = FALSE)

## Warning in lav_object_post_check(object): lavaan WARNING: some estimated ov variances are negative

lavaan::summary(model3fit, standardized = T, fit.measures = T)

## lavaan 0.6-11 ended normally after 94 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        21
##   Number of inequality constraints                   2
##                                                       
##   Number of observations                          1000
##   Number of clusters [clus]                        110
##                                                       
## Model Test User Model:
##                                                Standard      Robust
##   Test Statistic                                630.445     342.869
##   Degrees of freedom                                  3           3
##   P-value (Chi-square)                            0.000       0.000
##   Scaling correction factor                                   1.839
##        Yuan-Bentler correction (Mplus variant)                     
## 
## Model Test Baseline Model:
## 
##   Test statistic                              2619.397    2815.244
##   Degrees of freedom                                12          12
##   P-value                                        0.000       0.000
##   Scaling correction factor                                  0.930
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.759       0.879
##   Tucker-Lewis Index (TLI)                       0.037       0.515
##                                                                   
##   Robust Comparative Fit Index (CFI)                         0.760
##   Robust Tucker-Lewis Index (TLI)                            0.042
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -7189.561   -7189.561
##   Scaling correction factor                                  0.903
##       for the MLR correction                                      
##   Loglikelihood unrestricted model (H1)      -6874.338   -6874.338
##   Scaling correction factor                                  1.020
##       for the MLR correction                                      
##                                                                   
##   Akaike (AIC)                               14421.121   14421.121
##   Bayesian (BIC)                             14524.184   14524.184
##   Sample-size adjusted Bayesian (BIC)        14457.487   14457.487
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.457       0.337
##   90 Percent confidence interval - lower         0.428       0.315
##   90 Percent confidence interval - upper         0.488       0.359
##   P-value RMSEA <= 0.05                          0.000       0.000
##                                                                   
##   Robust RMSEA                                               0.456
##   90 Percent confidence interval - lower                     0.416
##   90 Percent confidence interval - upper                     0.498
## 
## Standardized Root Mean Square Residual (corr metric):
## 
##   SRMR (within covariance matrix)                0.266       0.266
##   SRMR (between covariance matrix)               0.010       0.010
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Observed
##   Observed information based on                Hessian
## 
## 
## Level 1 [within]:
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .y1 ~~                                                                 
##    .y2                1.289    0.106   12.112    0.000    1.289    0.480
##    .y3                1.417    0.076   18.696    0.000    1.417    0.580
##    .y4                1.209    0.106   11.409    0.000    1.209    0.440
##  .y2 ~~                                                                 
##    .y3                1.245    0.107   11.614    0.000    1.245    0.458
##  .y3 ~~                                                                 
##    .y4                1.228    0.108   11.332    0.000    1.228    0.442
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.000                               0.000    0.000
##    .y2                0.000                               0.000    0.000
##    .y3                0.000                               0.000    0.000
##    .y4                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                2.417    0.095   25.385    0.000    2.417    1.000
##    .y2                2.983    0.137   21.836    0.000    2.983    1.000
##    .y3                2.474    0.098   25.227    0.000    2.474    1.000
##    .y4                3.121    0.143   21.864    0.000    3.121    1.000
## 
## 
## Level 2 [clus]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fs =~                                                                 
##     y1                1.000                               0.772    1.001
##     y2                1.044    0.079   13.219    0.000    0.807    1.000
##     y3                0.963    0.066   14.530    0.000    0.744    1.008
##     y4                1.044    0.083   12.581    0.000    0.807    1.000
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1               -0.076    0.090   -0.853    0.394   -0.076   -0.099
##    .y2               -0.071    0.095   -0.747    0.455   -0.071   -0.088
##    .y3               -0.039    0.087   -0.445    0.656   -0.039   -0.053
##    .y4               -0.024    0.096   -0.248    0.804   -0.024   -0.030
##     fs                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y2        (v2)   -0.000       NA                     -0.000   -0.000
##    .y4        (v4)   -0.000                              -0.000   -0.000
##    .y1               -0.002    0.020   -0.077    0.939   -0.002   -0.003
##    .y3               -0.009    0.019   -0.455    0.649   -0.009   -0.016
##     fs                0.597    0.120    4.960    0.000    1.000    1.000
## 
## Constraints:
##                                                |Slack|
##     v2 - 0                                       0.000
##     v4 - 0                                       0.000

22.4 PART IV: Multilevel CFA 4: Configural construct

Model 4a specifies the latent variable both at the within and the between level;
The CFA at each level should have the same factor structure, but not necessarily the same parameter estimates;

model4a <- '
  level: 1
    fw =~ y1 + y2 + y3 + y4

  level: 2
    fb =~ y1 + y2 + y3 + y4
'

model4afit <- lavaan::sem(model4a, data = Data, 
                         cluster = "clus",
                         estimator = 'MLR',
                         fixed.x = FALSE)

## Warning in lav_object_post_check(object): lavaan WARNING: some estimated ov variances are negative

lavaan::summary(model4afit, fit.measures = T, standardized = T)

## lavaan 0.6-12 ended normally after 50 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        20
## 
##   Number of observations                          1000
##   Number of clusters [clus]                        110
## 
## Model Test User Model:
##                                               Standard      Robust
##   Test Statistic                                 0.000       0.000
##   Degrees of freedom                                 4           4
##   P-value (Chi-square)                           1.000       1.000
##   Scaling correction factor                                  0.977
##     Yuan-Bentler correction (Mplus variant)                       
## 
## Model Test Baseline Model:
## 
##   Test statistic                              2619.397    2815.244
##   Degrees of freedom                                12          12
##   P-value                                        0.000       0.000
##   Scaling correction factor                                  0.930
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    1.000       1.000
##   Tucker-Lewis Index (TLI)                       1.005       1.004
##                                                                   
##   Robust Comparative Fit Index (CFI)                         1.000
##   Robust Tucker-Lewis Index (TLI)                            1.004
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -6874.157   -6874.157
##   Scaling correction factor                                  1.028
##       for the MLR correction                                      
##   Loglikelihood unrestricted model (H1)      -6874.338   -6874.338
##   Scaling correction factor                                  1.020
##       for the MLR correction                                      
##                                                                   
##   Akaike (AIC)                               13788.314   13788.314
##   Bayesian (BIC)                             13886.469   13886.469
##   Sample-size adjusted Bayesian (BIC)        13822.948   13822.948
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.000       0.000
##   90 Percent confidence interval - lower         0.000       0.000
##   90 Percent confidence interval - upper         0.000       0.000
##   P-value RMSEA <= 0.05                          1.000       1.000
##                                                                   
##   Robust RMSEA                                               0.000
##   90 Percent confidence interval - lower                     0.000
##   90 Percent confidence interval - upper                     0.000
## 
## Standardized Root Mean Square Residual (corr metric):
## 
##   SRMR (within covariance matrix)                0.001       0.001
##   SRMR (between covariance matrix)               0.003       0.003
## 
## Parameter Estimates:
## 
##   Standard errors                             Sandwich
##   Information bread                           Observed
##   Observed information based on                Hessian
## 
## 
## Level 1 [within]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fw =~                                                                 
##     y1                1.000                               1.486    0.831
##     y2                0.992    0.035   28.208    0.000    1.474    0.836
##     y3                0.992    0.041   24.110    0.000    1.474    0.818
##     y4                1.007    0.035   28.619    0.000    1.496    0.830
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.000                               0.000    0.000
##    .y2                0.000                               0.000    0.000
##    .y3                0.000                               0.000    0.000
##    .y4                0.000                               0.000    0.000
##     fw                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.989    0.071   13.980    0.000    0.989    0.309
##    .y2                0.934    0.062   15.187    0.000    0.934    0.301
##    .y3                1.076    0.057   18.812    0.000    1.076    0.331
##    .y4                1.014    0.070   14.478    0.000    1.014    0.312
##     fw                2.207    0.165   13.371    0.000    1.000    1.000
## 
## 
## Level 2 [clus]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fb =~                                                                 
##     y1                1.000                               0.700    1.006
##     y2                1.005    0.086   11.732    0.000    0.704    0.998
##     y3                0.947    0.094   10.117    0.000    0.663    1.013
##     y4                0.984    0.080   12.240    0.000    0.689    0.983
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1               -0.072    0.089   -0.815    0.415   -0.072   -0.104
##    .y2               -0.067    0.088   -0.755    0.450   -0.067   -0.095
##    .y3               -0.035    0.086   -0.404    0.687   -0.035   -0.053
##    .y4               -0.018    0.089   -0.205    0.838   -0.018   -0.026
##     fb                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1               -0.006    0.018   -0.339    0.734   -0.006   -0.012
##    .y2                0.002    0.020    0.091    0.928    0.002    0.004
##    .y3               -0.011    0.020   -0.541    0.589   -0.011   -0.025
##    .y4                0.016    0.021    0.786    0.432    0.016    0.033
##     fb                0.490    0.120    4.099    0.000    1.000    1.000

Model 4b specifies the same CFA at each level and requires the same factor loadings;

model4b <- '
  level: 1
    fw =~ a*y1 + b*y2 + c*y3 + d*y4

  level: 2
    fb =~ a*y1 + b*y2 + c*y3 + d*y4
'

model4bfit <- lavaan::sem(model4b, data = Data, 
                         cluster = "clus",
                         estimator = 'MLR',
                         fixed.x = FALSE)

## Warning in lav_object_post_check(object): lavaan WARNING: some estimated ov variances are negative

lavaan::summary(model4bfit, fit.measures = T, standardized = T)

## lavaan 0.6-12 ended normally after 44 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        20
##   Number of equality constraints                     3
## 
##   Number of observations                          1000
##   Number of clusters [clus]                        110
## 
## Model Test User Model:
##                                               Standard      Robust
##   Test Statistic                                 0.057       0.059
##   Degrees of freedom                                 7           7
##   P-value (Chi-square)                           1.000       1.000
##   Scaling correction factor                                  0.967
##     Yuan-Bentler correction (Mplus variant)                       
## 
## Model Test Baseline Model:
## 
##   Test statistic                              2619.397    2815.244
##   Degrees of freedom                                12          12
##   P-value                                        0.000       0.000
##   Scaling correction factor                                  0.930
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    1.000       1.000
##   Tucker-Lewis Index (TLI)                       1.005       1.004
##                                                                   
##   Robust Comparative Fit Index (CFI)                         1.000
##   Robust Tucker-Lewis Index (TLI)                            1.004
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -6874.367   -6874.367
##   Scaling correction factor                                  0.885
##       for the MLR correction                                      
##   Loglikelihood unrestricted model (H1)      -6874.338   -6874.338
##   Scaling correction factor                                  1.020
##       for the MLR correction                                      
##                                                                   
##   Akaike (AIC)                               13782.733   13782.733
##   Bayesian (BIC)                             13866.165   13866.165
##   Sample-size adjusted Bayesian (BIC)        13812.172   13812.172
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.000       0.000
##   90 Percent confidence interval - lower         0.000       0.000
##   90 Percent confidence interval - upper         0.000       0.000
##   P-value RMSEA <= 0.05                          1.000       1.000
##                                                                   
##   Robust RMSEA                                               0.000
##   90 Percent confidence interval - lower                     0.000
##   90 Percent confidence interval - upper                     0.000
## 
## Standardized Root Mean Square Residual (corr metric):
## 
##   SRMR (within covariance matrix)                0.002       0.002
##   SRMR (between covariance matrix)               0.003       0.003
## 
## Parameter Estimates:
## 
##   Standard errors                             Sandwich
##   Information bread                           Observed
##   Observed information based on                Hessian
## 
## 
## Level 1 [within]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fw =~                                                                 
##     y1         (a)    1.000                               1.489    0.832
##     y2         (b)    0.995    0.031   31.649    0.000    1.482    0.838
##     y3         (c)    0.983    0.036   27.090    0.000    1.463    0.815
##     y4         (d)    1.003    0.029   34.668    0.000    1.494    0.829
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.000                               0.000    0.000
##    .y2                0.000                               0.000    0.000
##    .y3                0.000                               0.000    0.000
##    .y4                0.000                               0.000    0.000
##     fw                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.987    0.069   14.228    0.000    0.987    0.308
##    .y2                0.930    0.061   15.173    0.000    0.930    0.298
##    .y3                1.081    0.057   18.905    0.000    1.081    0.335
##    .y4                1.015    0.068   14.877    0.000    1.015    0.313
##     fw                2.218    0.158   14.021    0.000    1.000    1.000
## 
## 
## Level 2 [clus]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fb =~                                                                 
##     y1         (a)    1.000                               0.693    1.005
##     y2         (b)    0.995    0.031   31.649    0.000    0.690    0.997
##     y3         (c)    0.983    0.036   27.090    0.000    0.681    1.013
##     y4         (d)    1.003    0.029   34.668    0.000    0.695    0.984
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1               -0.071    0.089   -0.804    0.421   -0.071   -0.103
##    .y2               -0.066    0.088   -0.746    0.456   -0.066   -0.095
##    .y3               -0.034    0.086   -0.396    0.692   -0.034   -0.051
##    .y4               -0.017    0.089   -0.196    0.845   -0.017   -0.025
##     fb                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1               -0.005    0.016   -0.302    0.763   -0.005   -0.010
##    .y2                0.003    0.019    0.170    0.865    0.003    0.007
##    .y3               -0.012    0.020   -0.601    0.548   -0.012   -0.027
##    .y4                0.016    0.020    0.767    0.443    0.016    0.031
##     fb                0.480    0.119    4.027    0.000    1.000    1.000

22.5 PART V: Multilevel CFA 5: Shared + Configural construct

This model specifies the latent variable both at the within and the between level;
The CFA at each level should have the same factor structure, but not necessarily the same parameter estimates;

model5 <- '
  level: 1
    fw =~ a*y1 + b*y2 + c*y3 + d*y4

  level: 2
    fb =~ a*y1 + b*y2 + c*y3 + d*y4   # configural
    fs =~ y1 + y2 + y3 + y4           # shared
    
    # fb and fs must be orthogonal
    fs ~~ 0*fb
'

model5fit <- lavaan::sem(model5, data = Data, 
                         cluster = "clus",
                         estimator = 'MLR',
                         fixed.x = FALSE)

## Warning in lav_object_post_check(object): lavaan WARNING: some estimated ov variances are negative

lavaan::summary(model5fit, fit.measures = T, standardized = T)

## lavaan 0.6-12 ended normally after 56 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        24
##   Number of equality constraints                     3
## 
##   Number of observations                          1000
##   Number of clusters [clus]                        110
## 
## Model Test User Model:
##                                               Standard      Robust
##   Test Statistic                                 0.000       0.000
##   Degrees of freedom                                 3           3
##   P-value (Chi-square)                           1.000       1.000
##   Scaling correction factor                                  0.920
##     Yuan-Bentler correction (Mplus variant)                       
## 
## Model Test Baseline Model:
## 
##   Test statistic                              2619.397    2815.244
##   Degrees of freedom                                12          12
##   P-value                                        0.000       0.000
##   Scaling correction factor                                  0.930
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    1.000       1.000
##   Tucker-Lewis Index (TLI)                       1.005       1.004
##                                                                   
##   Robust Comparative Fit Index (CFI)                         1.000
##   Robust Tucker-Lewis Index (TLI)                            1.004
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -6874.156   -6874.156
##   Scaling correction factor                                  0.905
##       for the MLR correction                                      
##   Loglikelihood unrestricted model (H1)      -6874.338   -6874.338
##   Scaling correction factor                                  1.020
##       for the MLR correction                                      
##                                                                   
##   Akaike (AIC)                               13790.312   13790.312
##   Bayesian (BIC)                             13893.375   13893.375
##   Sample-size adjusted Bayesian (BIC)        13826.678   13826.678
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.000       0.000
##   90 Percent confidence interval - lower         0.000       0.000
##   90 Percent confidence interval - upper         0.000       0.000
##   P-value RMSEA <= 0.05                          1.000       1.000
##                                                                   
##   Robust RMSEA                                               0.000
##   90 Percent confidence interval - lower                     0.000
##   90 Percent confidence interval - upper                     0.000
## 
## Standardized Root Mean Square Residual (corr metric):
## 
##   SRMR (within covariance matrix)                0.001       0.001
##   SRMR (between covariance matrix)               0.003       0.003
## 
## Parameter Estimates:
## 
##   Standard errors                             Sandwich
##   Information bread                           Observed
##   Observed information based on                Hessian
## 
## 
## Level 1 [within]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fw =~                                                                 
##     y1         (a)    1.000                               1.486    0.831
##     y2         (b)    0.993    0.036   27.488    0.000    1.474    0.836
##     y3         (c)    0.992    0.042   23.900    0.000    1.473    0.818
##     y4         (d)    1.007    0.036   28.299    0.000    1.496    0.830
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.000                               0.000    0.000
##    .y2                0.000                               0.000    0.000
##    .y3                0.000                               0.000    0.000
##    .y4                0.000                               0.000    0.000
##     fw                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.989    0.071   13.979    0.000    0.989    0.309
##    .y2                0.934    0.061   15.212    0.000    0.934    0.301
##    .y3                1.076    0.057   18.809    0.000    1.076    0.331
##    .y4                1.014    0.070   14.414    0.000    1.014    0.312
##     fw                2.207    0.167   13.244    0.000    1.000    1.000
## 
## 
## Level 2 [clus]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fb =~                                                                 
##     y1         (a)    1.000                               0.492    0.706
##     y2         (b)    0.993    0.036   27.488    0.000    0.488    0.692
##     y3         (c)    0.992    0.042   23.900    0.000    0.488    0.745
##     y4         (d)    1.007    0.036   28.299    0.000    0.495    0.707
##   fs =~                                                                 
##     y1                1.000                               0.499    0.717
##     y2                1.018    0.288    3.533    0.000    0.508    0.720
##     y3                0.903    1.179    0.766    0.444    0.450    0.688
##     y4                0.960    0.690    1.391    0.164    0.479    0.684
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fb ~~                                                                 
##     fs                0.000                               0.000    0.000
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1               -0.073    0.089   -0.819    0.413   -0.073   -0.104
##    .y2               -0.067    0.088   -0.758    0.449   -0.067   -0.095
##    .y3               -0.035    0.086   -0.406    0.685   -0.035   -0.053
##    .y4               -0.018    0.089   -0.207    0.836   -0.018   -0.026
##     fb                0.000                               0.000    0.000
##     fs                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1               -0.006    0.018   -0.330    0.742   -0.006   -0.012
##    .y2                0.001    0.021    0.063    0.950    0.001    0.003
##    .y3               -0.012    0.027   -0.427    0.670   -0.012   -0.027
##    .y4                0.016    0.021    0.779    0.436    0.016    0.033
##     fb                0.242    3.053    0.079    0.937    1.000    1.000
##     fs                0.249    3.022    0.082    0.934    1.000    1.000

22.6 PART VI: Model Comparison

m1 = fitMeasures(model1fit, fit.measures = c("chisq.scaled", "df.scaled", "pvalue.scaled", "rmsea.robust", "rmsea.ci.lower.robust", "rmsea.ci.upper.robust", "cfi.robust", "tli.robust", "srmr_within", "srmr_between"))
m2 = fitMeasures(model2fit, fit.measures = c("chisq.scaled", "df.scaled", "pvalue.scaled", "rmsea.robust", "rmsea.ci.lower.robust", "rmsea.ci.upper.robust", "cfi.robust", "tli.robust", "srmr_within", "srmr_between"))
m3 = fitMeasures(model3fit, fit.measures = c("chisq.scaled", "df.scaled", "pvalue.scaled", "rmsea.robust", "rmsea.ci.lower.robust", "rmsea.ci.upper.robust", "cfi.robust", "tli.robust", "srmr_within", "srmr_between"))
m4a = fitMeasures(model4afit, fit.measures = c("chisq.scaled", "df.scaled", "pvalue.scaled", "rmsea.robust", "rmsea.ci.lower.robust", "rmsea.ci.upper.robust", "cfi.robust", "tli.robust", "srmr_within", "srmr_between"))
m4b = fitMeasures(model4bfit, fit.measures = c("chisq.scaled", "df.scaled", "pvalue.scaled", "rmsea.robust", "rmsea.ci.lower.robust", "rmsea.ci.upper.robust", "cfi.robust", "tli.robust", "srmr_within", "srmr_between"))
m5 = fitMeasures(model5fit, fit.measures = c("chisq.scaled", "df.scaled", "pvalue.scaled", "rmsea.robust", "rmsea.ci.lower.robust", "rmsea.ci.upper.robust", "cfi.robust", "tli.robust", "srmr_within", "srmr_between"))

#install.packages('qpcR')
library(qpcR)
round(qpcR:::cbind.na(m1, m2, m3, m4a, m4b, m5), 3)

##                           m1    m2      m3   m4a   m4b    m5
## chisq.scaled          64.738 0.107 342.869 0.000 0.059 0.000
## df.scaled              3.000 2.000   3.000 4.000 7.000 3.000
## pvalue.scaled          0.000 0.948   0.000 1.000 1.000 1.000
## rmsea.robust           0.129 0.000   0.456 0.000 0.000 0.000
## rmsea.ci.lower.robust  0.103 0.000   0.416 0.000 0.000 0.000
## rmsea.ci.upper.robust  0.157 0.006   0.498 0.000 0.000 0.000
## cfi.robust             0.981 1.000   0.760 1.000 1.000 1.000
## tli.robust             0.923 1.002   0.042 1.004 1.004 1.004
## srmr_within            0.022    NA   0.266 0.001 0.002 0.001
## srmr_between           0.334    NA   0.010 0.003 0.003 0.003

The final model goes to (drumroll)…model4b!

22.7 PART VII: Adding Covariates to Multilevel SEM

22.7.1 Model A: Adding a within-only covariate

model4wCovA <- '
  level: 1
    fw =~ a*y1 + b*y2 + c*y3 + d*y4
    fw ~ x1

  level: 2
    fb =~ a*y1 + b*y2 + c*y3 + d*y4
'

model4wCovAfit <- lavaan::sem(model4wCovA, data = Data, 
                              cluster = "clus", 
                              estimator = 'MLR',
                              fixed.x = FALSE)

## Warning in lav_object_post_check(object): lavaan WARNING: some estimated ov variances are negative

lavaan::summary(model4wCovAfit, fit.measures = T, standardized = T)

## lavaan 0.6-12 ended normally after 54 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        23
##   Number of equality constraints                     3
## 
##   Number of observations                          1000
##   Number of clusters [clus]                        110
## 
## Model Test User Model:
##                                               Standard      Robust
##   Test Statistic                                 1.793       1.874
##   Degrees of freedom                                10          10
##   P-value (Chi-square)                           0.998       0.997
##   Scaling correction factor                                  0.957
##     Yuan-Bentler correction (Mplus variant)                       
## 
## Model Test Baseline Model:
## 
##   Test statistic                              3085.335    3318.855
##   Degrees of freedom                                16          16
##   P-value                                        0.000       0.000
##   Scaling correction factor                                  0.930
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    1.000       1.000
##   Tucker-Lewis Index (TLI)                       1.004       1.004
##                                                                   
##   Robust Comparative Fit Index (CFI)                         1.000
##   Robust Tucker-Lewis Index (TLI)                            1.004
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -8053.772   -8053.772
##   Scaling correction factor                                  0.901
##       for the MLR correction                                      
##   Loglikelihood unrestricted model (H1)      -8052.875   -8052.875
##   Scaling correction factor                                  1.009
##       for the MLR correction                                      
##                                                                   
##   Akaike (AIC)                               16147.544   16147.544
##   Bayesian (BIC)                             16245.699   16245.699
##   Sample-size adjusted Bayesian (BIC)        16182.178   16182.178
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.000       0.000
##   90 Percent confidence interval - lower         0.000       0.000
##   90 Percent confidence interval - upper         0.000       0.000
##   P-value RMSEA <= 0.05                          1.000       1.000
##                                                                   
##   Robust RMSEA                                               0.000
##   90 Percent confidence interval - lower                     0.000
##   90 Percent confidence interval - upper                     0.000
## 
## Standardized Root Mean Square Residual (corr metric):
## 
##   SRMR (within covariance matrix)                0.004       0.004
##   SRMR (between covariance matrix)               0.004       0.004
## 
## Parameter Estimates:
## 
##   Standard errors                             Sandwich
##   Information bread                           Observed
##   Observed information based on                Hessian
## 
## 
## Level 1 [within]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fw =~                                                                 
##     y1         (a)    1.000                               1.501    0.834
##     y2         (b)    0.994    0.031   32.270    0.000    1.491    0.838
##     y3         (c)    0.984    0.036   27.550    0.000    1.476    0.818
##     y4         (d)    1.006    0.029   34.896    0.000    1.510    0.834
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fw ~                                                                  
##     x1                1.002    0.050   20.044    0.000    0.667    0.663
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.000                               0.000    0.000
##    .y2                0.000                               0.000    0.000
##    .y3                0.000                               0.000    0.000
##    .y4                0.000                               0.000    0.000
##     x1                0.007    0.036    0.195    0.846    0.007    0.007
##    .fw                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.989    0.065   15.126    0.000    0.989    0.305
##    .y2                0.942    0.055   17.114    0.000    0.942    0.297
##    .y3                1.080    0.058   18.491    0.000    1.080    0.331
##    .y4                1.001    0.065   15.438    0.000    1.001    0.305
##    .fw                1.264    0.096   13.119    0.000    0.561    0.561
##     x1                0.985    0.042   23.406    0.000    0.985    1.000
## 
## 
## Level 2 [clus]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fb =~                                                                 
##     y1         (a)    1.000                               0.636    1.004
##     y2         (b)    0.994    0.031   32.270    0.000    0.632    0.998
##     y3         (c)    0.984    0.036   27.550    0.000    0.626    1.016
##     y4         (d)    1.006    0.029   34.896    0.000    0.640    0.982
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1               -0.070    0.080   -0.873    0.383   -0.070   -0.110
##    .y2               -0.064    0.078   -0.823    0.410   -0.064   -0.101
##    .y3               -0.033    0.076   -0.430    0.667   -0.033   -0.053
##    .y4               -0.016    0.078   -0.203    0.839   -0.016   -0.024
##     fb                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1               -0.003    0.017   -0.177    0.860   -0.003   -0.007
##    .y2                0.002    0.019    0.094    0.925    0.002    0.005
##    .y3               -0.012    0.020   -0.624    0.533   -0.012   -0.033
##    .y4                0.015    0.020    0.766    0.444    0.015    0.036
##     fb                0.405    0.085    4.751    0.000    1.000    1.000

22.7.2 Model B: Adding a between-only covariate

model4wCovB <- '
  level: 1
    fw =~ a*y1 + b*y2 + c*y3 + d*y4

  level: 2
    fb =~ a*y1 + b*y2 + c*y3 + d*y4
    fb ~ w
'

model4wCovBfit <- lavaan::sem(model4wCovB, data = Data, 
                              cluster = "clus", 
                              estimator = 'MLR',
                              fixed.x = FALSE)

## Warning in lav_object_post_check(object): lavaan WARNING: some estimated ov variances are negative

lavaan::summary(model4wCovBfit, fit.measures = T, standardized = T)

## lavaan 0.6-12 ended normally after 35 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        23
##   Number of equality constraints                     3
## 
##   Number of observations                          1000
##   Number of clusters [clus]                        110
## 
## Model Test User Model:
##                                               Standard      Robust
##   Test Statistic                                 2.509       2.568
##   Degrees of freedom                                10          10
##   P-value (Chi-square)                           0.991       0.990
##   Scaling correction factor                                  0.977
##     Yuan-Bentler correction (Mplus variant)                       
## 
## Model Test Baseline Model:
## 
##   Test statistic                              2638.435    2788.623
##   Degrees of freedom                                16          16
##   P-value                                        0.000       0.000
##   Scaling correction factor                                  0.946
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    1.000       1.000
##   Tucker-Lewis Index (TLI)                       1.005       1.004
##                                                                   
##   Robust Comparative Fit Index (CFI)                         1.000
##   Robust Tucker-Lewis Index (TLI)                            1.004
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -7010.919   -7010.919
##   Scaling correction factor                                  0.897
##       for the MLR correction                                      
##   Loglikelihood unrestricted model (H1)      -7009.664   -7009.664
##   Scaling correction factor                                  1.013
##       for the MLR correction                                      
##                                                                   
##   Akaike (AIC)                               14061.838   14061.838
##   Bayesian (BIC)                             14159.993   14159.993
##   Sample-size adjusted Bayesian (BIC)        14096.472   14096.472
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.000       0.000
##   90 Percent confidence interval - lower         0.000       0.000
##   90 Percent confidence interval - upper         0.000       0.000
##   P-value RMSEA <= 0.05                          1.000       1.000
##                                                                   
##   Robust RMSEA                                               0.000
##   90 Percent confidence interval - lower                     0.000
##   90 Percent confidence interval - upper                     0.000
## 
## Standardized Root Mean Square Residual (corr metric):
## 
##   SRMR (within covariance matrix)                0.002       0.002
##   SRMR (between covariance matrix)               0.017       0.017
## 
## Parameter Estimates:
## 
##   Standard errors                             Sandwich
##   Information bread                           Observed
##   Observed information based on                Hessian
## 
## 
## Level 1 [within]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fw =~                                                                 
##     y1         (a)    1.000                               1.491    0.832
##     y2         (b)    0.994    0.031   31.642    0.000    1.482    0.838
##     y3         (c)    0.982    0.036   27.079    0.000    1.464    0.815
##     y4         (d)    1.002    0.029   34.699    0.000    1.494    0.829
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.000                               0.000    0.000
##    .y2                0.000                               0.000    0.000
##    .y3                0.000                               0.000    0.000
##    .y4                0.000                               0.000    0.000
##     fw                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.987    0.069   14.241    0.000    0.987    0.308
##    .y2                0.929    0.061   15.165    0.000    0.929    0.297
##    .y3                1.081    0.057   18.940    0.000    1.081    0.335
##    .y4                1.015    0.068   14.894    0.000    1.015    0.313
##     fw                2.223    0.159   14.004    0.000    1.000    1.000
## 
## 
## Level 2 [clus]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fb =~                                                                 
##     y1         (a)    1.000                               0.687    1.007
##     y2         (b)    0.994    0.031   31.642    0.000    0.683    0.994
##     y3         (c)    0.982    0.036   27.079    0.000    0.675    1.015
##     y4         (d)    1.002    0.029   34.699    0.000    0.688    0.983
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fb ~                                                                  
##     w                 0.363    0.103    3.533    0.000    0.529    0.477
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1               -0.073    0.081   -0.900    0.368   -0.073   -0.107
##    .y2               -0.068    0.084   -0.813    0.416   -0.068   -0.099
##    .y3               -0.036    0.079   -0.451    0.652   -0.036   -0.054
##    .y4               -0.019    0.083   -0.229    0.819   -0.019   -0.027
##     w                 0.006    0.086    0.070    0.944    0.006    0.007
##    .fb                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1               -0.007    0.016   -0.415    0.678   -0.007   -0.015
##    .y2                0.006    0.019    0.290    0.772    0.006    0.012
##    .y3               -0.013    0.020   -0.639    0.523   -0.013   -0.029
##    .y4                0.016    0.021    0.778    0.437    0.016    0.033
##    .fb                0.364    0.078    4.645    0.000    0.772    0.772
##     w                 0.815    0.107    7.638    0.000    0.815    1.000

22.7.3 Model C: Adding a covariate at both levels

model4wCovC <- '
  level: 1
    fw =~ a*y1 + b*y2 + c*y3 + d*y4
    fw ~ x1

  level: 2
    fb =~ a*y1 + b*y2 + c*y3 + d*y4
    fb ~ x1
'

model4wCovCfit <- lavaan::sem(model4wCovC, data = Data, 
                              cluster = "clus", 
                              estimator = 'MLR',
                              fixed.x = FALSE)

## Warning in lav_object_post_check(object): lavaan WARNING: some estimated ov variances are negative

lavaan::summary(model4wCovCfit, fit.measures = T, standardized = T)

## lavaan 0.6-12 ended normally after 83 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        25
##   Number of equality constraints                     3
## 
##   Number of observations                          1000
##   Number of clusters [clus]                        110
## 
## Model Test User Model:
##                                               Standard      Robust
##   Test Statistic                                 2.676       2.757
##   Degrees of freedom                                13          13
##   P-value (Chi-square)                           0.999       0.999
##   Scaling correction factor                                  0.971
##     Yuan-Bentler correction (Mplus variant)                       
## 
## Model Test Baseline Model:
## 
##   Test statistic                              3086.903    3246.984
##   Degrees of freedom                                20          20
##   P-value                                        0.000       0.000
##   Scaling correction factor                                  0.951
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    1.000       1.000
##   Tucker-Lewis Index (TLI)                       1.005       1.005
##                                                                   
##   Robust Comparative Fit Index (CFI)                         1.000
##   Robust Tucker-Lewis Index (TLI)                            1.005
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -8052.605   -8052.605
##   Scaling correction factor                                  0.898
##       for the MLR correction                                      
##   Loglikelihood unrestricted model (H1)      -8051.267   -8051.267
##   Scaling correction factor                                  1.002
##       for the MLR correction                                      
##                                                                   
##   Akaike (AIC)                               16149.210   16149.210
##   Bayesian (BIC)                             16257.180   16257.180
##   Sample-size adjusted Bayesian (BIC)        16187.307   16187.307
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.000       0.000
##   90 Percent confidence interval - lower         0.000       0.000
##   90 Percent confidence interval - upper         0.000       0.000
##   P-value RMSEA <= 0.05                          1.000       1.000
##                                                                   
##   Robust RMSEA                                               0.000
##   90 Percent confidence interval - lower                     0.000
##   90 Percent confidence interval - upper                     0.000
## 
## Standardized Root Mean Square Residual (corr metric):
## 
##   SRMR (within covariance matrix)                0.004       0.004
##   SRMR (between covariance matrix)               0.023       0.023
## 
## Parameter Estimates:
## 
##   Standard errors                             Sandwich
##   Information bread                           Observed
##   Observed information based on                Hessian
## 
## 
## Level 1 [within]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fw =~                                                                 
##     y1         (a)    1.000                               1.489    0.832
##     y2         (b)    0.994    0.031   32.214    0.000    1.480    0.836
##     y3         (c)    0.984    0.036   27.538    0.000    1.465    0.816
##     y4         (d)    1.006    0.029   34.847    0.000    1.499    0.832
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fw ~                                                                  
##     x1                0.995    0.052   18.985    0.000    0.668    0.656
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.000                               0.000    0.000
##    .y2                0.000                               0.000    0.000
##    .y3                0.000                               0.000    0.000
##    .y4                0.000                               0.000    0.000
##     x1                0.000                               0.000    0.000
##    .fw                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.989    0.065   15.129    0.000    0.989    0.308
##    .y2                0.942    0.055   17.105    0.000    0.942    0.301
##    .y3                1.080    0.058   18.488    0.000    1.080    0.335
##    .y4                1.001    0.065   15.435    0.000    1.001    0.308
##    .fw                1.263    0.097   13.085    0.000    0.570    0.570
##     x1                0.965    0.044   21.724    0.000    0.965    1.000
## 
## 
## Level 2 [clus]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fb =~                                                                 
##     y1         (a)    1.000                               0.686    1.003
##     y2         (b)    0.994    0.031   32.214    0.000    0.682    0.998
##     y3         (c)    0.984    0.036   27.538    0.000    0.675    1.014
##     y4         (d)    1.006    0.029   34.847    0.000    0.690    0.984
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fb ~                                                                  
##     x1                2.069    1.515    1.366    0.172    3.017    0.432
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1               -0.076    0.087   -0.865    0.387   -0.076   -0.110
##    .y2               -0.070    0.084   -0.839    0.402   -0.070   -0.103
##    .y3               -0.038    0.082   -0.470    0.638   -0.038   -0.058
##    .y4               -0.022    0.084   -0.259    0.795   -0.022   -0.031
##     x1                0.005    0.034    0.141    0.888    0.005    0.034
##    .fb                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1               -0.003    0.017   -0.158    0.875   -0.003   -0.006
##    .y2                0.002    0.019    0.089    0.929    0.002    0.004
##    .y3               -0.012    0.020   -0.624    0.533   -0.012   -0.028
##    .y4                0.015    0.020    0.756    0.450    0.015    0.031
##    .fb                0.383    0.101    3.800    0.000    0.813    0.813
##     x1                0.021    0.016    1.267    0.205    0.021    1.000

22.8 PART VII: Final Model

modelFinal <- '
  level: 1
    fw =~ a*y1 + b*y2 + c*y3 + d*y4
    fw ~ x1 + x2

  level: 2
    fb =~ a*y1 + b*y2 + c*y3 + d*y4
    fb ~ w
'

modelFinalfit <- lavaan::sem(modelFinal, data = Data, 
                              cluster = "clus", 
                              estimator = 'MLR',
                              fixed.x = FALSE)

## Warning in lav_object_post_check(object): lavaan WARNING: some estimated ov variances are negative

lavaan::summary(modelFinalfit, fit.measures = T, standardized = T)

## lavaan 0.6-12 ended normally after 48 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        30
##   Number of equality constraints                     3
## 
##   Number of observations                          1000
##   Number of clusters [clus]                        110
## 
## Model Test User Model:
##                                               Standard      Robust
##   Test Statistic                                 3.741       3.913
##   Degrees of freedom                                16          16
##   P-value (Chi-square)                           0.999       0.999
##   Scaling correction factor                                  0.956
##     Yuan-Bentler correction (Mplus variant)                       
## 
## Model Test Baseline Model:
## 
##   Test statistic                              3283.563    3481.228
##   Degrees of freedom                                24          24
##   P-value                                        0.000       0.000
##   Scaling correction factor                                  0.943
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    1.000       1.000
##   Tucker-Lewis Index (TLI)                       1.006       1.005
##                                                                   
##   Robust Comparative Fit Index (CFI)                         1.000
##   Robust Tucker-Lewis Index (TLI)                            1.005
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -9527.367   -9527.367
##   Scaling correction factor                                  0.941
##       for the MLR correction                                      
##   Loglikelihood unrestricted model (H1)      -9525.497   -9525.497
##   Scaling correction factor                                  1.012
##       for the MLR correction                                      
##                                                                   
##   Akaike (AIC)                               19108.735   19108.735
##   Bayesian (BIC)                             19241.244   19241.244
##   Sample-size adjusted Bayesian (BIC)        19155.491   19155.491
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.000       0.000
##   90 Percent confidence interval - lower         0.000       0.000
##   90 Percent confidence interval - upper         0.000       0.000
##   P-value RMSEA <= 0.05                          1.000       1.000
##                                                                   
##   Robust RMSEA                                               0.000
##   90 Percent confidence interval - lower                     0.000
##   90 Percent confidence interval - upper                     0.000
## 
## Standardized Root Mean Square Residual (corr metric):
## 
##   SRMR (within covariance matrix)                0.004       0.004
##   SRMR (between covariance matrix)               0.014       0.014
## 
## Parameter Estimates:
## 
##   Standard errors                             Sandwich
##   Information bread                           Observed
##   Observed information based on                Hessian
## 
## 
## Level 1 [within]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fw =~                                                                 
##     y1         (a)    1.000                               1.498    0.834
##     y2         (b)    0.992    0.031   31.996    0.000    1.485    0.837
##     y3         (c)    0.982    0.035   27.682    0.000    1.470    0.816
##     y4         (d)    1.005    0.029   34.760    0.000    1.505    0.833
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fw ~                                                                  
##     x1                0.980    0.047   20.863    0.000    0.654    0.649
##     x2                0.515    0.039   13.113    0.000    0.344    0.347
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   x1 ~~                                                                 
##     x2                0.032    0.033    0.972    0.331    0.032    0.032
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.000                               0.000    0.000
##    .y2                0.000                               0.000    0.000
##    .y3                0.000                               0.000    0.000
##    .y4                0.000                               0.000    0.000
##     x1                0.007    0.036    0.194    0.846    0.007    0.007
##     x2                0.014    0.038    0.372    0.710    0.014    0.014
##    .fw                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                0.983    0.063   15.631    0.000    0.983    0.305
##    .y2                0.944    0.056   16.896    0.000    0.944    0.300
##    .y3                1.084    0.057   18.963    0.000    1.084    0.334
##    .y4                1.001    0.064   15.764    0.000    1.001    0.307
##    .fw                0.995    0.081   12.358    0.000    0.444    0.444
##     x1                0.985    0.042   23.406    0.000    0.985    1.000
##     x2                1.017    0.049   20.603    0.000    1.017    1.000
## 
## 
## Level 2 [clus]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fb =~                                                                 
##     y1         (a)    1.000                               0.653    1.005
##     y2         (b)    0.992    0.031   31.996    0.000    0.647    0.996
##     y3         (c)    0.982    0.035   27.682    0.000    0.641    1.017
##     y4         (d)    1.005    0.029   34.760    0.000    0.656    0.982
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fb ~                                                                  
##     w                 0.333    0.076    4.389    0.000    0.510    0.461
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1               -0.083    0.072   -1.141    0.254   -0.083   -0.127
##    .y2               -0.077    0.073   -1.064    0.287   -0.077   -0.119
##    .y3               -0.045    0.069   -0.654    0.513   -0.045   -0.072
##    .y4               -0.029    0.072   -0.400    0.689   -0.029   -0.043
##     w                 0.006    0.086    0.070    0.944    0.006    0.007
##    .fb                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1               -0.004    0.016   -0.246    0.806   -0.004   -0.010
##    .y2                0.003    0.019    0.177    0.859    0.003    0.008
##    .y3               -0.013    0.020   -0.660    0.509   -0.013   -0.033
##    .y4                0.016    0.020    0.778    0.436    0.016    0.036
##    .fb                0.336    0.061    5.520    0.000    0.788    0.788
##     w                 0.815    0.107    7.638    0.000    0.815    1.000