Code
#install.packages("remotes")
#remotes::install_github("DevPsyLab/petersenlab")#install.packages("remotes")
#remotes::install_github("DevPsyLab/petersenlab")mydata_wide <- read.csv("./data/Bliese-Ployhart-2002-indicators-4.csv")mydata_wide$sad1 <- mydata_wide$JOBSAT11
mydata_wide$sad2 <- mydata_wide$JOBSAT21
mydata_wide$sad3 <- mydata_wide$JOBSAT31
mydata_wide$down1 <- mydata_wide$JOBSAT12
mydata_wide$down2 <- mydata_wide$JOBSAT22
mydata_wide$down3 <- mydata_wide$JOBSAT32
mydata_wide$depressed1 <- mydata_wide$JOBSAT13
mydata_wide$depressed2 <- mydata_wide$JOBSAT23
mydata_wide$depressed3 <- mydata_wide$JOBSAT33The example is adapted from the section on longitudinal measurement invariance.
Two primary options for identifying a latent factor include:
(Option #2 may be preferable in alignment optimization.)
configuralInvarianceStandardizeT1Latent_syntax <- '
# Factor Loadings
latent_t1 =~ NA*load1*sad1 + load12*down1 + load13*depressed1
latent_t2 =~ NA*load1*sad2 + load22*down2 + load23*depressed2
latent_t3 =~ NA*load1*sad3 + load32*down3 + load33*depressed3
# Factor Identification: Standardize Factors at T1
## Fix Factor Means at T1 to Zero
latent_t1 ~ 0*1
## Fix Factor Variances at T1 to One
latent_t1 ~~ 1*latent_t1
# Freely Estimate Factor Means at T2 and T3 (relative to T1)
latent_t2 ~ 1
latent_t3 ~ 1
# Freely Estimate Factor Variances at T2 and T3 (relative to T1)
latent_t2 ~~ latent_t2
latent_t3 ~~ latent_t3
# Fix Intercepts of Indicator 1 Across Time
sad1 ~ inta*1
sad2 ~ inta*1
sad3 ~ inta*1
# Free Intercepts of Remaining Manifest Variables
down1 ~ intb1*1
down2 ~ intb2*1
down3 ~ intb3*1
depressed1 ~ intc1*1
depressed2 ~ intc2*1
depressed3 ~ intc3*1
# Residual Covariances Within Indicator Across Time
sad1 ~~ sad2
sad2 ~~ sad3
sad3 ~~ sad3
down1 ~~ depressed2
down2 ~~ depressed3
down3 ~~ depressed3
depressed1 ~~ depressed2
depressed2 ~~ depressed3
depressed3 ~~ depressed3
'
configuralInvarianceStandardizeT1Latent_fit <- cfa(
configuralInvarianceStandardizeT1Latent_syntax,
data = mydata_wide,
missing = "ML",
estimator = "MLR",
meanstructure = TRUE,
#std.lv = TRUE,
fixed.x = FALSE)
summary(
configuralInvarianceStandardizeT1Latent_fit,
fit.measures = TRUE,
standardized = TRUE,
rsquare = TRUE)lavaan 0.6-21 ended normally after 48 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 41
Number of equality constraints 4
Number of observations 495
Number of missing patterns 1
Model Test User Model:
Standard Scaled
Test Statistic 88.966 94.183
Degrees of freedom 17 17
P-value (Chi-square) 0.000 0.000
Scaling correction factor 0.945
Yuan-Bentler correction (Mplus variant)
Model Test Baseline Model:
Test statistic 2120.731 2053.608
Degrees of freedom 36 36
P-value 0.000 0.000
Scaling correction factor 1.033
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.965 0.962
Tucker-Lewis Index (TLI) 0.927 0.919
Robust Comparative Fit Index (CFI) 0.965
Robust Tucker-Lewis Index (TLI) 0.926
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -6038.130 -6038.130
Scaling correction factor 0.870
for the MLR correction
Loglikelihood unrestricted model (H1) -5993.647 -5993.647
Scaling correction factor 0.958
for the MLR correction
Akaike (AIC) 12150.260 12150.260
Bayesian (BIC) 12305.828 12305.828
Sample-size adjusted Bayesian (SABIC) 12188.389 12188.389
Root Mean Square Error of Approximation:
RMSEA 0.092 0.096
90 Percent confidence interval - lower 0.074 0.077
90 Percent confidence interval - upper 0.112 0.116
P-value H_0: RMSEA <= 0.050 0.000 0.000
P-value H_0: RMSEA >= 0.080 0.872 0.917
Robust RMSEA 0.093
90 Percent confidence interval - lower 0.075
90 Percent confidence interval - upper 0.112
P-value H_0: Robust RMSEA <= 0.050 0.000
P-value H_0: Robust RMSEA >= 0.080 0.888
Standardized Root Mean Square Residual:
SRMR 0.029 0.029
Parameter Estimates:
Standard errors Sandwich
Information bread Observed
Observed information based on Hessian
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
latent_t1 =~
sad1 (lod1) 0.990 0.042 23.391 0.000 0.990 0.814
down1 (ld12) 1.030 0.043 24.231 0.000 1.030 0.822
dprssd1 (ld13) 0.950 0.041 22.988 0.000 0.950 0.794
latent_t2 =~
sad2 (lod1) 0.990 0.042 23.391 0.000 0.962 0.795
down2 (ld22) 1.025 0.061 16.793 0.000 0.995 0.830
dprssd2 (ld23) 0.967 0.063 15.445 0.000 0.939 0.782
latent_t3 =~
sad3 (lod1) 0.990 0.042 23.391 0.000 0.842 0.757
down3 (ld32) 1.016 0.125 8.110 0.000 0.864 0.788
dprssd3 (ld33) 1.128 0.132 8.537 0.000 0.959 0.832
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.sad1 ~~
.sad2 0.104 0.034 3.069 0.002 0.104 0.200
.sad2 ~~
.sad3 0.133 0.032 4.114 0.000 0.133 0.249
.down1 ~~
.depressed2 -0.086 0.034 -2.516 0.012 -0.086 -0.161
.down2 ~~
.depressed3 -0.085 0.033 -2.597 0.009 -0.085 -0.198
.down3 ~~
.depressed3 -0.034 0.093 -0.365 0.715 -0.034 -0.079
.depressed1 ~~
.depressed2 0.067 0.034 1.994 0.046 0.067 0.123
.depressed2 ~~
.depressed3 0.090 0.039 2.319 0.020 0.090 0.187
latent_t1 ~~
latent_t2 0.419 0.058 7.204 0.000 0.432 0.432
latent_t3 0.310 0.064 4.863 0.000 0.365 0.365
latent_t2 ~~
latent_t3 0.356 0.057 6.259 0.000 0.431 0.431
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
ltnt_1 0.000 0.000 0.000
ltnt_2 0.000 0.062 0.000 1.000 0.000 0.000
ltnt_3 0.126 0.060 2.097 0.036 0.149 0.149
.sad1 (inta) 3.319 0.055 60.662 0.000 3.319 2.727
.sad2 (inta) 3.319 0.055 60.662 0.000 3.319 2.744
.sad3 (inta) 3.319 0.055 60.662 0.000 3.319 2.986
.down1 (intb1) 3.311 0.056 58.833 0.000 3.311 2.641
.down2 (intb2) 4.331 0.068 63.566 0.000 4.331 3.611
.down3 (intb3) 3.276 0.068 47.867 0.000 3.276 2.987
.dprss1 (intc1) 3.321 0.054 61.701 0.000 3.321 2.776
.dprss2 (intc2) 2.329 0.067 34.561 0.000 2.329 1.940
.dprss3 (intc3) 3.243 0.074 43.744 0.000 3.243 2.812
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
latent_t1 1.000 1.000 1.000
latent_t2 0.943 0.100 9.409 0.000 1.000 1.000
latent_t3 0.722 0.107 6.761 0.000 1.000 1.000
.sad3 0.527 0.081 6.525 0.000 0.527 0.426
.depressed3 0.410 0.108 3.801 0.000 0.410 0.308
.sad1 0.501 0.055 9.034 0.000 0.501 0.338
.down1 0.510 0.052 9.878 0.000 0.510 0.324
.depressed1 0.530 0.048 11.116 0.000 0.530 0.370
.sad2 0.538 0.053 10.237 0.000 0.538 0.368
.down2 0.448 0.051 8.860 0.000 0.448 0.311
.depressed2 0.560 0.053 10.635 0.000 0.560 0.388
.down3 0.456 0.097 4.694 0.000 0.456 0.379
R-Square:
Estimate
sad3 0.574
depressed3 0.692
sad1 0.662
down1 0.676
depressed1 0.630
sad2 0.632
down2 0.689
depressed2 0.612
down3 0.621
metricInvarianceStandardizeT1Latent_syntax <- '
# Factor Loadings
latent_t1 =~ NA*load1*sad1 + load2*down1 + load3*depressed1
latent_t2 =~ NA*load1*sad2 + load2*down2 + load3*depressed2
latent_t3 =~ NA*load1*sad3 + load2*down3 + load3*depressed3
# Factor Identification: Standardize Factors at T1
## Fix Factor Means at T1 to Zero
latent_t1 ~ 0*1
## Fix Factor Variances at T1 to One
latent_t1 ~~ 1*latent_t1
# Freely Estimate Factor Means at T2 and T3 (relative to T1)
latent_t2 ~ 1
latent_t3 ~ 1
# Freely Estimate Factor Variances at T2 and T3 (relative to T1)
latent_t2 ~~ latent_t2
latent_t3 ~~ latent_t3
# Fix Intercepts of Indicator 1 Across Time
sad1 ~ inta*1
sad2 ~ inta*1
sad3 ~ inta*1
# Free Intercepts of Remaining Manifest Variables
down1 ~ intb1*1
down2 ~ intb2*1
down3 ~ intb3*1
depressed1 ~ intc1*1
depressed2 ~ intc2*1
depressed3 ~ intc3*1
# Residual Covariances Within Indicator Across Time
sad1 ~~ sad2
sad2 ~~ sad3
sad3 ~~ sad3
down1 ~~ depressed2
down2 ~~ depressed3
down3 ~~ depressed3
depressed1 ~~ depressed2
depressed2 ~~ depressed3
depressed3 ~~ depressed3
'
metricInvarianceStandardizeT1Latent_fit <- cfa(
metricInvarianceStandardizeT1Latent_syntax,
data = mydata_wide,
missing = "ML",
estimator = "MLR",
meanstructure = TRUE,
#std.lv = TRUE,
fixed.x = FALSE)
summary(
metricInvarianceStandardizeT1Latent_fit,
fit.measures = TRUE,
standardized = TRUE,
rsquare = TRUE)lavaan 0.6-21 ended normally after 41 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 41
Number of equality constraints 8
Number of observations 495
Number of missing patterns 1
Model Test User Model:
Standard Scaled
Test Statistic 95.528 103.806
Degrees of freedom 21 21
P-value (Chi-square) 0.000 0.000
Scaling correction factor 0.920
Yuan-Bentler correction (Mplus variant)
Model Test Baseline Model:
Test statistic 2120.731 2053.608
Degrees of freedom 36 36
P-value 0.000 0.000
Scaling correction factor 1.033
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.964 0.959
Tucker-Lewis Index (TLI) 0.939 0.930
Robust Comparative Fit Index (CFI) 0.964
Robust Tucker-Lewis Index (TLI) 0.937
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -6041.411 -6041.411
Scaling correction factor 0.791
for the MLR correction
Loglikelihood unrestricted model (H1) -5993.647 -5993.647
Scaling correction factor 0.958
for the MLR correction
Akaike (AIC) 12148.821 12148.821
Bayesian (BIC) 12287.572 12287.572
Sample-size adjusted Bayesian (SABIC) 12182.829 12182.829
Root Mean Square Error of Approximation:
RMSEA 0.085 0.089
90 Percent confidence interval - lower 0.068 0.072
90 Percent confidence interval - upper 0.102 0.107
P-value H_0: RMSEA <= 0.050 0.001 0.000
P-value H_0: RMSEA >= 0.080 0.691 0.817
Robust RMSEA 0.086
90 Percent confidence interval - lower 0.069
90 Percent confidence interval - upper 0.102
P-value H_0: Robust RMSEA <= 0.050 0.000
P-value H_0: Robust RMSEA >= 0.080 0.727
Standardized Root Mean Square Residual:
SRMR 0.035 0.035
Parameter Estimates:
Standard errors Sandwich
Information bread Observed
Observed information based on Hessian
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
latent_t1 =~
sad1 (lod1) 0.986 0.038 25.664 0.000 0.986 0.813
down1 (lod2) 1.001 0.036 27.623 0.000 1.001 0.809
dprssd1 (lod3) 0.983 0.037 26.579 0.000 0.983 0.808
latent_t2 =~
sad2 (lod1) 0.986 0.038 25.664 0.000 0.964 0.797
down2 (lod2) 1.001 0.036 27.623 0.000 0.979 0.821
dprssd2 (lod3) 0.983 0.037 26.579 0.000 0.961 0.792
latent_t3 =~
sad3 (lod1) 0.986 0.038 25.664 0.000 0.874 0.783
down3 (lod2) 1.001 0.036 27.623 0.000 0.888 0.792
dprssd3 (lod3) 0.983 0.037 26.579 0.000 0.871 0.775
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.sad1 ~~
.sad2 0.101 0.034 3.008 0.003 0.101 0.197
.sad2 ~~
.sad3 0.124 0.031 4.077 0.000 0.124 0.245
.down1 ~~
.depressed2 -0.087 0.034 -2.557 0.011 -0.087 -0.162
.down2 ~~
.depressed3 -0.071 0.031 -2.311 0.021 -0.071 -0.148
.down3 ~~
.depressed3 0.014 0.045 0.309 0.757 0.014 0.029
.depressed1 ~~
.depressed2 0.065 0.034 1.938 0.053 0.065 0.123
.depressed2 ~~
.depressed3 0.100 0.037 2.699 0.007 0.100 0.190
latent_t1 ~~
latent_t2 0.422 0.055 7.711 0.000 0.431 0.431
latent_t3 0.330 0.051 6.529 0.000 0.372 0.372
latent_t2 ~~
latent_t3 0.376 0.053 7.103 0.000 0.434 0.434
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
ltnt_1 0.000 0.000 0.000
ltnt_2 0.000 0.062 0.000 1.000 0.000 0.000
ltnt_3 0.127 0.061 2.098 0.036 0.143 0.143
.sad1 (inta) 3.319 0.055 60.662 0.000 3.319 2.737
.sad2 (inta) 3.319 0.055 60.662 0.000 3.319 2.743
.sad3 (inta) 3.319 0.055 60.662 0.000 3.319 2.973
.down1 (intb1) 3.311 0.056 58.833 0.000 3.311 2.674
.down2 (intb2) 4.331 0.067 64.302 0.000 4.331 3.634
.down3 (intb3) 3.277 0.065 50.314 0.000 3.277 2.922
.dprss1 (intc1) 3.321 0.054 61.701 0.000 3.321 2.731
.dprss2 (intc2) 2.329 0.068 34.149 0.000 2.329 1.920
.dprss3 (intc3) 3.261 0.066 49.639 0.000 3.261 2.902
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
latent_t1 1.000 1.000 1.000
latent_t2 0.956 0.071 13.430 0.000 1.000 1.000
latent_t3 0.786 0.068 11.518 0.000 1.000 1.000
.sad3 0.481 0.045 10.608 0.000 0.481 0.386
.depressed3 0.504 0.056 9.022 0.000 0.504 0.399
.sad1 0.498 0.052 9.616 0.000 0.498 0.338
.down1 0.530 0.047 11.394 0.000 0.530 0.346
.depressed1 0.513 0.046 11.088 0.000 0.513 0.347
.sad2 0.535 0.050 10.717 0.000 0.535 0.365
.down2 0.462 0.047 9.740 0.000 0.462 0.325
.depressed2 0.549 0.050 11.013 0.000 0.549 0.373
.down3 0.470 0.057 8.269 0.000 0.470 0.373
R-Square:
Estimate
sad3 0.614
depressed3 0.601
sad1 0.662
down1 0.654
depressed1 0.653
sad2 0.635
down2 0.675
depressed2 0.627
down3 0.627
anova(
configuralInvarianceStandardizeT1Latent_fit,
metricInvarianceStandardizeT1Latent_fit
)scalarInvarianceStandardizeT1Latent_syntax <- '
# Factor Loadings
latent_t1 =~ NA*load1*sad1 + load2*down1 + load3*depressed1
latent_t2 =~ NA*load1*sad2 + load2*down2 + load3*depressed2
latent_t3 =~ NA*load1*sad3 + load2*down3 + load3*depressed3
# Factor Identification: Standardize Factors at T1
## Fix Factor Means at T1 to Zero
latent_t1 ~ 0*1
## Fix Factor Variances at T1 to One
latent_t1 ~~ 1*latent_t1
# Freely Estimate Factor Means at T2 and T3 (relative to T1)
latent_t2 ~ 1
latent_t3 ~ 1
# Freely Estimate Factor Variances at T2 and T3 (relative to T1)
latent_t2 ~~ latent_t2
latent_t3 ~~ latent_t3
# Fix Intercepts Across Time
sad1 ~ inta*1
sad2 ~ inta*1
sad3 ~ inta*1
down1 ~ intb*1
down2 ~ intb*1
down3 ~ intb*1
depressed1 ~ intc*1
depressed2 ~ intc*1
depressed3 ~ intc*1
# Residual Covariances Within Indicator Across Time
sad1 ~~ sad2
sad2 ~~ sad3
sad3 ~~ sad3
down1 ~~ depressed2
down2 ~~ depressed3
down3 ~~ depressed3
depressed1 ~~ depressed2
depressed2 ~~ depressed3
depressed3 ~~ depressed3
'
scalarInvarianceStandardizeT1Latent_fit <- cfa(
scalarInvarianceStandardizeT1Latent_syntax,
data = mydata_wide,
missing = "ML",
estimator = "MLR",
meanstructure = TRUE,
#std.lv = TRUE,
fixed.x = FALSE)
summary(
scalarInvarianceStandardizeT1Latent_fit,
fit.measures = TRUE,
standardized = TRUE,
rsquare = TRUE)lavaan 0.6-21 ended normally after 41 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 41
Number of equality constraints 12
Number of observations 495
Number of missing patterns 1
Model Test User Model:
Standard Scaled
Test Statistic 932.999 996.637
Degrees of freedom 25 25
P-value (Chi-square) 0.000 0.000
Scaling correction factor 0.936
Yuan-Bentler correction (Mplus variant)
Model Test Baseline Model:
Test statistic 2120.731 2053.608
Degrees of freedom 36 36
P-value 0.000 0.000
Scaling correction factor 1.033
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.564 0.518
Tucker-Lewis Index (TLI) 0.373 0.307
Robust Comparative Fit Index (CFI) 0.564
Robust Tucker-Lewis Index (TLI) 0.372
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -6460.146 -6460.146
Scaling correction factor 0.691
for the MLR correction
Loglikelihood unrestricted model (H1) -5993.647 -5993.647
Scaling correction factor 0.958
for the MLR correction
Akaike (AIC) 12978.292 12978.292
Bayesian (BIC) 13100.224 13100.224
Sample-size adjusted Bayesian (SABIC) 13008.178 13008.178
Root Mean Square Error of Approximation:
RMSEA 0.271 0.280
90 Percent confidence interval - lower 0.256 0.265
90 Percent confidence interval - upper 0.286 0.296
P-value H_0: RMSEA <= 0.050 0.000 0.000
P-value H_0: RMSEA >= 0.080 1.000 1.000
Robust RMSEA 0.271
90 Percent confidence interval - lower 0.257
90 Percent confidence interval - upper 0.286
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.175 0.175
Parameter Estimates:
Standard errors Sandwich
Information bread Observed
Observed information based on Hessian
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
latent_t1 =~
sad1 (lod1) 1.102 0.039 27.933 0.000 1.102 0.887
down1 (lod2) 0.873 0.044 19.660 0.000 0.873 0.713
dprssd1 (lod3) 0.900 0.047 18.994 0.000 0.900 0.742
latent_t2 =~
sad2 (lod1) 1.102 0.039 27.933 0.000 1.084 0.911
down2 (lod2) 0.873 0.044 19.660 0.000 0.859 0.587
dprssd2 (lod3) 0.900 0.047 18.994 0.000 0.885 0.605
latent_t3 =~
sad3 (lod1) 1.102 0.039 27.933 0.000 0.973 0.873
down3 (lod2) 0.873 0.044 19.660 0.000 0.771 0.683
dprssd3 (lod3) 0.900 0.047 18.994 0.000 0.795 0.691
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.sad1 ~~
.sad2 0.041 0.039 1.054 0.292 0.041 0.144
.sad2 ~~
.sad3 0.051 0.035 1.441 0.150 0.051 0.192
.down1 ~~
.depressed2 0.199 0.056 3.555 0.000 0.199 0.199
.down2 ~~
.depressed3 0.241 0.054 4.450 0.000 0.241 0.245
.down3 ~~
.depressed3 0.108 0.050 2.185 0.029 0.108 0.158
.depressed1 ~~
.depressed2 -0.099 0.050 -2.001 0.045 -0.099 -0.105
.depressed2 ~~
.depressed3 0.064 0.054 1.188 0.235 0.064 0.066
latent_t1 ~~
latent_t2 0.416 0.056 7.490 0.000 0.423 0.423
latent_t3 0.343 0.049 6.957 0.000 0.388 0.388
latent_t2 ~~
latent_t3 0.357 0.055 6.449 0.000 0.411 0.411
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
ltnt_t1 0.000 0.000 0.000
ltnt_t2 -0.022 0.059 -0.380 0.704 -0.023 -0.023
ltnt_t3 0.072 0.055 1.326 0.185 0.082 0.082
.sad1 (inta) 3.343 0.055 60.503 0.000 3.343 2.690
.sad2 (inta) 3.343 0.055 60.503 0.000 3.343 2.809
.sad3 (inta) 3.343 0.055 60.503 0.000 3.343 3.000
.down1 (intb) 3.549 0.053 67.381 0.000 3.549 2.897
.down2 (intb) 3.549 0.053 67.381 0.000 3.549 2.425
.down3 (intb) 3.549 0.053 67.381 0.000 3.549 3.143
.dprssd1 (intc) 3.093 0.055 55.832 0.000 3.093 2.549
.dprssd2 (intc) 3.093 0.055 55.832 0.000 3.093 2.113
.dprssd3 (intc) 3.093 0.055 55.832 0.000 3.093 2.689
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
latent_t1 1.000 1.000 1.000
latent_t2 0.968 0.088 11.012 0.000 1.000 1.000
latent_t3 0.780 0.068 11.528 0.000 1.000 1.000
.sad3 0.294 0.065 4.524 0.000 0.294 0.237
.depressed3 0.692 0.059 11.691 0.000 0.692 0.523
.sad1 0.330 0.070 4.694 0.000 0.330 0.214
.down1 0.739 0.063 11.642 0.000 0.739 0.492
.depressed1 0.662 0.063 10.524 0.000 0.662 0.450
.sad2 0.241 0.088 2.744 0.006 0.241 0.170
.down2 1.405 0.099 14.119 0.000 1.405 0.656
.depressed2 1.360 0.118 11.535 0.000 1.360 0.634
.down3 0.681 0.064 10.705 0.000 0.681 0.534
R-Square:
Estimate
sad3 0.763
depressed3 0.477
sad1 0.786
down1 0.508
depressed1 0.550
sad2 0.830
down2 0.344
depressed2 0.366
down3 0.466
anova(
scalarInvarianceStandardizeT1Latent_fit,
metricInvarianceStandardizeT1Latent_fit
)partialStrongInvarianceStandardizeT1Latent_syntax <- '
# Factor Loadings
latent_t1 =~ NA*load1*sad1 + load2*down1 + load3*depressed1
latent_t2 =~ NA*load1*sad2 + load2*down2 + load3*depressed2
latent_t3 =~ NA*load1*sad3 + load2*down3 + load33*depressed3
# Factor Identification: Standardize Factors at T1
## Fix Factor Means at T1 to Zero
latent_t1 ~ 0*1
## Fix Factor Variances at T1 to One
latent_t1 ~~ 1*latent_t1
# Freely Estimate Factor Means at T2 and T3 (relative to T1)
latent_t2 ~ 1
latent_t3 ~ 1
# Freely Estimate Factor Variances at T2 and T3 (relative to T1)
latent_t2 ~~ latent_t2
latent_t3 ~~ latent_t3
# Fix Some Intercepts Across Time
sad1 ~ inta*1
sad2 ~ inta*1
sad3 ~ inta*1
down1 ~ intb*1
down2 ~ intb2*1
down3 ~ intb*1
depressed1 ~ intc*1
depressed2 ~ intc2*1
depressed3 ~ intc*1
# Residual Covariances Within Indicator Across Time
sad1 ~~ sad2
sad2 ~~ sad3
sad3 ~~ sad3
down1 ~~ depressed2
down2 ~~ depressed3
down3 ~~ depressed3
depressed1 ~~ depressed2
depressed2 ~~ depressed3
depressed3 ~~ depressed3
'
partialStrongInvarianceStandardizeT1Latent_fit <- cfa(
partialStrongInvarianceStandardizeT1Latent_syntax,
data = mydata_wide,
missing = "ML",
estimator = "MLR",
meanstructure = TRUE,
#std.lv = TRUE,
fixed.x = FALSE)
summary(
partialStrongInvarianceStandardizeT1Latent_fit,
fit.measures = TRUE,
standardized = TRUE,
rsquare = TRUE)lavaan 0.6-21 ended normally after 40 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 41
Number of equality constraints 9
Number of observations 495
Number of missing patterns 1
Model Test User Model:
Standard Scaled
Test Statistic 90.407 98.968
Degrees of freedom 22 22
P-value (Chi-square) 0.000 0.000
Scaling correction factor 0.913
Yuan-Bentler correction (Mplus variant)
Model Test Baseline Model:
Test statistic 2120.731 2053.608
Degrees of freedom 36 36
P-value 0.000 0.000
Scaling correction factor 1.033
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.967 0.962
Tucker-Lewis Index (TLI) 0.946 0.938
Robust Comparative Fit Index (CFI) 0.967
Robust Tucker-Lewis Index (TLI) 0.945
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -6038.850 -6038.850
Scaling correction factor 0.772
for the MLR correction
Loglikelihood unrestricted model (H1) -5993.647 -5993.647
Scaling correction factor 0.958
for the MLR correction
Akaike (AIC) 12141.700 12141.700
Bayesian (BIC) 12276.246 12276.246
Sample-size adjusted Bayesian (SABIC) 12174.677 12174.677
Root Mean Square Error of Approximation:
RMSEA 0.079 0.084
90 Percent confidence interval - lower 0.063 0.067
90 Percent confidence interval - upper 0.097 0.102
P-value H_0: RMSEA <= 0.050 0.002 0.001
P-value H_0: RMSEA >= 0.080 0.493 0.668
Robust RMSEA 0.080
90 Percent confidence interval - lower 0.064
90 Percent confidence interval - upper 0.097
P-value H_0: Robust RMSEA <= 0.050 0.001
P-value H_0: Robust RMSEA >= 0.080 0.522
Standardized Root Mean Square Residual:
SRMR 0.029 0.029
Parameter Estimates:
Standard errors Sandwich
Information bread Observed
Observed information based on Hessian
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
latent_t1 =~
sad1 (lod1) 0.991 0.039 25.578 0.000 0.991 0.814
down1 (lod2) 1.025 0.037 27.496 0.000 1.025 0.819
dprssd1 (lod3) 0.956 0.038 25.293 0.000 0.956 0.796
latent_t2 =~
sad2 (lod1) 0.991 0.039 25.578 0.000 0.965 0.797
down2 (lod2) 1.025 0.037 27.496 0.000 0.998 0.831
dprssd2 (lod3) 0.956 0.038 25.293 0.000 0.932 0.779
latent_t3 =~
sad3 (lod1) 0.991 0.039 25.578 0.000 0.840 0.756
down3 (lod2) 1.025 0.037 27.496 0.000 0.869 0.791
dprssd3 (ld33) 1.128 0.073 15.393 0.000 0.957 0.831
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.sad1 ~~
.sad2 0.103 0.034 3.054 0.002 0.103 0.199
.sad2 ~~
.sad3 0.133 0.032 4.215 0.000 0.133 0.250
.down1 ~~
.depressed2 -0.085 0.034 -2.494 0.013 -0.085 -0.157
.down2 ~~
.depressed3 -0.085 0.032 -2.683 0.007 -0.085 -0.199
.down3 ~~
.depressed3 -0.037 0.053 -0.697 0.486 -0.037 -0.085
.depressed1 ~~
.depressed2 0.067 0.033 2.003 0.045 0.067 0.123
.depressed2 ~~
.depressed3 0.090 0.038 2.390 0.017 0.090 0.188
latent_t1 ~~
latent_t2 0.420 0.054 7.749 0.000 0.431 0.431
latent_t3 0.309 0.049 6.313 0.000 0.364 0.364
latent_t2 ~~
latent_t3 0.356 0.051 7.026 0.000 0.431 0.431
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
ltnt_1 0.000 0.000 0.000
ltnt_2 -0.024 0.060 -0.391 0.696 -0.024 -0.024
ltnt_3 0.083 0.051 1.626 0.104 0.098 0.098
.sad1 (inta) 3.344 0.052 64.836 0.000 3.344 2.747
.sad2 (inta) 3.344 0.052 64.836 0.000 3.344 2.760
.sad3 (inta) 3.344 0.052 64.836 0.000 3.344 3.008
.down1 (intb) 3.318 0.051 64.881 0.000 3.318 2.654
.down2 (intb2) 4.355 0.066 65.915 0.000 4.355 3.626
.down3 (intb) 3.318 0.051 64.881 0.000 3.318 3.021
.dprss1 (intc) 3.309 0.053 62.234 0.000 3.309 2.756
.dprss2 (intc2) 2.355 0.064 36.629 0.000 2.355 1.968
.dprss3 (intc) 3.309 0.053 62.234 0.000 3.309 2.873
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
latent_t1 1.000 1.000 1.000
latent_t2 0.949 0.070 13.531 0.000 1.000 1.000
latent_t3 0.720 0.069 10.494 0.000 1.000 1.000
.sad3 0.529 0.049 10.797 0.000 0.529 0.428
.depressed3 0.411 0.075 5.480 0.000 0.411 0.310
.sad1 0.501 0.052 9.589 0.000 0.501 0.338
.down1 0.513 0.048 10.755 0.000 0.513 0.328
.depressed1 0.527 0.046 11.458 0.000 0.527 0.366
.sad2 0.536 0.051 10.609 0.000 0.536 0.365
.down2 0.446 0.049 9.132 0.000 0.446 0.309
.depressed2 0.564 0.050 11.274 0.000 0.564 0.394
.down3 0.451 0.057 7.973 0.000 0.451 0.374
R-Square:
Estimate
sad3 0.572
depressed3 0.690
sad1 0.662
down1 0.672
depressed1 0.634
sad2 0.635
down2 0.691
depressed2 0.606
down3 0.626
anova(
configuralInvarianceStandardizeT1Latent_fit,
partialStrongInvarianceStandardizeT1Latent_fit
)partialStrongInvarianceGrowthCurve_syntax <- '
# Factor Loadings
latent_t1 =~ NA*load1*sad1 + load2*down1 + load3*depressed1
latent_t2 =~ NA*load1*sad2 + load2*down2 + load3*depressed2
latent_t3 =~ NA*load1*sad3 + load2*down3 + load33*depressed3
# Factor Identification: Standardize Factors at T1
## Fix Factor Means at T1 to Zero
latent_t1 ~ 0*1
## Fix Factor Variances at T1 to One
latent_t1 ~~ 1*latent_t1
# Fix Some Intercepts Across Time
sad1 ~ inta*1
sad2 ~ inta*1
sad3 ~ inta*1
down1 ~ intb*1
down2 ~ intb2*1
down3 ~ intb*1
depressed1 ~ intc*1
depressed2 ~ intc2*1
depressed3 ~ intc*1
# Residual Covariances Within Indicator Across Time
sad1 ~~ sad2
sad2 ~~ sad3
sad3 ~~ sad3
down1 ~~ depressed2
down2 ~~ depressed3
down3 ~~ depressed3
depressed1 ~~ depressed2
depressed2 ~~ depressed3
depressed3 ~~ depressed3
# Linear Growth Model
i =~ 1*latent_t1 + 1*latent_t2 + 1*latent_t3
s =~ 0*latent_t1 + 1*latent_t2 + 2*latent_t3
# Variance-covariances of intercepts and slopes
i ~~ i
s ~~ s
i ~~ s
# Means of level and slope
i ~ 1
s ~ 1
'partialStrongInvarianceGrowthCurve_fit <- cfa(
partialStrongInvarianceGrowthCurve_syntax,
data = mydata_wide,
missing = "ML",
estimator = "MLR",
meanstructure = TRUE,
#std.lv = TRUE,
fixed.x = FALSE)summary(
partialStrongInvarianceGrowthCurve_fit,
fit.measures = TRUE,
standardized = TRUE,
rsquare = TRUE)lavaan 0.6-21 ended normally after 59 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 41
Number of equality constraints 9
Number of observations 495
Number of missing patterns 1
Model Test User Model:
Standard Scaled
Test Statistic 90.911 98.285
Degrees of freedom 22 22
P-value (Chi-square) 0.000 0.000
Scaling correction factor 0.925
Yuan-Bentler correction (Mplus variant)
Model Test Baseline Model:
Test statistic 2120.731 2053.608
Degrees of freedom 36 36
P-value 0.000 0.000
Scaling correction factor 1.033
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.967 0.962
Tucker-Lewis Index (TLI) 0.946 0.938
Robust Comparative Fit Index (CFI) 0.966
Robust Tucker-Lewis Index (TLI) 0.945
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -6039.102 -6039.102
Scaling correction factor 0.766
for the MLR correction
Loglikelihood unrestricted model (H1) -5993.647 -5993.647
Scaling correction factor 0.958
for the MLR correction
Akaike (AIC) 12142.205 12142.205
Bayesian (BIC) 12276.750 12276.750
Sample-size adjusted Bayesian (SABIC) 12175.182 12175.182
Root Mean Square Error of Approximation:
RMSEA 0.080 0.084
90 Percent confidence interval - lower 0.063 0.067
90 Percent confidence interval - upper 0.097 0.102
P-value H_0: RMSEA <= 0.050 0.002 0.001
P-value H_0: RMSEA >= 0.080 0.504 0.656
Robust RMSEA 0.080
90 Percent confidence interval - lower 0.064
90 Percent confidence interval - upper 0.097
P-value H_0: Robust RMSEA <= 0.050 0.001
P-value H_0: Robust RMSEA >= 0.080 0.534
Standardized Root Mean Square Residual:
SRMR 0.030 0.030
Parameter Estimates:
Standard errors Sandwich
Information bread Observed
Observed information based on Hessian
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
latent_t1 =~
sad1 (lod1) 0.679 0.081 8.362 0.000 0.991 0.814
down1 (lod2) 0.702 0.080 8.807 0.000 1.025 0.819
dprssd1 (lod3) 0.655 0.081 8.138 0.000 0.957 0.797
latent_t2 =~
sad2 (lod1) 0.679 0.081 8.362 0.000 0.966 0.796
down2 (lod2) 0.702 0.080 8.807 0.000 0.999 0.831
dprssd2 (lod3) 0.655 0.081 8.138 0.000 0.932 0.779
latent_t3 =~
sad3 (lod1) 0.679 0.081 8.362 0.000 0.839 0.755
down3 (lod2) 0.702 0.080 8.807 0.000 0.868 0.791
dprssd3 (ld33) 0.777 0.095 8.168 0.000 0.961 0.834
i =~
ltnt_t1 1.000 0.729 0.729
ltnt_t2 1.000 0.747 0.747
ltnt_t3 1.000 0.860 0.860
s =~
ltnt_t1 0.000 0.000 0.000
ltnt_t2 1.000 0.288 0.288
ltnt_t3 2.000 0.664 0.664
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.sad1 ~~
.sad2 0.103 0.034 3.059 0.002 0.103 0.200
.sad2 ~~
.sad3 0.133 0.032 4.210 0.000 0.133 0.249
.down1 ~~
.depressed2 -0.085 0.034 -2.494 0.013 -0.085 -0.157
.down2 ~~
.depressed3 -0.086 0.032 -2.695 0.007 -0.086 -0.201
.down3 ~~
.depressed3 -0.039 0.053 -0.739 0.460 -0.039 -0.091
.depressed1 ~~
.depressed2 0.067 0.033 2.002 0.045 0.067 0.123
.depressed2 ~~
.depressed3 0.090 0.038 2.374 0.018 0.090 0.188
i ~~
s -0.237 0.178 -1.330 0.184 -0.543 -0.543
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.ltnt_1 0.000 0.000 0.000
.sad1 (inta) 3.646 0.278 13.092 0.000 3.646 2.995
.sad2 (inta) 3.646 0.278 13.092 0.000 3.646 3.007
.sad3 (inta) 3.646 0.278 13.092 0.000 3.646 3.280
.down1 (intb) 3.637 0.286 12.721 0.000 3.637 2.908
.down2 (intb2) 4.649 0.287 16.201 0.000 4.649 3.870
.down3 (intb) 3.637 0.286 12.721 0.000 3.637 3.315
.dprss1 (intc) 3.634 0.292 12.451 0.000 3.634 3.027
.dprss2 (intc2) 2.627 0.270 9.741 0.000 2.627 2.195
.dprss3 (intc) 3.634 0.292 12.451 0.000 3.634 3.153
i -0.494 0.429 -1.151 0.250 -0.465 -0.465
s 0.081 0.042 1.944 0.052 0.198 0.198
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.latent_t1 1.000 0.469 0.469
.sad3 0.531 0.049 10.844 0.000 0.531 0.430
.depressed3 0.405 0.075 5.382 0.000 0.405 0.305
i 1.131 0.472 2.396 0.017 1.000 1.000
s 0.168 0.129 1.301 0.193 1.000 1.000
.sad1 0.500 0.052 9.607 0.000 0.500 0.338
.down1 0.514 0.048 10.781 0.000 0.514 0.329
.depressed1 0.526 0.046 11.443 0.000 0.526 0.365
.sad2 0.537 0.051 10.591 0.000 0.537 0.366
.down2 0.446 0.049 9.132 0.000 0.446 0.309
.depressed2 0.563 0.050 11.261 0.000 0.563 0.393
.down3 0.451 0.057 7.965 0.000 0.451 0.375
.latent_t2 1.199 0.276 4.345 0.000 0.592 0.592
.latent_t3 0.672 0.216 3.115 0.002 0.440 0.440
R-Square:
Estimate
latent_t1 0.531
sad3 0.570
depressed3 0.695
sad1 0.662
down1 0.671
depressed1 0.635
sad2 0.634
down2 0.691
depressed2 0.607
down3 0.625
latent_t2 0.408
latent_t3 0.560
Parameter estimates:
The example is adapted from the section on longitudinal measurement invariance. Two primary options for identifying a latent factor include:
(Option #1 may be preferable in alignment optimization.)
configuralInvarianceStandardizeAllLatent_syntax <- '
# Factor Loadings
latent_t1 =~ NA*load11*sad1 + load12*down1 + load13*depressed1
latent_t2 =~ NA*load21*sad2 + load22*down2 + load23*depressed2
latent_t3 =~ NA*load31*sad3 + load32*down3 + load33*depressed3
# Factor Identification: Standardize Factors at All Timepoints (using std.lv = TRUE)
# Free Intercepts of Manifest Variables
sad1 ~ inta1*1
sad2 ~ inta2*1
sad3 ~ inta3*1
down1 ~ intb1*1
down2 ~ intb2*1
down3 ~ intb3*1
depressed1 ~ intc1*1
depressed2 ~ intc2*1
depressed3 ~ intc3*1
# Residual Covariances Within Indicator Across Time
sad1 ~~ sad2
sad2 ~~ sad3
sad3 ~~ sad3
down1 ~~ depressed2
down2 ~~ depressed3
down3 ~~ depressed3
depressed1 ~~ depressed2
depressed2 ~~ depressed3
depressed3 ~~ depressed3
'configuralInvarianceStandardizeAllLatent_fit <- cfa(
configuralInvarianceStandardizeAllLatent_syntax,
data = mydata_wide,
missing = "ML",
estimator = "MLR",
meanstructure = TRUE,
std.lv = TRUE,
fixed.x = FALSE)summary(
configuralInvarianceStandardizeAllLatent_fit,
fit.measures = TRUE,
standardized = TRUE,
rsquare = TRUE)lavaan 0.6-21 ended normally after 23 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 37
Number of observations 495
Number of missing patterns 1
Model Test User Model:
Standard Scaled
Test Statistic 88.966 94.183
Degrees of freedom 17 17
P-value (Chi-square) 0.000 0.000
Scaling correction factor 0.945
Yuan-Bentler correction (Mplus variant)
Model Test Baseline Model:
Test statistic 2120.731 2053.608
Degrees of freedom 36 36
P-value 0.000 0.000
Scaling correction factor 1.033
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.965 0.962
Tucker-Lewis Index (TLI) 0.927 0.919
Robust Comparative Fit Index (CFI) 0.965
Robust Tucker-Lewis Index (TLI) 0.926
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -6038.130 -6038.130
Scaling correction factor 0.964
for the MLR correction
Loglikelihood unrestricted model (H1) -5993.647 -5993.647
Scaling correction factor 0.958
for the MLR correction
Akaike (AIC) 12150.260 12150.260
Bayesian (BIC) 12305.828 12305.828
Sample-size adjusted Bayesian (SABIC) 12188.389 12188.389
Root Mean Square Error of Approximation:
RMSEA 0.092 0.096
90 Percent confidence interval - lower 0.074 0.077
90 Percent confidence interval - upper 0.112 0.116
P-value H_0: RMSEA <= 0.050 0.000 0.000
P-value H_0: RMSEA >= 0.080 0.872 0.917
Robust RMSEA 0.093
90 Percent confidence interval - lower 0.075
90 Percent confidence interval - upper 0.112
P-value H_0: Robust RMSEA <= 0.050 0.000
P-value H_0: Robust RMSEA >= 0.080 0.888
Standardized Root Mean Square Residual:
SRMR 0.029 0.029
Parameter Estimates:
Standard errors Sandwich
Information bread Observed
Observed information based on Hessian
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
latent_t1 =~
sad1 (ld11) 0.990 0.042 23.391 0.000 0.990 0.814
down1 (ld12) 1.030 0.043 24.231 0.000 1.030 0.822
dprssd1 (ld13) 0.950 0.041 22.988 0.000 0.950 0.794
latent_t2 =~
sad2 (ld21) 0.962 0.042 22.715 0.000 0.962 0.795
down2 (ld22) 0.995 0.040 24.899 0.000 0.995 0.830
dprssd2 (ld23) 0.939 0.044 21.471 0.000 0.939 0.782
latent_t3 =~
sad3 (ld31) 0.842 0.059 14.292 0.000 0.842 0.757
down3 (ld32) 0.864 0.067 12.938 0.000 0.864 0.788
dprssd3 (ld33) 0.959 0.065 14.858 0.000 0.959 0.832
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.sad1 ~~
.sad2 0.104 0.034 3.069 0.002 0.104 0.200
.sad2 ~~
.sad3 0.133 0.032 4.114 0.000 0.133 0.249
.down1 ~~
.depressed2 -0.086 0.034 -2.516 0.012 -0.086 -0.161
.down2 ~~
.depressed3 -0.085 0.033 -2.597 0.009 -0.085 -0.198
.down3 ~~
.depressed3 -0.034 0.093 -0.365 0.715 -0.034 -0.079
.depressed1 ~~
.depressed2 0.067 0.034 1.994 0.046 0.067 0.123
.depressed2 ~~
.depressed3 0.090 0.039 2.319 0.020 0.090 0.187
latent_t1 ~~
latent_t2 0.432 0.054 8.035 0.000 0.432 0.432
latent_t3 0.365 0.060 6.062 0.000 0.365 0.365
latent_t2 ~~
latent_t3 0.431 0.052 8.316 0.000 0.431 0.431
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.sad1 (inta1) 3.319 0.055 60.662 0.000 3.319 2.727
.sad2 (inta2) 3.319 0.054 61.164 0.000 3.319 2.744
.sad3 (inta3) 3.444 0.050 69.021 0.000 3.444 3.099
.down1 (intb1) 3.311 0.056 58.833 0.000 3.311 2.641
.down2 (intb2) 4.331 0.054 80.256 0.000 4.331 3.611
.down3 (intb3) 3.404 0.049 69.067 0.000 3.404 3.104
.dprss1 (intc1) 3.321 0.054 61.701 0.000 3.321 2.776
.dprss2 (intc2) 2.329 0.053 43.663 0.000 2.329 1.940
.dprss3 (intc3) 3.386 0.052 65.572 0.000 3.386 2.936
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.sad3 0.527 0.081 6.525 0.000 0.527 0.426
.depressed3 0.410 0.108 3.801 0.000 0.410 0.308
.sad1 0.501 0.055 9.034 0.000 0.501 0.338
.down1 0.510 0.052 9.878 0.000 0.510 0.324
.depressed1 0.530 0.048 11.116 0.000 0.530 0.370
.sad2 0.538 0.053 10.237 0.000 0.538 0.368
.down2 0.448 0.051 8.860 0.000 0.448 0.311
.depressed2 0.560 0.053 10.635 0.000 0.560 0.388
.down3 0.456 0.097 4.694 0.000 0.456 0.379
latent_t1 1.000 1.000 1.000
latent_t2 1.000 1.000 1.000
latent_t3 1.000 1.000 1.000
R-Square:
Estimate
sad3 0.574
depressed3 0.692
sad1 0.662
down1 0.676
depressed1 0.630
sad2 0.632
down2 0.689
depressed2 0.612
down3 0.621
Adapted from Lai (2013; code):
# Extract loadings and intercepts for alignment
lam_mat <- lavInspect(configuralInvarianceStandardizeAllLatent_fit, what = "est")$lambda
nu_vec <- lavInspect(configuralInvarianceStandardizeAllLatent_fit, what = "est")$nu
# Put them into T x p matrices
num_items <- 3
num_waves <- 3
lam_config <- crossprod(lam_mat, rep(1, num_waves) %x% diag(num_items))
nu_config <- matrix(nu_vec, nrow = num_waves, ncol = num_items, byrow = TRUE)
# Add indicator names
colnames(lam_config) <- colnames(nu_config) <- c("sad", "down", "depressed")
# Alignment optimization
aligned_pars <- sirt::invariance.alignment(
lambda = lam_config,
nu = nu_config,
fixed = TRUE
)The aligned loadings and intercepts are below:
alignedParameters <- rbind(
"Loadings" = rep(NA, num_items),
`rownames<-`(aligned_pars$lambda.aligned,
paste0("Time", 1:num_waves)),
"Intercepts" = rep(NA, num_items),
`rownames<-`(aligned_pars$nu.aligned,
paste0("Time", 1:num_waves))
)
#opts <- options(knitr.kable.NA = "")
#
#knitr::kable(
# alignedParameters,
# format = "simple",
# digits = 3L
#)
print(alignedParameters, na.print = "") sad down depressed
Loadings
Time1 0.9904750 1.030416 0.9495119
Time2 0.9880715 1.022116 0.9644043
Time3 0.9863225 1.012066 1.1235955
Intercepts
Time1 3.3191919 3.311111 3.3212121
Time2 3.3322764 4.344848 2.3420640
Time3 3.3456313 3.302648 3.2732930
aligned_pars$lambda.aligned sad down depressed
latent_t1 0.9904750 1.030416 0.9495119
latent_t2 0.9880715 1.022116 0.9644043
latent_t3 0.9863225 1.012066 1.1235955
aligned_pars$nu.aligned sad down depressed
G1 3.319192 3.311111 3.321212
G2 3.332276 4.344848 2.342064
G3 3.345631 3.302648 3.273293
Lai (2023) suggests using a “50/30/20” rule of thumb for using \(d_{MACS}\) to evaluate the appropriateness of the alignment optimization method: Alignment optimization within CFA “is trustworthy when (a) no more than 50% of items have one or more \(d_{MACS}\) > .20 and (b) no more than 30% of the pairwise \(d_{MACS}\) > .20. Lai (2023) suggests reporting the range of \(d_{MACS}\) values, the proportion of \(d_{MACS}\) > .20, and the proportion of items with at least one \(d_{MACS}\) > .20.
# Function for dMACS
dmacs <- function(loadings, intercepts, pooled_item_sd,
latent_mean = 0, latent_var = 1) {
dloading <- diff(loadings)
dintercept <- diff(intercepts)
integral <- dintercept^2 + 2 * dintercept * dloading * latent_mean +
dloading^2 * (latent_var + latent_mean^2)
sqrt(integral) / pooled_item_sd
}
# Use item SDs at first time point
item_sds_wave1 <-
apply(
mydata_wide[c("sad1", "down1", "depressed1")],
2, sd
)
dmacs_pairwise <- function(loading_mat, intercept_mat, pooled_item_sd,
latent_mean = 0, latent_var = 1) {
ngroups <- nrow(loading_mat)
pairs <- combn(ngroups, 2)
out <- matrix(NA, nrow = ncol(pairs), ncol = ncol(loading_mat))
for (i in seq_len(ncol(pairs))) {
out[i, ] <- dmacs(loading_mat[pairs[, i], ],
intercepts = intercept_mat[pairs[, i], ],
pooled_item_sd,
latent_mean,
latent_var
)
}
rownames(out) <- apply(pairs, 2, paste, collapse = " vs ")
colnames(out) <- colnames(loading_mat)
out
}
# All pairwise dMACS
dmacs_pairwise(aligned_pars$lambda.aligned,
intercept_mat = aligned_pars$nu.aligned,
pooled_item_sd = item_sds_wave1,
latent_mean = 0,
latent_var = 1
) sad down depressed
1 vs 2 0.01091710 0.82476828 0.8168742
1 vs 3 0.02196274 0.01612141 0.1506171
2 vs 3 0.01105291 0.83153231 0.7880754
For the reference indicator (for which the loadings and intercepts are fixed), it is preferable to select an indicator with a strong factor loading (Lai, 2023).
equivalentConfiguralModelAligned_syntax <- '
# Factor Loadings (Loadings of first indicator fixed to alignment solution)
latent_t1 =~ 0.9904750*load11*sad1 + load12*down1 + load13*depressed1
latent_t2 =~ 0.9880715*load21*sad2 + load22*down2 + load23*depressed2
latent_t3 =~ 0.9863225*load31*sad3 + load32*down3 + load33*depressed3
# Intercepts of Manifest Variables (Intercepts of first indicator fixed to alignment solution)
sad1 ~ 3.319192*inta1*1
sad2 ~ 3.332276*inta2*1
sad3 ~ 3.345631*inta3*1
down1 ~ intb1*1
down2 ~ intb2*1
down3 ~ intb3*1
depressed1 ~ intc1*1
depressed2 ~ intc2*1
depressed3 ~ intc3*1
# Residual Covariances Within Indicator Across Time
sad1 ~~ sad2
sad2 ~~ sad3
sad3 ~~ sad3
down1 ~~ depressed2
down2 ~~ depressed3
down3 ~~ depressed3
depressed1 ~~ depressed2
depressed2 ~~ depressed3
depressed3 ~~ depressed3
# Free latent means
latent_t1 + latent_t2 + latent_t3 ~ NA*1
'
equivalentConfiguralModelAligned_fit <- cfa(
equivalentConfiguralModelAligned_syntax,
data = mydata_wide,
missing = "ML",
estimator = "MLR",
meanstructure = TRUE,
#std.lv = TRUE,
fixed.x = FALSE)Compare model fit (they’re equivalent):
anova(
configuralInvarianceStandardizeAllLatent_fit,
equivalentConfiguralModelAligned_fit
)For the reference indicator (for which the loadings and intercepts are fixed), it is preferable to select an indicator with a strong factor loading (Lai, 2023).
awcGrowthModel_syntax <- '
# Factor Loadings (Loadings of first indicator fixed to alignment solution)
latent_t1 =~ 0.9904750*load11*sad1 + load12*down1 + load13*depressed1
latent_t2 =~ 0.9880715*load21*sad2 + load22*down2 + load23*depressed2
latent_t3 =~ 0.9863225*load31*sad3 + load32*down3 + load33*depressed3
# Intercepts of Manifest Variables (Intercepts of first indicator fixed to alignment solution)
sad1 ~ 3.319192*inta1*1
sad2 ~ 3.332276*inta2*1
sad3 ~ 3.345631*inta3*1
down1 ~ intb1*1
down2 ~ intb2*1
down3 ~ intb3*1
depressed1 ~ intc1*1
depressed2 ~ intc2*1
depressed3 ~ intc3*1
# Residual Covariances Within Indicator Across Time
sad1 ~~ sad2
sad2 ~~ sad3
sad3 ~~ sad3
down1 ~~ depressed2
down2 ~~ depressed3
down3 ~~ depressed3
depressed1 ~~ depressed2
depressed2 ~~ depressed3
depressed3 ~~ depressed3
# Linear Growth Model
i =~ 1*latent_t1 + 1*latent_t2 + 1*latent_t3
s =~ 0*latent_t1 + 1*latent_t2 + 2*latent_t3
# Variance-covariances of intercepts and slopes
i ~~ i
s ~~ s
i ~~ s
# Means of level and slope
i ~ 1
s ~ 1
# Fixed disturbances of latent outcomes to zero
latent_t1 ~ 0*1
latent_t2 ~ 0*1
latent_t3 ~ 0*1
'awcGrowthModel_fit <- cfa(
awcGrowthModel_syntax,
data = mydata_wide,
missing = "ML",
estimator = "MLR",
meanstructure = TRUE,
#std.lv = TRUE,
fixed.x = FALSE)summary(
configuralInvarianceStandardizeAllLatent_fit,
fit.measures = TRUE,
standardized = TRUE,
rsquare = TRUE)lavaan 0.6-21 ended normally after 23 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 37
Number of observations 495
Number of missing patterns 1
Model Test User Model:
Standard Scaled
Test Statistic 88.966 94.183
Degrees of freedom 17 17
P-value (Chi-square) 0.000 0.000
Scaling correction factor 0.945
Yuan-Bentler correction (Mplus variant)
Model Test Baseline Model:
Test statistic 2120.731 2053.608
Degrees of freedom 36 36
P-value 0.000 0.000
Scaling correction factor 1.033
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.965 0.962
Tucker-Lewis Index (TLI) 0.927 0.919
Robust Comparative Fit Index (CFI) 0.965
Robust Tucker-Lewis Index (TLI) 0.926
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -6038.130 -6038.130
Scaling correction factor 0.964
for the MLR correction
Loglikelihood unrestricted model (H1) -5993.647 -5993.647
Scaling correction factor 0.958
for the MLR correction
Akaike (AIC) 12150.260 12150.260
Bayesian (BIC) 12305.828 12305.828
Sample-size adjusted Bayesian (SABIC) 12188.389 12188.389
Root Mean Square Error of Approximation:
RMSEA 0.092 0.096
90 Percent confidence interval - lower 0.074 0.077
90 Percent confidence interval - upper 0.112 0.116
P-value H_0: RMSEA <= 0.050 0.000 0.000
P-value H_0: RMSEA >= 0.080 0.872 0.917
Robust RMSEA 0.093
90 Percent confidence interval - lower 0.075
90 Percent confidence interval - upper 0.112
P-value H_0: Robust RMSEA <= 0.050 0.000
P-value H_0: Robust RMSEA >= 0.080 0.888
Standardized Root Mean Square Residual:
SRMR 0.029 0.029
Parameter Estimates:
Standard errors Sandwich
Information bread Observed
Observed information based on Hessian
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
latent_t1 =~
sad1 (ld11) 0.990 0.042 23.391 0.000 0.990 0.814
down1 (ld12) 1.030 0.043 24.231 0.000 1.030 0.822
dprssd1 (ld13) 0.950 0.041 22.988 0.000 0.950 0.794
latent_t2 =~
sad2 (ld21) 0.962 0.042 22.715 0.000 0.962 0.795
down2 (ld22) 0.995 0.040 24.899 0.000 0.995 0.830
dprssd2 (ld23) 0.939 0.044 21.471 0.000 0.939 0.782
latent_t3 =~
sad3 (ld31) 0.842 0.059 14.292 0.000 0.842 0.757
down3 (ld32) 0.864 0.067 12.938 0.000 0.864 0.788
dprssd3 (ld33) 0.959 0.065 14.858 0.000 0.959 0.832
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.sad1 ~~
.sad2 0.104 0.034 3.069 0.002 0.104 0.200
.sad2 ~~
.sad3 0.133 0.032 4.114 0.000 0.133 0.249
.down1 ~~
.depressed2 -0.086 0.034 -2.516 0.012 -0.086 -0.161
.down2 ~~
.depressed3 -0.085 0.033 -2.597 0.009 -0.085 -0.198
.down3 ~~
.depressed3 -0.034 0.093 -0.365 0.715 -0.034 -0.079
.depressed1 ~~
.depressed2 0.067 0.034 1.994 0.046 0.067 0.123
.depressed2 ~~
.depressed3 0.090 0.039 2.319 0.020 0.090 0.187
latent_t1 ~~
latent_t2 0.432 0.054 8.035 0.000 0.432 0.432
latent_t3 0.365 0.060 6.062 0.000 0.365 0.365
latent_t2 ~~
latent_t3 0.431 0.052 8.316 0.000 0.431 0.431
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.sad1 (inta1) 3.319 0.055 60.662 0.000 3.319 2.727
.sad2 (inta2) 3.319 0.054 61.164 0.000 3.319 2.744
.sad3 (inta3) 3.444 0.050 69.021 0.000 3.444 3.099
.down1 (intb1) 3.311 0.056 58.833 0.000 3.311 2.641
.down2 (intb2) 4.331 0.054 80.256 0.000 4.331 3.611
.down3 (intb3) 3.404 0.049 69.067 0.000 3.404 3.104
.dprss1 (intc1) 3.321 0.054 61.701 0.000 3.321 2.776
.dprss2 (intc2) 2.329 0.053 43.663 0.000 2.329 1.940
.dprss3 (intc3) 3.386 0.052 65.572 0.000 3.386 2.936
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.sad3 0.527 0.081 6.525 0.000 0.527 0.426
.depressed3 0.410 0.108 3.801 0.000 0.410 0.308
.sad1 0.501 0.055 9.034 0.000 0.501 0.338
.down1 0.510 0.052 9.878 0.000 0.510 0.324
.depressed1 0.530 0.048 11.116 0.000 0.530 0.370
.sad2 0.538 0.053 10.237 0.000 0.538 0.368
.down2 0.448 0.051 8.860 0.000 0.448 0.311
.depressed2 0.560 0.053 10.635 0.000 0.560 0.388
.down3 0.456 0.097 4.694 0.000 0.456 0.379
latent_t1 1.000 1.000 1.000
latent_t2 1.000 1.000 1.000
latent_t3 1.000 1.000 1.000
R-Square:
Estimate
sad3 0.574
depressed3 0.692
sad1 0.662
down1 0.676
depressed1 0.630
sad2 0.632
down2 0.689
depressed2 0.612
down3 0.621
Parameter estimates:
| lhs | op | rhs | est | se | z | pvalue | ci.lower | ci.upper | |
|---|---|---|---|---|---|---|---|---|---|
| 34 | i | i | 0.528 | 0.116 | 4.547 | 0.000 | 0.300 | 0.756 | |
| 35 | s | s | 0.077 | 0.056 | 1.384 | 0.166 | -0.032 | 0.186 | |
| 36 | i | s | -0.108 | 0.072 | -1.500 | 0.134 | -0.249 | 0.033 | |
| 37 | i | ~1 | -0.018 | 0.053 | -0.349 | 0.727 | -0.121 | 0.085 | |
| 38 | s | ~1 | 0.052 | 0.031 | 1.710 | 0.087 | -0.008 | 0.112 |
It can be potentially helpful to consider regularization/Bayesian penalty methods (e.g., using priors such as Horseshoe, Lasso, or Spike-and-slab).
For examples of IRT models, see here.
---
title: "Developmental Scaling"
---
# Preamble
## Install Libraries
```{r}
#install.packages("remotes")
#remotes::install_github("DevPsyLab/petersenlab")
```
## Load Libraries
```{r}
library("lavaan")
library("semTools")
library("semPlot")
library("lavaanPlot")
library("lavaangui")
library("sirt")
library("kableExtra")
library("tidyverse")
```
## Import data
```{r}
mydata_wide <- read.csv("./data/Bliese-Ployhart-2002-indicators-4.csv")
```
## Data Processing
```{r}
mydata_wide$sad1 <- mydata_wide$JOBSAT11
mydata_wide$sad2 <- mydata_wide$JOBSAT21
mydata_wide$sad3 <- mydata_wide$JOBSAT31
mydata_wide$down1 <- mydata_wide$JOBSAT12
mydata_wide$down2 <- mydata_wide$JOBSAT22
mydata_wide$down3 <- mydata_wide$JOBSAT32
mydata_wide$depressed1 <- mydata_wide$JOBSAT13
mydata_wide$depressed2 <- mydata_wide$JOBSAT23
mydata_wide$depressed3 <- mydata_wide$JOBSAT33
```
# Approaches
1. [Partial strong longitudinal factorial invariance](#sec-partialStrongInvariance) (paired with second-order growth functions—i.e., a second-order growth curve model)
1. [Alignment optimization](#sec-alignment) (paired with second-order growth functions)
1. [Moderated nonlinear factor analysis](#sec-mnlfa) (MNLFA) (paired with second-order growth functions)
1. [Item response theory (IRT) with anchor items](#sec-irt) (potentially with regularization/Bayesian penalty methods on the priors)
- [separate calibration](#sec-irtSeparateCalibration)
- [concurrent calibration](#sec-irtConcurrentCalibration)
# Partial Strong Longitudinal Factorial Invariance {#sec-partialStrongInvariance}
The example is adapted from the section on [longitudinal measurement invariance](#sec-longitudinalMI).
Two primary options for identifying a latent factor include:
1. standardizing the latent factor(s) at all timepoints (fixing their intercepts to zero and their variances to one)
1. standardizing the latent factor(s) at T1 (fixing their intercepts to zero and their variances to one) and constraining the factor loadings and intercepts of one item to be equal across time
(Option #2 may be preferable in alignment optimization.)
## Standardizing the Latent Factors at All Timepoints
## Standardizing the Latent Factors at T1 while Constraining Loadings and Intercepts of One Item Across Time
### Configural Invariance With Correlated Residuals Within-Indicator Across Time
1. Standardize the latent factor(s) at T1 (i.e., fix the mean to zero and the variance to one)
1. For each latent construct, constrain the first indicator's factor loading to be the same across time
1. For each latent construct, constrain the first indicator's intercept to be the same across time
1. Allow within-indicator residuals to covary across time
```{r}
configuralInvarianceStandardizeT1Latent_syntax <- '
# Factor Loadings
latent_t1 =~ NA*load1*sad1 + load12*down1 + load13*depressed1
latent_t2 =~ NA*load1*sad2 + load22*down2 + load23*depressed2
latent_t3 =~ NA*load1*sad3 + load32*down3 + load33*depressed3
# Factor Identification: Standardize Factors at T1
## Fix Factor Means at T1 to Zero
latent_t1 ~ 0*1
## Fix Factor Variances at T1 to One
latent_t1 ~~ 1*latent_t1
# Freely Estimate Factor Means at T2 and T3 (relative to T1)
latent_t2 ~ 1
latent_t3 ~ 1
# Freely Estimate Factor Variances at T2 and T3 (relative to T1)
latent_t2 ~~ latent_t2
latent_t3 ~~ latent_t3
# Fix Intercepts of Indicator 1 Across Time
sad1 ~ inta*1
sad2 ~ inta*1
sad3 ~ inta*1
# Free Intercepts of Remaining Manifest Variables
down1 ~ intb1*1
down2 ~ intb2*1
down3 ~ intb3*1
depressed1 ~ intc1*1
depressed2 ~ intc2*1
depressed3 ~ intc3*1
# Residual Covariances Within Indicator Across Time
sad1 ~~ sad2
sad2 ~~ sad3
sad3 ~~ sad3
down1 ~~ depressed2
down2 ~~ depressed3
down3 ~~ depressed3
depressed1 ~~ depressed2
depressed2 ~~ depressed3
depressed3 ~~ depressed3
'
configuralInvarianceStandardizeT1Latent_fit <- cfa(
configuralInvarianceStandardizeT1Latent_syntax,
data = mydata_wide,
missing = "ML",
estimator = "MLR",
meanstructure = TRUE,
#std.lv = TRUE,
fixed.x = FALSE)
summary(
configuralInvarianceStandardizeT1Latent_fit,
fit.measures = TRUE,
standardized = TRUE,
rsquare = TRUE)
```
### Metric ("Weak Factorial") Invariance
1. Standardize the latent factor(s) at T1 (i.e., fix the mean to zero and the variance to one)
1. For each latent construct, constrain the first indicator's intercept to be the same across time
1. Allow within-indicator residuals to covary across time
1. **For each indicator, constrain its factor loading to be the same across time**
```{r}
metricInvarianceStandardizeT1Latent_syntax <- '
# Factor Loadings
latent_t1 =~ NA*load1*sad1 + load2*down1 + load3*depressed1
latent_t2 =~ NA*load1*sad2 + load2*down2 + load3*depressed2
latent_t3 =~ NA*load1*sad3 + load2*down3 + load3*depressed3
# Factor Identification: Standardize Factors at T1
## Fix Factor Means at T1 to Zero
latent_t1 ~ 0*1
## Fix Factor Variances at T1 to One
latent_t1 ~~ 1*latent_t1
# Freely Estimate Factor Means at T2 and T3 (relative to T1)
latent_t2 ~ 1
latent_t3 ~ 1
# Freely Estimate Factor Variances at T2 and T3 (relative to T1)
latent_t2 ~~ latent_t2
latent_t3 ~~ latent_t3
# Fix Intercepts of Indicator 1 Across Time
sad1 ~ inta*1
sad2 ~ inta*1
sad3 ~ inta*1
# Free Intercepts of Remaining Manifest Variables
down1 ~ intb1*1
down2 ~ intb2*1
down3 ~ intb3*1
depressed1 ~ intc1*1
depressed2 ~ intc2*1
depressed3 ~ intc3*1
# Residual Covariances Within Indicator Across Time
sad1 ~~ sad2
sad2 ~~ sad3
sad3 ~~ sad3
down1 ~~ depressed2
down2 ~~ depressed3
down3 ~~ depressed3
depressed1 ~~ depressed2
depressed2 ~~ depressed3
depressed3 ~~ depressed3
'
metricInvarianceStandardizeT1Latent_fit <- cfa(
metricInvarianceStandardizeT1Latent_syntax,
data = mydata_wide,
missing = "ML",
estimator = "MLR",
meanstructure = TRUE,
#std.lv = TRUE,
fixed.x = FALSE)
summary(
metricInvarianceStandardizeT1Latent_fit,
fit.measures = TRUE,
standardized = TRUE,
rsquare = TRUE)
anova(
configuralInvarianceStandardizeT1Latent_fit,
metricInvarianceStandardizeT1Latent_fit
)
```
### Scalar ("Strong Factorial") Invariance
1. Standardize the latent factor(s) at T1 (i.e., fix the mean to zero and the variance to one)
1. Allow within-indicator residuals to covary across time
1. For each indicator, constrain its factor loading to be the same across time
1. **For each indicator, constrain its intercept to be the same across time**
```{r}
scalarInvarianceStandardizeT1Latent_syntax <- '
# Factor Loadings
latent_t1 =~ NA*load1*sad1 + load2*down1 + load3*depressed1
latent_t2 =~ NA*load1*sad2 + load2*down2 + load3*depressed2
latent_t3 =~ NA*load1*sad3 + load2*down3 + load3*depressed3
# Factor Identification: Standardize Factors at T1
## Fix Factor Means at T1 to Zero
latent_t1 ~ 0*1
## Fix Factor Variances at T1 to One
latent_t1 ~~ 1*latent_t1
# Freely Estimate Factor Means at T2 and T3 (relative to T1)
latent_t2 ~ 1
latent_t3 ~ 1
# Freely Estimate Factor Variances at T2 and T3 (relative to T1)
latent_t2 ~~ latent_t2
latent_t3 ~~ latent_t3
# Fix Intercepts Across Time
sad1 ~ inta*1
sad2 ~ inta*1
sad3 ~ inta*1
down1 ~ intb*1
down2 ~ intb*1
down3 ~ intb*1
depressed1 ~ intc*1
depressed2 ~ intc*1
depressed3 ~ intc*1
# Residual Covariances Within Indicator Across Time
sad1 ~~ sad2
sad2 ~~ sad3
sad3 ~~ sad3
down1 ~~ depressed2
down2 ~~ depressed3
down3 ~~ depressed3
depressed1 ~~ depressed2
depressed2 ~~ depressed3
depressed3 ~~ depressed3
'
scalarInvarianceStandardizeT1Latent_fit <- cfa(
scalarInvarianceStandardizeT1Latent_syntax,
data = mydata_wide,
missing = "ML",
estimator = "MLR",
meanstructure = TRUE,
#std.lv = TRUE,
fixed.x = FALSE)
summary(
scalarInvarianceStandardizeT1Latent_fit,
fit.measures = TRUE,
standardized = TRUE,
rsquare = TRUE)
anova(
scalarInvarianceStandardizeT1Latent_fit,
metricInvarianceStandardizeT1Latent_fit
)
```
### Partial Strong Factorial Invariance
1. Standardize the latent factor(s) at T1 (i.e., fix the mean to zero and the variance to one)
1. Allow within-indicator residuals to covary across time
1. For some indicators (i.e., the indicators with invariant factor loadings, called "anchor items"), constrain its factor loading to be the same across time
1. For some indicators (i.e., the indicators with invariant intercepts, called "anchor items"), constrain its intercept to be the same across time
```{r}
partialStrongInvarianceStandardizeT1Latent_syntax <- '
# Factor Loadings
latent_t1 =~ NA*load1*sad1 + load2*down1 + load3*depressed1
latent_t2 =~ NA*load1*sad2 + load2*down2 + load3*depressed2
latent_t3 =~ NA*load1*sad3 + load2*down3 + load33*depressed3
# Factor Identification: Standardize Factors at T1
## Fix Factor Means at T1 to Zero
latent_t1 ~ 0*1
## Fix Factor Variances at T1 to One
latent_t1 ~~ 1*latent_t1
# Freely Estimate Factor Means at T2 and T3 (relative to T1)
latent_t2 ~ 1
latent_t3 ~ 1
# Freely Estimate Factor Variances at T2 and T3 (relative to T1)
latent_t2 ~~ latent_t2
latent_t3 ~~ latent_t3
# Fix Some Intercepts Across Time
sad1 ~ inta*1
sad2 ~ inta*1
sad3 ~ inta*1
down1 ~ intb*1
down2 ~ intb2*1
down3 ~ intb*1
depressed1 ~ intc*1
depressed2 ~ intc2*1
depressed3 ~ intc*1
# Residual Covariances Within Indicator Across Time
sad1 ~~ sad2
sad2 ~~ sad3
sad3 ~~ sad3
down1 ~~ depressed2
down2 ~~ depressed3
down3 ~~ depressed3
depressed1 ~~ depressed2
depressed2 ~~ depressed3
depressed3 ~~ depressed3
'
partialStrongInvarianceStandardizeT1Latent_fit <- cfa(
partialStrongInvarianceStandardizeT1Latent_syntax,
data = mydata_wide,
missing = "ML",
estimator = "MLR",
meanstructure = TRUE,
#std.lv = TRUE,
fixed.x = FALSE)
summary(
partialStrongInvarianceStandardizeT1Latent_fit,
fit.measures = TRUE,
standardized = TRUE,
rsquare = TRUE)
anova(
configuralInvarianceStandardizeT1Latent_fit,
partialStrongInvarianceStandardizeT1Latent_fit
)
```
### Partial Strong Factorial Invariance With Second-Order Growth Curve {#sec-secondOrderLGCM}
#### Model Syntax
```{r}
partialStrongInvarianceGrowthCurve_syntax <- '
# Factor Loadings
latent_t1 =~ NA*load1*sad1 + load2*down1 + load3*depressed1
latent_t2 =~ NA*load1*sad2 + load2*down2 + load3*depressed2
latent_t3 =~ NA*load1*sad3 + load2*down3 + load33*depressed3
# Factor Identification: Standardize Factors at T1
## Fix Factor Means at T1 to Zero
latent_t1 ~ 0*1
## Fix Factor Variances at T1 to One
latent_t1 ~~ 1*latent_t1
# Fix Some Intercepts Across Time
sad1 ~ inta*1
sad2 ~ inta*1
sad3 ~ inta*1
down1 ~ intb*1
down2 ~ intb2*1
down3 ~ intb*1
depressed1 ~ intc*1
depressed2 ~ intc2*1
depressed3 ~ intc*1
# Residual Covariances Within Indicator Across Time
sad1 ~~ sad2
sad2 ~~ sad3
sad3 ~~ sad3
down1 ~~ depressed2
down2 ~~ depressed3
down3 ~~ depressed3
depressed1 ~~ depressed2
depressed2 ~~ depressed3
depressed3 ~~ depressed3
# Linear Growth Model
i =~ 1*latent_t1 + 1*latent_t2 + 1*latent_t3
s =~ 0*latent_t1 + 1*latent_t2 + 2*latent_t3
# Variance-covariances of intercepts and slopes
i ~~ i
s ~~ s
i ~~ s
# Means of level and slope
i ~ 1
s ~ 1
'
```
#### Fit Model
```{r}
partialStrongInvarianceGrowthCurve_fit <- cfa(
partialStrongInvarianceGrowthCurve_syntax,
data = mydata_wide,
missing = "ML",
estimator = "MLR",
meanstructure = TRUE,
#std.lv = TRUE,
fixed.x = FALSE)
```
#### Model Summary
```{r}
summary(
partialStrongInvarianceGrowthCurve_fit,
fit.measures = TRUE,
standardized = TRUE,
rsquare = TRUE)
```
Parameter estimates:
```{r}
lavaan::parameterEstimates(partialStrongInvarianceGrowthCurve_fit) %>%
subset(
lhs %in% c("i", "s") & substr(rhs, 1, 6) != "latent",
-label
)
```
# Alignment Optimization {#sec-alignment}
The example is adapted from the section on [longitudinal measurement invariance](#sec-longitudinalMI).
Two primary options for identifying a latent factor include:
1. standardizing the latent factor(s) at all timepoints (fixing their intercepts to zero and their variances to one)
1. standardizing the latent factor(s) at T1 (fixing their intercepts to zero and their variances to one) and constraining the factor loadings and intercepts of one item to be equal across time
(Option #1 may be preferable in alignment optimization.)
## Standardizing the Latent Factors at All Timepoints
### Configural Invariance With Correlated Residuals Within-Indicator Across Time
1. Standardize the latent factor(s) at all timepoints (i.e., fix the mean to zero and the variance to one)
1. Allow within-indicator residuals to covary across time
#### Model Syntax
```{r}
configuralInvarianceStandardizeAllLatent_syntax <- '
# Factor Loadings
latent_t1 =~ NA*load11*sad1 + load12*down1 + load13*depressed1
latent_t2 =~ NA*load21*sad2 + load22*down2 + load23*depressed2
latent_t3 =~ NA*load31*sad3 + load32*down3 + load33*depressed3
# Factor Identification: Standardize Factors at All Timepoints (using std.lv = TRUE)
# Free Intercepts of Manifest Variables
sad1 ~ inta1*1
sad2 ~ inta2*1
sad3 ~ inta3*1
down1 ~ intb1*1
down2 ~ intb2*1
down3 ~ intb3*1
depressed1 ~ intc1*1
depressed2 ~ intc2*1
depressed3 ~ intc3*1
# Residual Covariances Within Indicator Across Time
sad1 ~~ sad2
sad2 ~~ sad3
sad3 ~~ sad3
down1 ~~ depressed2
down2 ~~ depressed3
down3 ~~ depressed3
depressed1 ~~ depressed2
depressed2 ~~ depressed3
depressed3 ~~ depressed3
'
```
#### Fit Model
```{r}
configuralInvarianceStandardizeAllLatent_fit <- cfa(
configuralInvarianceStandardizeAllLatent_syntax,
data = mydata_wide,
missing = "ML",
estimator = "MLR",
meanstructure = TRUE,
std.lv = TRUE,
fixed.x = FALSE)
```
#### Model Summary
```{r}
summary(
configuralInvarianceStandardizeAllLatent_fit,
fit.measures = TRUE,
standardized = TRUE,
rsquare = TRUE)
```
### Longitudinal Alignment Optimization
Adapted from Lai (2013; [code](https://github.com/marklhc/awc-growth-supp/blob/main/ex_MIDUS_growth.Rmd)):
```{r}
# Extract loadings and intercepts for alignment
lam_mat <- lavInspect(configuralInvarianceStandardizeAllLatent_fit, what = "est")$lambda
nu_vec <- lavInspect(configuralInvarianceStandardizeAllLatent_fit, what = "est")$nu
# Put them into T x p matrices
num_items <- 3
num_waves <- 3
lam_config <- crossprod(lam_mat, rep(1, num_waves) %x% diag(num_items))
nu_config <- matrix(nu_vec, nrow = num_waves, ncol = num_items, byrow = TRUE)
# Add indicator names
colnames(lam_config) <- colnames(nu_config) <- c("sad", "down", "depressed")
# Alignment optimization
aligned_pars <- sirt::invariance.alignment(
lambda = lam_config,
nu = nu_config,
fixed = TRUE
)
```
The aligned loadings and intercepts are below:
```{r}
alignedParameters <- rbind(
"Loadings" = rep(NA, num_items),
`rownames<-`(aligned_pars$lambda.aligned,
paste0("Time", 1:num_waves)),
"Intercepts" = rep(NA, num_items),
`rownames<-`(aligned_pars$nu.aligned,
paste0("Time", 1:num_waves))
)
#opts <- options(knitr.kable.NA = "")
#
#knitr::kable(
# alignedParameters,
# format = "simple",
# digits = 3L
#)
print(alignedParameters, na.print = "")
aligned_pars$lambda.aligned
aligned_pars$nu.aligned
```
#### Effect Size of Noninvariance ($d_{MACS}$)
Lai (2023) suggests using a "50/30/20" rule of thumb for using $d_{MACS}$ to evaluate the appropriateness of the alignment optimization method:
Alignment optimization within CFA "is trustworthy when (a) no more than 50% of items have one or more $d_{MACS}$ > .20 and (b) no more than 30% of the pairwise $d_{MACS}$ > .20.
Lai (2023) suggests reporting the range of $d_{MACS}$ values, the proportion of $d_{MACS}$ > .20, and the proportion of items with at least one $d_{MACS}$ > .20.
```{r}
# Function for dMACS
dmacs <- function(loadings, intercepts, pooled_item_sd,
latent_mean = 0, latent_var = 1) {
dloading <- diff(loadings)
dintercept <- diff(intercepts)
integral <- dintercept^2 + 2 * dintercept * dloading * latent_mean +
dloading^2 * (latent_var + latent_mean^2)
sqrt(integral) / pooled_item_sd
}
# Use item SDs at first time point
item_sds_wave1 <-
apply(
mydata_wide[c("sad1", "down1", "depressed1")],
2, sd
)
dmacs_pairwise <- function(loading_mat, intercept_mat, pooled_item_sd,
latent_mean = 0, latent_var = 1) {
ngroups <- nrow(loading_mat)
pairs <- combn(ngroups, 2)
out <- matrix(NA, nrow = ncol(pairs), ncol = ncol(loading_mat))
for (i in seq_len(ncol(pairs))) {
out[i, ] <- dmacs(loading_mat[pairs[, i], ],
intercepts = intercept_mat[pairs[, i], ],
pooled_item_sd,
latent_mean,
latent_var
)
}
rownames(out) <- apply(pairs, 2, paste, collapse = " vs ")
colnames(out) <- colnames(loading_mat)
out
}
# All pairwise dMACS
dmacs_pairwise(aligned_pars$lambda.aligned,
intercept_mat = aligned_pars$nu.aligned,
pooled_item_sd = item_sds_wave1,
latent_mean = 0,
latent_var = 1
)
```
### An Equivalent Configural Model with Aligned Loadings
For the reference indicator (for which the loadings and intercepts are fixed), it is preferable to select an indicator with a strong factor loading (Lai, 2023).
```{r}
equivalentConfiguralModelAligned_syntax <- '
# Factor Loadings (Loadings of first indicator fixed to alignment solution)
latent_t1 =~ 0.9904750*load11*sad1 + load12*down1 + load13*depressed1
latent_t2 =~ 0.9880715*load21*sad2 + load22*down2 + load23*depressed2
latent_t3 =~ 0.9863225*load31*sad3 + load32*down3 + load33*depressed3
# Intercepts of Manifest Variables (Intercepts of first indicator fixed to alignment solution)
sad1 ~ 3.319192*inta1*1
sad2 ~ 3.332276*inta2*1
sad3 ~ 3.345631*inta3*1
down1 ~ intb1*1
down2 ~ intb2*1
down3 ~ intb3*1
depressed1 ~ intc1*1
depressed2 ~ intc2*1
depressed3 ~ intc3*1
# Residual Covariances Within Indicator Across Time
sad1 ~~ sad2
sad2 ~~ sad3
sad3 ~~ sad3
down1 ~~ depressed2
down2 ~~ depressed3
down3 ~~ depressed3
depressed1 ~~ depressed2
depressed2 ~~ depressed3
depressed3 ~~ depressed3
# Free latent means
latent_t1 + latent_t2 + latent_t3 ~ NA*1
'
equivalentConfiguralModelAligned_fit <- cfa(
equivalentConfiguralModelAligned_syntax,
data = mydata_wide,
missing = "ML",
estimator = "MLR",
meanstructure = TRUE,
#std.lv = TRUE,
fixed.x = FALSE)
```
Compare model fit (they're equivalent):
```{r}
anova(
configuralInvarianceStandardizeAllLatent_fit,
equivalentConfiguralModelAligned_fit
)
```
### Alignment-Within-CFA (AwC) Approach for Growth Modeling
#### Model Syntax
For the reference indicator (for which the loadings and intercepts are fixed), it is preferable to select an indicator with a strong factor loading (Lai, 2023).
```{r}
awcGrowthModel_syntax <- '
# Factor Loadings (Loadings of first indicator fixed to alignment solution)
latent_t1 =~ 0.9904750*load11*sad1 + load12*down1 + load13*depressed1
latent_t2 =~ 0.9880715*load21*sad2 + load22*down2 + load23*depressed2
latent_t3 =~ 0.9863225*load31*sad3 + load32*down3 + load33*depressed3
# Intercepts of Manifest Variables (Intercepts of first indicator fixed to alignment solution)
sad1 ~ 3.319192*inta1*1
sad2 ~ 3.332276*inta2*1
sad3 ~ 3.345631*inta3*1
down1 ~ intb1*1
down2 ~ intb2*1
down3 ~ intb3*1
depressed1 ~ intc1*1
depressed2 ~ intc2*1
depressed3 ~ intc3*1
# Residual Covariances Within Indicator Across Time
sad1 ~~ sad2
sad2 ~~ sad3
sad3 ~~ sad3
down1 ~~ depressed2
down2 ~~ depressed3
down3 ~~ depressed3
depressed1 ~~ depressed2
depressed2 ~~ depressed3
depressed3 ~~ depressed3
# Linear Growth Model
i =~ 1*latent_t1 + 1*latent_t2 + 1*latent_t3
s =~ 0*latent_t1 + 1*latent_t2 + 2*latent_t3
# Variance-covariances of intercepts and slopes
i ~~ i
s ~~ s
i ~~ s
# Means of level and slope
i ~ 1
s ~ 1
# Fixed disturbances of latent outcomes to zero
latent_t1 ~ 0*1
latent_t2 ~ 0*1
latent_t3 ~ 0*1
'
```
#### Fit Model
```{r}
awcGrowthModel_fit <- cfa(
awcGrowthModel_syntax,
data = mydata_wide,
missing = "ML",
estimator = "MLR",
meanstructure = TRUE,
#std.lv = TRUE,
fixed.x = FALSE)
```
#### Model Summary
```{r}
summary(
configuralInvarianceStandardizeAllLatent_fit,
fit.measures = TRUE,
standardized = TRUE,
rsquare = TRUE)
```
Parameter estimates:
```{r}
lavaan::parameterEstimates(awcGrowthModel_fit) %>%
subset(
lhs %in% c("i", "s") & substr(rhs, 1, 6) != "latent",
-label
) %>%
knitr::kable(format = "simple", digits = 3L)
```
## Standardizing the Latent Factors at T1 while Constraining Loadings and Intercepts of One Item Across Time
### Configural Invariance With Correlated Residuals Within-Indicator Across Time
1. Standardize the latent factor(s) at T1 (i.e., fix the mean to zero and the variance to one)
1. For each latent construct, constrain the first indicator's factor loading to be the same across time
1. For each latent construct, constrain the first indicator's intercept to be the same across time
1. **Allow within-indicator residuals to covary across time**
# Moderated Nonlinear Factor Analysis (MNLFA) {#sec-mnlfa}
# Item Response Theory (IRT) with Anchor Items {#sec-irt}
It can be potentially helpful to consider regularization/Bayesian penalty methods (e.g., using priors such as Horseshoe, Lasso, or Spike-and-slab).
For examples of IRT models, see [here](irt.qmd).
## Separate Calibration {#sec-irtSeparateCalibration}
## Concurrent Calibration {#sec-irtConcurrentCalibration}