Longitudinal Factor Analysis - Measurement Invariance (WISC)
Overview
In this tutorial we walk through the very basics of testing measurement invariance in the context of a longitudinal factor model. This is not a comprehensive coverage, just an illustration of the basic process.
Outline
In this session we cover …
A. Introduction to Factorial Invariance
B. Configural Invariance Model
C. Weak Invariance Model
D. Strong Invariance Model
E. Strict Invariance Model
F. Comparisons
Prelim - Loading libraries used in this script.
## Warning: package 'psych' was built under R version 4.3.3Introduction to the Common Factor Model
Recall from last time that the basic factor analysis model can be written as a matrix equation …
\[ \boldsymbol{Y_{i}} = \boldsymbol{\Lambda}\boldsymbol{F_{i}} + \boldsymbol{U_{i}} \]
where \(\boldsymbol{Y_{i}}\) is a \(p\) x 1 vector of observed variable scores, \(\boldsymbol{\Lambda}\) is a p x q matrix of factor loadings, \(\boldsymbol{F_{i}}\) is a \(q\) x 1 vector of common factor scores, and \(\boldsymbol{U_{i}}\) is a p x 1 vector of unique factor scores.
We can rewrite the model in terms of variance-covariance and mean expectations. For example, the expected covariance matrix, \(\boldsymbol{\Sigma} = \boldsymbol{Y}'\boldsymbol{Y}\), becomes
\[ \boldsymbol{\Sigma} = \boldsymbol{\Lambda}\boldsymbol{\Psi}\boldsymbol{\Lambda}' + \boldsymbol{\Theta} \] where \(\boldsymbol{\Sigma}\) is a p x p covariance (or correlation) matrix of the observed variables, \(\boldsymbol{\Lambda}\) is a p x q matrix of factor loadings, \(\boldsymbol{\Psi}\) is a q x q covariance matrix of the latent factor variables, and \(\boldsymbol{\Theta}\) is a diagonal matrix of unique factor variances.
and the expected p x 1 mean vector, $ $ becomes \[ \boldsymbol{\mu} = \boldsymbol{\tau} + \boldsymbol{\Lambda}\boldsymbol{\alpha} \]
where \(\boldsymbol{\tau}\) is a p x 1 vector of manifest variable means, \(\boldsymbol{\Lambda}\) is a p x q matrix of factor loadings, and \(\boldsymbol{\alpha}\) is a q x 1 vector of latent variable means.
We can then extend the model to multipe-group or multiple-occasion settings with a group or occasion-specific subscript, so that (e.g., for the two-occasion case)
\[ \boldsymbol{\Sigma_g} = \boldsymbol{\Lambda_g}\boldsymbol{\Psi_g}\boldsymbol{\Lambda_g}' + \boldsymbol{\Theta_g} \] and \[ \boldsymbol{\mu_g} = \boldsymbol{\tau_g} + \boldsymbol{\Lambda_g}\boldsymbol{\alpha_g} \]
Different levels of measurement invariance are established by testing (or requiring) that various matrices are equal.
Specifically, …
For Configural Invaraince we establish that the structure of
\(\boldsymbol{\Lambda_g}\) is
equivalent across occasions.
For Weak Invariance we additionally establish that the factor loadings are equivalent across occasions, \(\boldsymbol{\Lambda_g} = \boldsymbol{\Lambda}\).
For Strong Invariance we additionally establish that the manifest means are equivalent across occasions, \(\boldsymbol{\tau_g} = \boldsymbol{\tau}\). and
For Strict Invariance we additionally establish that the residual/unique variances are also equivalent across occasions \(\boldsymbol{\Theta_g} = \boldsymbol{\Theta}\).
Measurement Invariance testing is usually conducted within a
Structural Equation Modeling (SEM) framework. Here we illustrate how
this may be done using the lavaan package.
Lots of good information and instruction (including about invariance testing - in a multiple-group setting) can be found on the package website … http://lavaan.ugent.be.
Prelim - Reading in Repeated Measures Data
For this example, we use classic data from Osborne & Sudick
(1972) that has been the basis for many examples in Jack McArdle and
colleagues’ papers, courses, and workshops.
The data were obtained from 204 participants (rows) who completed a
variety of cognitive ability tests that are part of the WISC.
Prelim - Descriptives
Lets have a quick look at the data file and the descriptives.
## id info1 comp1 simi1 voca1 info6 comp6 simi6 voca6
## 1 1 31.287 25.627 22.932 22.215 69.883 44.424 68.045 51.162
## 2 2 13.801 14.787 7.581 15.373 41.871 44.862 33.897 37.741
## 3 3 34.970 34.675 28.052 26.841 60.424 50.260 35.844 55.477
## 4 4 24.795 31.391 8.208 20.197 52.865 42.669 45.802 35.987
## 5 5 25.263 30.263 15.977 35.417 67.368 86.654 72.368 60.417
## 6 6 15.402 23.399 11.453 20.560 46.437 52.956 22.537 47.716
## 7 7 15.380 -1.253 2.318 13.004 53.977 36.341 39.912 35.373
## 8 8 19.883 6.704 14.160 14.868 26.901 33.020 21.679 22.763
## 9 9 12.632 13.847 10.276 23.224 54.737 40.163 36.591 54.803
## 10 10 23.690 25.446 11.416 21.786 45.119 52.232 53.508 53.929## vars n mean sd median trimmed mad min max range skew
## info1 1 204 19.78 6.12 19.13 19.77 5.53 1.04 34.97 33.93 -0.06
## comp1 2 204 21.80 9.74 21.90 21.98 10.32 -1.25 46.49 47.74 -0.13
## simi1 3 204 14.90 7.56 14.54 14.57 7.83 -4.52 37.76 42.28 0.39
## voca1 4 204 20.40 6.29 19.84 20.11 6.52 2.77 39.66 36.89 0.37
## info6 5 204 48.51 12.79 46.34 47.80 13.27 24.76 80.00 55.24 0.47
## comp6 6 204 45.17 12.97 44.53 45.22 12.32 0.19 86.65 86.47 -0.01
## simi6 7 204 41.30 14.52 39.44 40.68 14.42 14.29 77.96 63.67 0.39
## voca6 8 204 44.45 11.05 43.80 44.43 12.48 18.17 69.87 51.70 0.01
## kurtosis se
## info1 0.17 0.43
## comp1 -0.49 0.68
## simi1 -0.14 0.53
## voca1 0.02 0.44
## info6 -0.65 0.90
## comp6 0.65 0.91
## simi6 -0.38 1.02
## voca6 -0.63 0.77## info1 comp1 simi1 voca1 info6 comp6 simi6 voca6
## info1 1.00 0.51 0.54 0.56 0.47 0.36 0.44 0.49
## comp1 0.51 1.00 0.45 0.58 0.40 0.36 0.40 0.46
## simi1 0.54 0.45 1.00 0.44 0.35 0.29 0.33 0.38
## voca1 0.56 0.58 0.44 1.00 0.55 0.44 0.54 0.60
## info6 0.47 0.40 0.35 0.55 1.00 0.62 0.67 0.75
## comp6 0.36 0.36 0.29 0.44 0.62 1.00 0.53 0.70
## simi6 0.44 0.40 0.33 0.54 0.67 0.53 1.00 0.66
## voca6 0.49 0.46 0.38 0.60 0.75 0.70 0.66 1.00
Configural Invariance Model
Model Specification
configural <- '
# Define the latent factors.
verb_1 =~ NA*info1 + lambda1*info1 + voca1 + comp1 + simi1
verb_6 =~ NA*info6 + lambda1*info6 + voca6 + comp6 + simi6
#NA undoes the default identification constriant
#lamda1 names the parameter
# Intercepts
info1 ~ i1*1 #constraint on first intercept
voca1 ~ 1
comp1 ~ 1
simi1 ~ 1
info6 ~ i1*1
voca6 ~ 1
comp6 ~ 1
simi6 ~ 1
# Unique Variances
info1 ~~ info1
voca1 ~~ voca1
comp1 ~~ comp1
simi1 ~~ simi1
info6 ~~ info6
voca6 ~~ voca6
comp6 ~~ comp6
simi6 ~~ simi6
# Latent Variable Means
verb_1 ~ 0*1
verb_6 ~ 1
# Latent Variable Variances and Covariance
verb_1 ~~ 1*verb_1 #identification constraint on variance
verb_6 ~~ verb_6
verb_1 ~~ verb_6
'Model Estimation and Interpretation
#Model fitting
fit_configural <- cfa(configural, data = dat, mimic = "mplus")
summary(fit_configural, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6.16 ended normally after 78 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 27
## Number of equality constraints 2
##
## Number of observations 204
## Number of missing patterns 1
##
## Model Test User Model:
##
## Test statistic 25.968
## Degrees of freedom 19
## P-value (Chi-square) 0.131
##
## Model Test Baseline Model:
##
## Test statistic 847.740
## Degrees of freedom 28
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.991
## Tucker-Lewis Index (TLI) 0.987
##
## Robust Comparative Fit Index (CFI) 0.991
## Robust Tucker-Lewis Index (TLI) 0.987
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -5601.115
## Loglikelihood unrestricted model (H1) -5588.131
##
## Akaike (AIC) 11252.229
## Bayesian (BIC) 11335.182
## Sample-size adjusted Bayesian (SABIC) 11255.975
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.042
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.079
## P-value H_0: RMSEA <= 0.050 0.588
## P-value H_0: RMSEA >= 0.080 0.047
##
## Robust RMSEA 0.042
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.079
## P-value H_0: Robust RMSEA <= 0.050 0.588
## P-value H_0: Robust RMSEA >= 0.080 0.047
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.031
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## verb_1 =~
## info1 (lmb1) 4.451 0.400 11.137 0.000 4.451 0.729
## voca1 5.039 0.396 12.728 0.000 5.039 0.803
## comp1 6.850 0.637 10.745 0.000 6.850 0.705
## simi1 4.590 0.520 8.821 0.000 4.590 0.609
## verb_6 =~
## info6 (lmb1) 4.451 0.400 11.137 0.000 10.752 0.843
## voca6 4.102 0.453 9.057 0.000 9.907 0.899
## comp6 4.006 0.489 8.194 0.000 9.675 0.748
## simi6 4.551 0.545 8.346 0.000 10.992 0.759
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## verb_1 ~~
## verb_6 1.837 0.215 8.558 0.000 0.760 0.760
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .info1 (i1) 19.776 0.427 46.273 0.000 19.776 3.240
## .voca1 20.396 0.439 46.416 0.000 20.396 3.250
## .comp1 21.797 0.680 32.036 0.000 21.797 2.243
## .simi1 14.903 0.528 28.223 0.000 14.903 1.976
## .info6 (i1) 19.776 0.427 46.273 0.000 19.776 1.550
## .voca6 17.970 1.844 9.747 0.000 17.970 1.631
## .comp6 19.317 2.299 8.404 0.000 19.317 1.493
## .simi6 11.922 2.516 4.738 0.000 11.922 0.823
## verb_1 0.000 0.000 0.000
## verb_6 6.455 0.606 10.649 0.000 2.672 2.672
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .info1 17.448 2.240 7.789 0.000 17.448 0.468
## .voca1 13.999 2.133 6.563 0.000 13.999 0.355
## .comp1 47.511 5.754 8.257 0.000 47.511 0.503
## .simi1 35.810 4.031 8.884 0.000 35.810 0.630
## .info6 47.096 6.432 7.322 0.000 47.096 0.289
## .voca6 23.267 4.182 5.564 0.000 23.267 0.192
## .comp6 73.851 8.388 8.805 0.000 73.851 0.441
## .simi6 88.920 10.354 8.588 0.000 88.920 0.424
## verb_1 1.000 1.000 1.000
## verb_6 5.834 1.167 4.997 0.000 1.000 1.000
We can also add the possibility that each of the unique factors are correlated across the two occasions
Configural invariance model with correlated uniquenesses
Model Specification
corrunique <- '
# Define the latent factors.
verb_1 =~ NA*info1 + lambda1*info1 + voca1 + comp1 + simi1
verb_6 =~ NA*info6 + lambda1*info6 + voca6 + comp6 + simi6
# Intercepts
info1 ~ i1*1
voca1 ~ 1
comp1 ~ 1
simi1 ~ 1
info6 ~ i1*1
voca6 ~ 1
comp6 ~ 1
simi6 ~ 1
# Unique Variances and Covariances
info1 ~~ info1
voca1 ~~ voca1
comp1 ~~ comp1
simi1 ~~ simi1
info6 ~~ info6
voca6 ~~ voca6
comp6 ~~ comp6
simi6 ~~ simi6
#Cross-time unique Covariances
info1 ~~ info6
voca1 ~~ voca6
comp1 ~~ comp6
simi1 ~~ simi6
# Latent Variable Means
verb_1 ~ 0*1
verb_6 ~ 1
# Latent Variable Variances and Covariance
verb_1 ~~ 1*verb_1
verb_6 ~~ verb_6
verb_1 ~~ verb_6
'Model Estimation and Interpretation
#Model estimation
fit_corrunique <- cfa(corrunique, data = dat, mimic = "mplus")
summary(fit_corrunique, fit.measures = TRUE, standardized=TRUE)## lavaan 0.6.16 ended normally after 111 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 31
## Number of equality constraints 2
##
## Number of observations 204
## Number of missing patterns 1
##
## Model Test User Model:
##
## Test statistic 24.882
## Degrees of freedom 15
## P-value (Chi-square) 0.052
##
## Model Test Baseline Model:
##
## Test statistic 847.740
## Degrees of freedom 28
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.988
## Tucker-Lewis Index (TLI) 0.977
##
## Robust Comparative Fit Index (CFI) 0.988
## Robust Tucker-Lewis Index (TLI) 0.977
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -5600.572
## Loglikelihood unrestricted model (H1) -5588.131
##
## Akaike (AIC) 11259.143
## Bayesian (BIC) 11355.369
## Sample-size adjusted Bayesian (SABIC) 11263.488
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.057
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.095
## P-value H_0: RMSEA <= 0.050 0.350
## P-value H_0: RMSEA >= 0.080 0.177
##
## Robust RMSEA 0.057
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.095
## P-value H_0: Robust RMSEA <= 0.050 0.350
## P-value H_0: Robust RMSEA >= 0.080 0.177
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.030
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## verb_1 =~
## info1 (lmb1) 4.470 0.401 11.140 0.000 4.470 0.732
## voca1 5.000 0.400 12.501 0.000 5.000 0.797
## comp1 6.868 0.639 10.749 0.000 6.868 0.707
## simi1 4.636 0.524 8.847 0.000 4.636 0.614
## verb_6 =~
## info6 (lmb1) 4.470 0.401 11.140 0.000 10.762 0.844
## voca6 4.090 0.447 9.144 0.000 9.849 0.896
## comp6 4.026 0.487 8.258 0.000 9.693 0.749
## simi6 4.564 0.543 8.406 0.000 10.989 0.759
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .info1 ~~
## .info6 0.980 2.601 0.377 0.706 0.980 0.034
## .voca1 ~~
## .voca6 1.741 2.069 0.841 0.400 1.741 0.094
## .comp1 ~~
## .comp6 -0.223 4.828 -0.046 0.963 -0.223 -0.004
## .simi1 ~~
## .simi6 -1.445 4.484 -0.322 0.747 -1.445 -0.026
## verb_1 ~~
## verb_6 1.817 0.212 8.555 0.000 0.755 0.755
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .info1 (i1) 19.776 0.427 46.267 0.000 19.776 3.239
## .voca1 20.396 0.439 46.428 0.000 20.396 3.251
## .comp1 21.797 0.680 32.039 0.000 21.797 2.243
## .simi1 14.903 0.528 28.207 0.000 14.903 1.975
## .info6 (i1) 19.776 0.427 46.267 0.000 19.776 1.550
## .voca6 18.151 1.859 9.762 0.000 18.151 1.651
## .comp6 19.294 2.301 8.386 0.000 19.294 1.491
## .simi6 11.957 2.515 4.754 0.000 11.957 0.826
## verb_1 0.000 0.000 0.000
## verb_6 6.429 0.603 10.658 0.000 2.670 2.670
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .info1 17.294 2.265 7.636 0.000 17.294 0.464
## .voca1 14.374 2.201 6.530 0.000 14.374 0.365
## .comp1 47.245 5.785 8.166 0.000 47.245 0.500
## .simi1 35.454 4.035 8.786 0.000 35.454 0.623
## .info6 46.958 6.519 7.203 0.000 46.958 0.288
## .voca6 23.858 4.329 5.512 0.000 23.858 0.197
## .comp6 73.487 8.419 8.729 0.000 73.487 0.439
## .simi6 88.773 10.405 8.532 0.000 88.773 0.424
## verb_1 1.000 1.000 1.000
## verb_6 5.797 1.150 5.043 0.000 1.000 1.000
The cross-time correlations of the uniquenesses are very small (-.02 to +.09). Because they only contribute a bit of misfit, and because our pedagogical focus is on factorial invariance we remove them from the models going forward.
Weak Invariance Model
Model Specification (additional constraints that \(\boldsymbol{\Lambda_g} = \boldsymbol{\Lambda}\))
weak <- '
# Define the latent factors (all parameters have names for constraints)
verb_1 =~ NA*info1 + lambda1*info1 + lambda2*voca1 + lambda3*comp1 + lambda4*simi1
verb_6 =~ NA*info6 + lambda1*info6 + lambda2*voca6 + lambda3*comp6 + lambda4*simi6
# Intercepts
info1 ~ i1*1
voca1 ~ 1
comp1 ~ 1
simi1 ~ 1
info6 ~ i1*1
voca6 ~ 1
comp6 ~ 1
simi6 ~ 1
# Unique Variances
info1 ~~ info1
voca1 ~~ voca1
comp1 ~~ comp1
simi1 ~~ simi1
info6 ~~ info6
voca6 ~~ voca6
comp6 ~~ comp6
simi6 ~~ simi6
# Latent Variable Means
verb_1 ~ 0*1
verb_6 ~ 1
# Latent Variable Variances and Covariance
verb_1 ~~ 1*verb_1
verb_6 ~~ verb_6
verb_1 ~~ verb_6
'Model Estimation and Interpretation
#Model Estimation
fit_weak <- cfa(weak, data = dat, mimic = "mplus")
summary(fit_weak, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6.16 ended normally after 54 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 27
## Number of equality constraints 5
##
## Number of observations 204
## Number of missing patterns 1
##
## Model Test User Model:
##
## Test statistic 41.897
## Degrees of freedom 22
## P-value (Chi-square) 0.006
##
## Model Test Baseline Model:
##
## Test statistic 847.740
## Degrees of freedom 28
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.976
## Tucker-Lewis Index (TLI) 0.969
##
## Robust Comparative Fit Index (CFI) 0.976
## Robust Tucker-Lewis Index (TLI) 0.969
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -5609.079
## Loglikelihood unrestricted model (H1) -5588.131
##
## Akaike (AIC) 11262.158
## Bayesian (BIC) 11335.157
## Sample-size adjusted Bayesian (SABIC) 11265.454
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.067
## 90 Percent confidence interval - lower 0.035
## 90 Percent confidence interval - upper 0.097
## P-value H_0: RMSEA <= 0.050 0.173
## P-value H_0: RMSEA >= 0.080 0.253
##
## Robust RMSEA 0.067
## 90 Percent confidence interval - lower 0.035
## 90 Percent confidence interval - upper 0.097
## P-value H_0: Robust RMSEA <= 0.050 0.173
## P-value H_0: Robust RMSEA >= 0.080 0.253
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.076
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## verb_1 =~
## info1 (lmb1) 4.933 0.339 14.562 0.000 4.933 0.776
## voca1 (lmb2) 4.865 0.336 14.492 0.000 4.865 0.782
## comp1 (lmb3) 5.172 0.429 12.052 0.000 5.172 0.575
## simi1 (lmb4) 5.072 0.397 12.762 0.000 5.072 0.653
## verb_6 =~
## info6 (lmb1) 4.933 0.339 14.562 0.000 10.158 0.821
## voca6 (lmb2) 4.865 0.336 14.492 0.000 10.017 0.906
## comp6 (lmb3) 5.172 0.429 12.052 0.000 10.650 0.782
## simi6 (lmb4) 5.072 0.397 12.762 0.000 10.444 0.737
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## verb_1 ~~
## verb_6 1.558 0.136 11.439 0.000 0.757 0.757
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .info1 (i1) 19.776 0.445 44.430 0.000 19.776 3.111
## .voca1 20.396 0.436 46.803 0.000 20.396 3.277
## .comp1 21.797 0.629 34.633 0.000 21.797 2.425
## .simi1 14.903 0.544 27.411 0.000 14.903 1.919
## .info6 (i1) 19.776 0.445 44.430 0.000 19.776 1.599
## .voca6 16.111 1.765 9.130 0.000 16.111 1.456
## .comp6 15.049 2.299 6.547 0.000 15.049 1.104
## .simi6 11.756 2.212 5.315 0.000 11.756 0.830
## verb_1 0.000 0.000 0.000
## verb_6 5.824 0.429 13.591 0.000 2.829 2.829
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .info1 16.079 2.188 7.349 0.000 16.079 0.398
## .voca1 15.075 2.109 7.149 0.000 15.075 0.389
## .comp1 54.055 6.012 8.992 0.000 54.055 0.669
## .simi1 34.578 3.941 8.775 0.000 34.578 0.573
## .info6 49.748 6.494 7.661 0.000 49.748 0.325
## .voca6 22.022 4.032 5.462 0.000 22.022 0.180
## .comp6 72.254 8.353 8.650 0.000 72.254 0.389
## .simi6 91.610 10.487 8.736 0.000 91.610 0.457
## verb_1 1.000 1.000 1.000
## verb_6 4.240 0.546 7.759 0.000 1.000 1.000
Strong Invariance Model
Model Specification (additional constraints on manifest intercepts, \(\boldsymbol{\tau_g} = \boldsymbol{\tau}\))
strong <- '
# Define the latent factors (all parameters have names for constraints)
verb_1 =~ NA*info1 + lambda1*info1 + lambda2*voca1 + lambda3*comp1 + lambda4*simi1
verb_6 =~ NA*info6 + lambda1*info6 + lambda2*voca6 + lambda3*comp6 + lambda4*simi6
# Intercepts (all parameters have names for constraints)
info1 ~ i1*1
voca1 ~ i2*1
comp1 ~ i3*1
simi1 ~ i4*1
info6 ~ i1*1
voca6 ~ i2*1
comp6 ~ i3*1
simi6 ~ i4*1
# Unique Variances
info1 ~~ info1
voca1 ~~ voca1
comp1 ~~ comp1
simi1 ~~ simi1
info6 ~~ info6
voca6 ~~ voca6
comp6 ~~ comp6
simi6 ~~ simi6
# Latent Variable Means
verb_1 ~ 0*1
verb_6 ~ 1
# Latent Variable Variances and Covariance
verb_1 ~~ 1*verb_1
verb_6 ~~ verb_6
verb_1 ~~ verb_6
'Model Estimation & Interpretation
fit_strong <- cfa(strong, data = dat, mimic = "mplus")
summary(fit_strong, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6.16 ended normally after 52 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 27
## Number of equality constraints 8
##
## Number of observations 204
## Number of missing patterns 1
##
## Model Test User Model:
##
## Test statistic 53.723
## Degrees of freedom 25
## P-value (Chi-square) 0.001
##
## Model Test Baseline Model:
##
## Test statistic 847.740
## Degrees of freedom 28
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.965
## Tucker-Lewis Index (TLI) 0.961
##
## Robust Comparative Fit Index (CFI) 0.965
## Robust Tucker-Lewis Index (TLI) 0.961
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -5614.992
## Loglikelihood unrestricted model (H1) -5588.131
##
## Akaike (AIC) 11267.984
## Bayesian (BIC) 11331.028
## Sample-size adjusted Bayesian (SABIC) 11270.830
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.075
## 90 Percent confidence interval - lower 0.047
## 90 Percent confidence interval - upper 0.103
## P-value H_0: RMSEA <= 0.050 0.067
## P-value H_0: RMSEA >= 0.080 0.410
##
## Robust RMSEA 0.075
## 90 Percent confidence interval - lower 0.047
## 90 Percent confidence interval - upper 0.103
## P-value H_0: Robust RMSEA <= 0.050 0.067
## P-value H_0: Robust RMSEA >= 0.080 0.410
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.087
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## verb_1 =~
## info1 (lmb1) 5.270 0.331 15.917 0.000 5.270 0.805
## voca1 (lmb2) 4.547 0.291 15.634 0.000 4.547 0.747
## comp1 (lmb3) 4.525 0.312 14.487 0.000 4.525 0.511
## simi1 (lmb4) 4.960 0.328 15.112 0.000 4.960 0.647
## verb_6 =~
## info6 (lmb1) 5.270 0.331 15.917 0.000 11.249 0.858
## voca6 (lmb2) 4.547 0.291 15.634 0.000 9.706 0.891
## comp6 (lmb3) 4.525 0.312 14.487 0.000 9.660 0.746
## simi6 (lmb4) 4.960 0.328 15.112 0.000 10.589 0.746
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## verb_1 ~~
## verb_6 1.608 0.142 11.316 0.000 0.753 0.753
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .info1 (i1) 19.929 0.457 43.623 0.000 19.929 3.043
## .voca1 (i2) 20.311 0.421 48.252 0.000 20.311 3.339
## .comp1 (i3) 21.459 0.598 35.901 0.000 21.459 2.425
## .simi1 (i4) 14.882 0.529 28.146 0.000 14.882 1.941
## .info6 (i1) 19.929 0.457 43.623 0.000 19.929 1.520
## .voca6 (i2) 20.311 0.421 48.252 0.000 20.311 1.863
## .comp6 (i3) 21.459 0.598 35.901 0.000 21.459 1.658
## .simi6 (i4) 14.882 0.529 28.146 0.000 14.882 1.048
## verb_1 0.000 0.000 0.000
## verb_6 5.337 0.354 15.077 0.000 2.500 2.500
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .info1 15.124 2.238 6.757 0.000 15.124 0.353
## .voca1 16.339 2.114 7.729 0.000 16.339 0.441
## .comp1 57.819 6.150 9.402 0.000 57.819 0.738
## .simi1 34.209 3.869 8.841 0.000 34.209 0.582
## .info6 45.360 6.338 7.157 0.000 45.360 0.264
## .voca6 24.586 3.951 6.223 0.000 24.586 0.207
## .comp6 74.229 8.352 8.887 0.000 74.229 0.443
## .simi6 89.572 10.110 8.860 0.000 89.572 0.444
## verb_1 1.000 1.000 1.000
## verb_6 4.557 0.585 7.791 0.000 1.000 1.000
Strict Invariance Model
Model Specification (additional constraints on manifest intercepts, \(\boldsymbol{\Theta_g} = \boldsymbol{\Theta}\))
strict <- '
# Define the latent factors (all parameters have names for constraints)
verb_1 =~ NA*info1 + lambda1*info1 + lambda2*voca1 + lambda3*comp1 + lambda4*simi1
verb_6 =~ NA*info6 + lambda1*info6 + lambda2*voca6 + lambda3*comp6 + lambda4*simi6
# Intercepts (all parameters have names for constraints)
info1 ~ i1*1
voca1 ~ i2*1
comp1 ~ i3*1
simi1 ~ i4*1
info6 ~ i1*1
voca6 ~ i2*1
comp6 ~ i3*1
simi6 ~ i4*1
# Unique Variances (all parameters have names for constraints)
info1 ~~ u1*info1
voca1 ~~ u2*voca1
comp1 ~~ u3*comp1
simi1 ~~ u4*simi1
info6 ~~ u1*info6
voca6 ~~ u2*voca6
comp6 ~~ u3*comp6
simi6 ~~ u4*simi6
# Latent Variable Means
verb_1 ~ 0*1
verb_6 ~ 1
# Latent Variable Variances and Covariance
verb_1 ~~ 1*verb_1
verb_6 ~~ verb_6
verb_1 ~~ verb_6
'Model Estimation & Interpretation
#Model estimation
fit_strict <- cfa(strict, data = dat, mimic = "mplus")
summary(fit_strict, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6.16 ended normally after 57 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 27
## Number of equality constraints 12
##
## Number of observations 204
## Number of missing patterns 1
##
## Model Test User Model:
##
## Test statistic 134.559
## Degrees of freedom 29
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 847.740
## Degrees of freedom 28
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.871
## Tucker-Lewis Index (TLI) 0.876
##
## Robust Comparative Fit Index (CFI) 0.871
## Robust Tucker-Lewis Index (TLI) 0.876
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -5655.410
## Loglikelihood unrestricted model (H1) -5588.131
##
## Akaike (AIC) 11340.820
## Bayesian (BIC) 11390.592
## Sample-size adjusted Bayesian (SABIC) 11343.067
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.134
## 90 Percent confidence interval - lower 0.111
## 90 Percent confidence interval - upper 0.157
## P-value H_0: RMSEA <= 0.050 0.000
## P-value H_0: RMSEA >= 0.080 1.000
##
## Robust RMSEA 0.134
## 90 Percent confidence interval - lower 0.111
## 90 Percent confidence interval - upper 0.157
## P-value H_0: Robust RMSEA <= 0.050 0.000
## P-value H_0: Robust RMSEA >= 0.080 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.169
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## verb_1 =~
## info1 (lmb1) 5.083 0.364 13.978 0.000 5.083 0.689
## voca1 (lmb2) 4.358 0.317 13.765 0.000 4.358 0.690
## comp1 (lmb3) 4.309 0.333 12.944 0.000 4.309 0.463
## simi1 (lmb4) 4.785 0.356 13.447 0.000 4.785 0.521
## verb_6 =~
## info6 (lmb1) 5.083 0.364 13.978 0.000 11.430 0.906
## voca6 (lmb2) 4.358 0.317 13.765 0.000 9.799 0.906
## comp6 (lmb3) 4.309 0.333 12.944 0.000 9.689 0.761
## simi6 (lmb4) 4.785 0.356 13.447 0.000 10.759 0.808
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## verb_1 ~~
## verb_6 1.812 0.168 10.812 0.000 0.806 0.806
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .info1 (i1) 20.019 0.509 39.339 0.000 20.019 2.712
## .voca1 (i2) 20.313 0.433 46.915 0.000 20.313 3.217
## .comp1 (i3) 21.513 0.621 34.648 0.000 21.513 2.312
## .simi1 (i4) 14.805 0.617 23.988 0.000 14.805 1.613
## .info6 (i1) 20.019 0.509 39.339 0.000 20.019 1.586
## .voca6 (i2) 20.313 0.433 46.915 0.000 20.313 1.879
## .comp6 (i3) 21.513 0.621 34.648 0.000 21.513 1.691
## .simi6 (i4) 14.805 0.617 23.988 0.000 14.805 1.112
## verb_1 0.000 0.000 0.000
## verb_6 5.557 0.415 13.396 0.000 2.471 2.471
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .info1 (u1) 28.657 2.958 9.689 0.000 28.657 0.526
## .voca1 (u2) 20.878 2.142 9.746 0.000 20.878 0.524
## .comp1 (u3) 68.013 5.253 12.946 0.000 68.013 0.786
## .simi1 (u4) 61.387 4.897 12.535 0.000 61.387 0.728
## .info6 (u1) 28.657 2.958 9.689 0.000 28.657 0.180
## .voca6 (u2) 20.878 2.142 9.746 0.000 20.878 0.179
## .comp6 (u3) 68.013 5.253 12.946 0.000 68.013 0.420
## .simi6 (u4) 61.387 4.897 12.535 0.000 61.387 0.347
## verb_1 1.000 1.000 1.000
## verb_6 5.056 0.692 7.304 0.000 1.000 1.000
Evaluation of Invariance Models
Comparison of Model fits side-by-side
Model comparison…
Looking at fit statistics
#Compare model fit statistics
round(cbind(configural=inspect(fit_configural, 'fit.measures'),
weak=inspect(fit_weak, 'fit.measures'),
strong=inspect(fit_strong, 'fit.measures'),
strict=inspect(fit_strict, 'fit.measures')),3)## configural weak strong strict
## npar 25.000 22.000 19.000 15.000
## fmin 0.064 0.103 0.132 0.330
## chisq 25.968 41.897 53.723 134.559
## df 19.000 22.000 25.000 29.000
## pvalue 0.131 0.006 0.001 0.000
## baseline.chisq 847.740 847.740 847.740 847.740
## baseline.df 28.000 28.000 28.000 28.000
## baseline.pvalue 0.000 0.000 0.000 0.000
## cfi 0.991 0.976 0.965 0.871
## tli 0.987 0.969 0.961 0.876
## cfi.robust 0.991 0.976 0.965 0.871
## tli.robust 0.987 0.969 0.961 0.876
## nnfi 0.987 0.969 0.961 0.876
## rfi 0.955 0.937 0.929 NA
## nfi 0.969 0.951 0.937 NA
## pnfi 0.658 0.747 0.836 0.871
## ifi 0.992 0.976 0.965 0.871
## rni 0.991 0.976 0.965 0.871
## nnfi.robust 0.987 0.969 0.961 0.876
## rni.robust 0.991 0.976 0.965 0.871
## logl -5601.115 -5609.079 -5614.992 -5655.410
## unrestricted.logl -5588.131 -5588.131 -5588.131 -5588.131
## aic 11252.229 11262.158 11267.984 11340.820
## bic 11335.182 11335.157 11331.028 11390.592
## ntotal 204.000 204.000 204.000 204.000
## bic2 11255.975 11265.454 11270.830 11343.067
## rmsea 0.042 0.067 0.075 0.134
## rmsea.ci.lower 0.000 0.035 0.047 0.111
## rmsea.ci.upper 0.079 0.097 0.103 0.157
## rmsea.ci.level 0.900 0.900 0.900 0.900
## rmsea.pvalue 0.588 0.173 0.067 0.000
## rmsea.close.h0 0.050 0.050 0.050 0.050
## rmsea.notclose.pvalue 0.047 0.253 0.410 1.000
## rmsea.notclose.h0 0.080 0.080 0.080 0.080
## rmsea.robust 0.042 0.067 0.075 0.134
## rmsea.ci.lower.robust 0.000 0.035 0.047 0.111
## rmsea.ci.upper.robust 0.079 0.097 0.103 0.157
## rmsea.pvalue.robust 0.588 0.173 0.067 0.000
## rmsea.notclose.pvalue.robust 0.047 0.253 0.410 1.000
## rmr 2.860 7.497 7.601 10.063
## rmr_nomean 3.162 8.288 8.402 11.125
## srmr 0.031 0.076 0.087 0.169
## srmr_bentler 0.031 0.072 0.086 0.133
## srmr_bentler_nomean 0.034 0.080 0.094 0.147
## crmr 0.034 0.054 0.069 0.089
## crmr_nomean 0.038 0.061 0.077 0.100
## srmr_mplus 0.031 0.076 0.087 0.169
## srmr_mplus_nomean 0.034 0.065 0.079 0.145
## cn_05 237.800 166.180 143.977 65.519
## cn_01 285.306 197.171 169.273 76.178
## gfi 0.995 0.991 0.989 0.976
## agfi 0.988 0.983 0.981 0.963
## pgfi 0.430 0.496 0.562 0.643
## mfi 0.983 0.952 0.932 0.772
## ecvi 0.372 0.421 0.450 0.807Overall, a Strong Invariance model fits the data reasonably well.
We can also do formal model tests.
##
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
## fit_configural 19 11252 11335 25.968
## fit_weak 22 11262 11335 41.897 15.929 0.14535 3 0.001173 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
## fit_weak 22 11262 11335 41.897
## fit_strong 25 11268 11331 53.723 11.826 0.12009 3 0.008005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
## fit_strong 25 11268 11331 53.723
## fit_strict 29 11341 11391 134.559 80.836 0.30686 4 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In all cases, we reject the null hypothesis that the models fit the same.
Thus we might alternatively conclude, based on the statistical information, that there is NOT measurement invariance. Remember that we often have incentives to go on towards the examination of change. Given our goals, we may accept the strong invariance model as sufficient and carry on with the analysis of factor-level change.
Be careful and Thanks for playing!