Level-Specific Fit

Statistics

Multilevel structural equation modeling

Fit statistics

An exploration of alternative ways to obtain level-specific fit indices in multilevel SEM

Author

Mark Lai

Published

February 16, 2024

Code

library(lavaan)

This is lavaan 0.6-20
lavaan is FREE software! Please report any bugs.

Example 1: Single level model

Code

hs_mod <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
hs_fit <- cfa(hs_mod, data = HolzingerSwineford1939)
# Chi-square
anova(hs_fit)

Chi-Squared Test Statistic (unscaled)

          Df    AIC    BIC  Chisq Chisq diff Df diff Pr(>Chisq)    
Saturated  0                0.000                                  
Model     24 7517.5 7595.3 85.305     85.306      24  8.503e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We can obtain the fit index by substituting the model-implied matrix into the fit function

Code

# Model-implied matrix
sigmahat <- lavInspect(hs_fit, what = "implied")$cov
# Saturated matrix
smat <- hs_fit@h1$implied$cov[[1]]
# Fit function
ml_discrepancy <- function(sigmahat, smat) {
    determinant(sigmahat)$modulus - determinant(smat)$modulus + sum(diag(solve(sigmahat, smat))) - nrow(sigmahat)
}
f1 <- ml_discrepancy(sigmahat, smat)
f1 * lavInspect(hs_fit, what = "nobs")

[1] 85.30552
attr(,"logarithm")
[1] TRUE

Multilevel SEM

Built-in example from lavaan

Code

msem_mod <- "
  level: 1
    fw =~ y1 + y2 + y3
    fw ~ x1 + x2 + x3
  level: 2
    fb =~ y1 + y2 + y3
    fb ~ w1 + w2
"
msem_fit <- sem(model = msem_mod, data = Demo.twolevel,
                cluster = "cluster")
anova(msem_fit)

Chi-Squared Test Statistic (unscaled)

          Df   AIC   BIC Chisq Chisq diff Df diff Pr(>Chisq)
Saturated  0             0.000                              
Model     10 24193 24310 8.092      8.092      10     0.6198

Code

# Model-implied matrix
msem_implied <- lavInspect(msem_fit, what = "implied")
sigmawhat <- msem_implied$within$cov
sigmabhat <- msem_implied$cluster$cov
# Saturated matrix
smat <- msem_fit@h1$implied
swmat <- smat$cov[[1]]
sbmat <- smat$cov[[2]]
# Fit function (within)
fw <- ml_discrepancy(sigmawhat, swmat)
fw * lavInspect(msem_fit, what = "nobs")

[1] 1.429024
attr(,"logarithm")
[1] TRUE

Code

# Fit function (between)
fb <- ml_discrepancy(sigmabhat, sbmat)
fb * lavInspect(msem_fit, what = "nclusters")

[1] 11.45465
attr(,"logarithm")
[1] TRUE

Code

# Partially saturated withtin
msem_psb_mod <- "
  level: 1
    y1 ~~ y2 + y3 + x1 + x2 + x3
    y2 ~~ y3 + x1 + x2 + x3
    y3 ~~ x1 + x2 + x3
    x1 ~~ x2 + x3
    x2 ~~ x3
  level: 2
    fb =~ y1 + y2 + y3
    fb ~ w1 + w2
"
msem_psb_fit <- sem(model = msem_psb_mod,
                    data = Demo.twolevel,
                    cluster = "cluster")
anova(msem_psb_fit)

Chi-Squared Test Statistic (unscaled)

          Df   AIC   BIC  Chisq Chisq diff Df diff Pr(>Chisq)
Saturated  0             0.0000                              
Model      4 45599 45803 6.6562     6.6562       4     0.1552

Code

# Model-implied matrix
psb_implied <- lavInspect(msem_psb_fit, what = "implied")
sigmawhat <- psb_implied$within$cov
sigmabhat <- psb_implied$cluster$cov
# Saturated matrix
smat <- msem_psb_fit@h1$implied
swmat <- smat$cov[[1]]
sbmat <- smat$cov[[2]]
# Sample sizes
n_lv1 <- lavInspect(msem_psb_fit, what = "nobs")
n_lv2 <- lavInspect(msem_psb_fit, what = "nclusters")
# Fit function (within)
fw <- ml_discrepancy(sigmawhat, swmat)
fw * (n_lv1 - n_lv2)

[1] 0.003649577
attr(,"logarithm")
[1] TRUE

Code

# Fit function (between)
fb <- ml_discrepancy(sigmabhat, sbmat)
fb * n_lv2

[1] 11.45071
attr(,"logarithm")
[1] TRUE

Code

sigmabhat2 <- sigmabhat
sigmabhat2[1:3, 1:3] <- sigmabhat2[1:3, 1:3] - sigmawhat[1:3, 1:3] / 12.5
sbmat2 <- sbmat
sbmat2[1:3, 1:3] <- sbmat2[1:3, 1:3] - swmat[1:3, 1:3] / 12.5
ml_discrepancy(sigmabhat2, sbmat2)

[1] 0.1803347
attr(,"logarithm")
[1] TRUE

Code

# Partially saturated within
msem_psw_mod <- "
  level: 1
    fw =~ y1 + y2 + y3
    fw ~ x1 + x2 + x3
  level: 2
    y1 ~~ y2 + y3 + w1 + w2
    y2 ~~ y3 + w1 + w2
    y3 ~~ w1 + w2
    w1 ~~ w2
"
msem_psw_fit <- sem(model = msem_psw_mod,
                    data = Demo.twolevel,
                    cluster = "cluster")
anova(msem_psw_fit)

Chi-Squared Test Statistic (unscaled)

          Df   AIC   BIC  Chisq Chisq diff Df diff Pr(>Chisq)
Saturated  0             0.0000                              
Model      6 25293 25462 1.3514     1.3514       6     0.9687

Code

# Model-implied matrix
psw_implied <- lavInspect(msem_psw_fit, what = "implied")
sigmawhat <- psw_implied$within$cov
sigmabhat <- psw_implied$cluster$cov
# Saturated matrix
smat <- msem_psw_fit@h1$implied
swmat <- smat$cov[[1]]
sbmat <- smat$cov[[2]]
n_lv1 <- lavInspect(msem_psw_fit, what = "nobs")
n_lv2 <- lavInspect(msem_psw_fit, what = "nclusters")
# Fit function (within)
fw <- ml_discrepancy(sigmawhat, swmat)
fw * (n_lv1 - n_lv2)

[1] 1.313822
attr(,"logarithm")
[1] TRUE

Code

# Fit function (between)
fb <- ml_discrepancy(sigmabhat, sbmat)
fb * n_lv2

[1] 0.02658882
attr(,"logarithm")
[1] TRUE

Code

fw * (n_lv1 - n_lv2) + fb * n_lv2

[1] 1.34041
attr(,"logarithm")
[1] TRUE