This tutorial provides companion lavaan code for Chapter 14 of Growth Modeling (Grimm, Ram & Estabrook, 2017). We walk through set up a series of increasingly constrained common factor models for 3-occasion repeated measures data, compare these models to test for longitudinal factorial invariance, and then use a second order latent growth model to examine change in the invaraint factor.
The tutorial includes line-by-line code for the following models:
Load required libraries and read in the data. Data is in wide format for SEM approach.
The ECLSK dataset is a wide-form dataset that contain science, reading, and math aptitude scores, obtained in 3rd, 5th, and 8th grade. These data will be used to illustrate procedures for examining change in an “academic achievement” latent factor.
library(lavaan)
## This is lavaan 0.5-22
## lavaan is BETA software! Please report any bugs.
library(ggplot2)
#set filepath
filepath<-"https://quantdev.ssri.psu.edu/sites/qdev/files/ECLS_Science.csv"
#read in the .csv file using the url() function
eclsk <- read.csv(file=url(filepath),header=TRUE)
#looking at top of data file
head(eclsk)
## id s_g3 r_g3 m_g3 s_g5 r_g5 m_g5 s_g8 r_g8 m_g8 st_g3
## 1 1 NA NA NA NA NA NA NA NA NA NA
## 2 3 NA NA NA NA NA NA NA NA NA NA
## 3 8 NA NA NA NA NA NA 103.9 204.10 166.67 NA
## 4 16 51.57 142.18 115.59 65.94 141.02 133.67 86.9 169.83 156.67 -0.241
## 5 28 NA NA NA NA NA NA NA NA NA NA
## 6 44 NA NA NA NA NA NA NA NA NA NA
## rt_g3 mt_g3 st_g5 rt_g5 mt_g5 st_g8 rt_g8 mt_g8
## 1 NA NA NA NA NA NA NA NA
## 2 NA NA NA NA NA NA NA NA
## 3 NA NA NA NA NA 2.236 1.987 2.041
## 4 1.413 1.601 0.476 1.161 1.512 0.950 1.243 1.704
## 5 NA NA NA NA NA NA NA NA
## 6 NA NA NA NA NA NA NA NA
The configural invariance model imposes the least constraints on the factor structure across time. The only constraint is that the number of factors and pattern of zero and nonzero loadings of the common factor model are identical across measurement occasions. We define an “academic achievement” factor that is indicated by the math, science, and reading variables, seprately at each of the three measurement occasions.
While the common factor at each of the measurement occasions can be interpreted similarly (e.g., “academic achievement”), this model does not impose or assume that the factors measures the same construct across time, or that the factors are measured on the same scale.
The configural invariance model is specified as …
config.invar<-'#factor loadings
eta1 =~ lambda_S*s_g3+
lambda_R3*r_g3+
lambda_M3*m_g3
eta2 =~ lambda_S*s_g5+
lambda_R5*r_g5+
lambda_M5*m_g5
eta3 =~ lambda_S*s_g8+
lambda_R8*r_g8+
lambda_M8*m_g8
#latent variable variances
eta1~~1*eta1
eta2~~eta2
eta3~~eta3
#latent variable covariances
eta1~~eta2
eta1~~eta3
eta2~~eta3
#unique variances
s_g3~~s_g3
s_g5~~s_g5
s_g8~~s_g8
r_g3~~r_g3
r_g5~~r_g5
r_g8~~r_g8
m_g3~~m_g3
m_g5~~m_g5
m_g8~~m_g8
#unique covariances
s_g3~~s_g5
s_g3~~s_g8
s_g5~~s_g8
r_g3~~r_g5
r_g3~~r_g8
r_g5~~r_g8
m_g3~~m_g5
m_g3~~m_g8
m_g5~~m_g8
#latent variable intercepts
eta1~0*1
eta2~1
eta3~1
#observed variable intercepts
s_g3~tau_S*1
s_g5~tau_S*1
s_g8~tau_S*1
r_g3~tau_R3*1
r_g5~tau_R5*1
r_g8~tau_R8*1
m_g3~tau_M3*1
m_g5~tau_M5*1
m_g8~tau_M8*1'
lavaan.fit1<-lavaan(config.invar,
data = eclsk,
meanstructure = TRUE,
estimator = "ML",
missing = "fiml",
fixed.x = FALSE)
summary(lavaan.fit1, fit.measures=TRUE)
## lavaan (0.5-22) converged normally after 298 iterations
##
## Used Total
## Number of observations 1478 2108
##
## Number of missing patterns 24
##
## Estimator ML
## Minimum Function Test Statistic 35.522
## Degrees of freedom 15
## P-value (Chi-square) 0.002
##
## Model test baseline model:
##
## Minimum Function Test Statistic 11669.413
## Degrees of freedom 36
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.998
## Tucker-Lewis Index (TLI) 0.996
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -41918.095
## Loglikelihood unrestricted model (H1) -41900.334
##
## Number of free parameters 39
## Akaike (AIC) 83914.190
## Bayesian (BIC) 84120.830
## Sample-size adjusted Bayesian (BIC) 83996.938
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.030
## 90 Percent Confidence Interval 0.018 0.043
## P-value RMSEA <= 0.05 0.994
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.007
##
## Parameter Estimates:
##
## Information Observed
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## eta1 =~
## s_g3 (lm_S) 13.269 0.339 39.163 0.000
## r_g3 (l_R3) 26.301 0.627 41.921 0.000
## m_g3 (l_M3) 21.601 0.553 39.074 0.000
## eta2 =~
## s_g5 (lm_S) 13.269 0.339 39.163 0.000
## r_g5 (l_R5) 22.745 0.698 32.598 0.000
## m_g5 (l_M5) 20.066 0.631 31.789 0.000
## eta3 =~
## s_g8 (lm_S) 13.269 0.339 39.163 0.000
## r_g8 (l_R8) 21.663 0.753 28.763 0.000
## m_g8 (l_M8) 17.259 0.593 29.129 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## eta1 ~~
## eta2 1.034 0.021 49.169 0.000
## eta3 1.050 0.029 36.246 0.000
## eta2 ~~
## eta3 1.176 0.048 24.703 0.000
## .s_g3 ~~
## .s_g5 34.822 3.060 11.379 0.000
## .s_g8 14.743 3.006 4.904 0.000
## .s_g5 ~~
## .s_g8 18.619 3.177 5.861 0.000
## .r_g3 ~~
## .r_g5 82.590 9.155 9.021 0.000
## .r_g8 54.235 9.416 5.760 0.000
## .r_g5 ~~
## .r_g8 60.822 9.129 6.663 0.000
## .m_g3 ~~
## .m_g5 123.727 8.135 15.208 0.000
## .m_g8 82.025 6.914 11.864 0.000
## .m_g5 ~~
## .m_g8 91.527 7.139 12.821 0.000
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## eta1 0.000
## eta2 1.040 0.033 31.645 0.000
## eta3 2.435 0.068 35.961 0.000
## .s_g3 (ta_S) 51.013 0.407 125.219 0.000
## .s_g5 (ta_S) 51.013 0.407 125.219 0.000
## .s_g8 (ta_S) 51.013 0.407 125.219 0.000
## .r_g3 (t_R3) 127.260 0.772 164.876 0.000
## .r_g5 (t_R5) 126.654 0.997 127.087 0.000
## .r_g8 (t_R8) 116.432 1.678 69.373 0.000
## .m_g3 (t_M3) 99.744 0.668 149.395 0.000
## .m_g5 (t_M5) 102.902 0.908 113.381 0.000
## .m_g8 (t_M8) 98.540 1.315 74.929 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## eta1 1.000
## eta2 1.160 0.046 25.220 0.000
## eta3 1.329 0.069 19.227 0.000
## .s_g3 67.122 3.535 18.989 0.000
## .s_g5 63.195 3.842 16.447 0.000
## .s_g8 53.716 4.113 13.059 0.000
## .r_g3 180.361 11.527 15.647 0.000
## .r_g5 162.710 10.586 15.371 0.000
## .r_g8 185.198 12.274 15.089 0.000
## .m_g3 187.384 9.448 19.833 0.000
## .m_g5 176.465 9.534 18.510 0.000
## .m_g8 126.219 7.760 16.265 0.000
The weak (or metric) invariance model imposes the additional constraint that the factor loading matrix (lambda) is identical at all measurement occasions. The equality constraints on the factor laodings impose that the covariance structures across time are proportional. However, the weak invariance model allows that the observed variable intercepts may vary across time and thus does not provide for examination of longitudinal change in the factor.
The weak invariance model is specified as …
weak.invar<-'#factor loadings (with constraints)
eta1 =~ lambda_S*s_g3+
lambda_R*r_g3+
lambda_M*m_g3
eta2 =~ lambda_S*s_g5+
lambda_R*r_g5+
lambda_M*m_g5
eta3 =~ lambda_S*s_g8+
lambda_R*r_g8+
lambda_M*m_g8
#latent variable variances
eta1~~1*eta1
eta2~~eta2
eta3~~eta3
#latent variable covariances
eta1~~eta2
eta1~~eta3
eta2~~eta3
#unique variances
s_g3~~s_g3
s_g5~~s_g5
s_g8~~s_g8
r_g3~~r_g3
r_g5~~r_g5
r_g8~~r_g8
m_g3~~m_g3
m_g5~~m_g5
m_g8~~m_g8
#unique covariances
s_g3~~s_g5
s_g3~~s_g8
s_g5~~s_g8
r_g3~~r_g5
r_g3~~r_g8
r_g5~~r_g8
m_g3~~m_g5
m_g3~~m_g8
m_g5~~m_g8
#latent variable intercepts
eta1~0*1
eta2~1
eta3~1
#observed variable intercepts
s_g3~tau_S*1
s_g5~tau_S*1
s_g8~tau_S*1
r_g3~tau_R3*1
r_g5~tau_R5*1
r_g8~tau_R8*1
m_g3~tau_M3*1
m_g5~tau_M5*1
m_g8~tau_M8*1'
lavaan.fit2<-lavaan(weak.invar,
data = eclsk,
meanstructure = TRUE,
estimator = "ML",
missing = "fiml",
fixed.x = FALSE)
summary(lavaan.fit2, fit.measures=TRUE)
## lavaan (0.5-22) converged normally after 240 iterations
##
## Used Total
## Number of observations 1478 2108
##
## Number of missing patterns 24
##
## Estimator ML
## Minimum Function Test Statistic 116.826
## Degrees of freedom 19
## P-value (Chi-square) 0.000
##
## Model test baseline model:
##
## Minimum Function Test Statistic 11669.413
## Degrees of freedom 36
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.992
## Tucker-Lewis Index (TLI) 0.984
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -41958.747
## Loglikelihood unrestricted model (H1) -41900.334
##
## Number of free parameters 35
## Akaike (AIC) 83987.495
## Bayesian (BIC) 84172.940
## Sample-size adjusted Bayesian (BIC) 84061.756
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.059
## 90 Percent Confidence Interval 0.049 0.070
## P-value RMSEA <= 0.05 0.069
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.042
##
## Parameter Estimates:
##
## Information Observed
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## eta1 =~
## s_g3 (lm_S) 14.148 0.328 43.128 0.000
## r_g3 (lm_R) 25.418 0.590 43.112 0.000
## m_g3 (lm_M) 20.876 0.524 39.821 0.000
## eta2 =~
## s_g5 (lm_S) 14.148 0.328 43.128 0.000
## r_g5 (lm_R) 25.418 0.590 43.112 0.000
## m_g5 (lm_M) 20.876 0.524 39.821 0.000
## eta3 =~
## s_g8 (lm_S) 14.148 0.328 43.128 0.000
## r_g8 (lm_R) 25.418 0.590 43.112 0.000
## m_g8 (lm_M) 20.876 0.524 39.821 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## eta1 ~~
## eta2 0.961 0.014 71.113 0.000
## eta3 0.902 0.019 47.928 0.000
## eta2 ~~
## eta3 0.936 0.027 34.664 0.000
## .s_g3 ~~
## .s_g5 33.360 3.070 10.866 0.000
## .s_g8 14.355 3.054 4.700 0.000
## .s_g5 ~~
## .s_g8 22.119 3.248 6.811 0.000
## .r_g3 ~~
## .r_g5 78.937 9.183 8.596 0.000
## .r_g8 52.347 9.433 5.549 0.000
## .r_g5 ~~
## .r_g8 54.402 9.102 5.977 0.000
## .m_g3 ~~
## .m_g5 129.396 8.261 15.664 0.000
## .m_g8 82.602 7.027 11.754 0.000
## .m_g5 ~~
## .m_g8 92.414 7.292 12.673 0.000
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## eta1 0.000
## eta2 0.977 0.029 33.565 0.000
## eta3 2.293 0.059 38.788 0.000
## .s_g3 (ta_S) 51.012 0.424 120.209 0.000
## .s_g5 (ta_S) 51.012 0.424 120.209 0.000
## .s_g8 (ta_S) 51.012 0.424 120.209 0.000
## .r_g3 (t_R3) 127.284 0.755 168.686 0.000
## .r_g5 (t_R5) 125.448 0.985 127.310 0.000
## .r_g8 (t_R8) 110.884 1.494 74.231 0.000
## .m_g3 (t_M3) 99.743 0.656 152.160 0.000
## .m_g5 (t_M5) 103.390 0.860 120.251 0.000
## .m_g8 (t_M8) 92.568 1.288 71.896 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## eta1 1.000
## eta2 1.000 0.026 37.899 0.000
## eta3 0.979 0.036 27.053 0.000
## .s_g3 63.797 3.533 18.060 0.000
## .s_g5 64.271 3.837 16.751 0.000
## .s_g8 62.052 4.178 14.851 0.000
## .r_g3 187.090 11.494 16.277 0.000
## .r_g5 153.835 10.547 14.585 0.000
## .r_g8 179.482 12.171 14.746 0.000
## .m_g3 194.484 9.513 20.443 0.000
## .m_g5 183.051 9.672 18.926 0.000
## .m_g8 122.637 7.860 15.603 0.000
The strong invariance model constrains observed variable intercepts (tau) to be equal across time, while means of the latent variables are allowed to vary across occasions. Since all mean-level change in the observed variables can be attributed to the factors, the scale of the latent variables can be assumed to be equall across measurement occasions.
The strong invariance model is specified as …
strong.invar<-'#factor loadings (with constraints)
eta1 =~ lambda_S*s_g3+
lambda_R*r_g3+
lambda_M*m_g3
eta2 =~ lambda_S*s_g5+
lambda_R*r_g5+
lambda_M*m_g5
eta3 =~ lambda_S*s_g8+
lambda_R*r_g8+
lambda_M*m_g8
#latent variable variances
eta1~~1*eta1
eta2~~eta2
eta3~~eta3
#latent variable covariances
eta1~~eta2
eta1~~eta3
eta2~~eta3
#unique variances
s_g3~~s_g3
s_g5~~s_g5
s_g8~~s_g8
r_g3~~r_g3
r_g5~~r_g5
r_g8~~r_g8
m_g3~~m_g3
m_g5~~m_g5
m_g8~~m_g8
#unique covariances
s_g3~~s_g5
s_g3~~s_g8
s_g5~~s_g8
r_g3~~r_g5
r_g3~~r_g8
r_g5~~r_g8
m_g3~~m_g5
m_g3~~m_g8
m_g5~~m_g8
#latent variable intercepts
eta1~0*1
eta2~1
eta3~1
#observed variable intercepts (with constraints)
s_g3~tau_S*1
s_g5~tau_S*1
s_g8~tau_S*1
r_g3~tau_R*1
r_g5~tau_R*1
r_g8~tau_R*1
m_g3~tau_M*1
m_g5~tau_M*1
m_g8~tau_M*1'
lavaan.fit3<-lavaan(strong.invar,
data = eclsk,
meanstructure = TRUE,
estimator = "ML",
missing = "fiml",
fixed.x = FALSE)
summary(lavaan.fit3, fit.measures=TRUE)
## lavaan (0.5-22) converged normally after 226 iterations
##
## Used Total
## Number of observations 1478 2108
##
## Number of missing patterns 24
##
## Estimator ML
## Minimum Function Test Statistic 540.750
## Degrees of freedom 23
## P-value (Chi-square) 0.000
##
## Model test baseline model:
##
## Minimum Function Test Statistic 11669.413
## Degrees of freedom 36
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.955
## Tucker-Lewis Index (TLI) 0.930
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -42170.709
## Loglikelihood unrestricted model (H1) -41900.334
##
## Number of free parameters 31
## Akaike (AIC) 84403.419
## Bayesian (BIC) 84567.670
## Sample-size adjusted Bayesian (BIC) 84469.193
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.123
## 90 Percent Confidence Interval 0.115 0.133
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.075
##
## Parameter Estimates:
##
## Information Observed
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## eta1 =~
## s_g3 (lm_S) 15.088 0.322 46.858 0.000
## r_g3 (lm_R) 22.740 0.524 43.414 0.000
## m_g3 (lm_M) 20.807 0.473 44.030 0.000
## eta2 =~
## s_g5 (lm_S) 15.088 0.322 46.858 0.000
## r_g5 (lm_R) 22.740 0.524 43.414 0.000
## m_g5 (lm_M) 20.807 0.473 44.030 0.000
## eta3 =~
## s_g8 (lm_S) 15.088 0.322 46.858 0.000
## r_g8 (lm_R) 22.740 0.524 43.414 0.000
## m_g8 (lm_M) 20.807 0.473 44.030 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## eta1 ~~
## eta2 0.976 0.014 70.707 0.000
## eta3 0.910 0.019 47.377 0.000
## eta2 ~~
## eta3 0.961 0.028 34.465 0.000
## .s_g3 ~~
## .s_g5 29.487 3.487 8.456 0.000
## .s_g8 7.394 3.251 2.274 0.023
## .s_g5 ~~
## .s_g8 7.867 3.332 2.361 0.018
## .r_g3 ~~
## .r_g5 109.044 9.413 11.585 0.000
## .r_g8 73.172 9.982 7.330 0.000
## .r_g5 ~~
## .r_g8 71.480 9.209 7.762 0.000
## .m_g3 ~~
## .m_g5 128.408 8.511 15.088 0.000
## .m_g8 85.517 7.421 11.523 0.000
## .m_g5 ~~
## .m_g8 90.389 7.594 11.903 0.000
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## eta1 0.000
## eta2 1.032 0.026 40.382 0.000
## eta3 1.972 0.046 42.668 0.000
## .s_g3 (ta_S) 51.422 0.450 114.287 0.000
## .s_g5 (ta_S) 51.422 0.450 114.287 0.000
## .s_g8 (ta_S) 51.422 0.450 114.287 0.000
## .r_g3 (ta_R) 126.310 0.706 179.023 0.000
## .r_g5 (ta_R) 126.310 0.706 179.023 0.000
## .r_g8 (ta_R) 126.310 0.706 179.023 0.000
## .m_g3 (ta_M) 100.035 0.658 152.097 0.000
## .m_g5 (ta_M) 100.035 0.658 152.097 0.000
## .m_g8 (ta_M) 100.035 0.658 152.097 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## eta1 1.000
## eta2 1.025 0.027 37.549 0.000
## eta3 0.990 0.037 26.783 0.000
## .s_g3 58.999 3.875 15.224 0.000
## .s_g5 63.402 4.314 14.698 0.000
## .s_g8 58.256 4.688 12.427 0.000
## .r_g3 230.033 11.800 19.494 0.000
## .r_g5 181.922 10.581 17.193 0.000
## .r_g8 208.269 12.238 17.019 0.000
## .m_g3 197.179 9.868 19.981 0.000
## .m_g5 188.328 10.123 18.603 0.000
## .m_g8 126.741 8.199 15.458 0.000
Finally, the strict invariance model imposes additional constraints on the longitudinal factor model. Here, factor loadings, observed variable intercepts, and unique variances are all constrained to be equal across the repeated measures. When strict factorial invariance is in place, longitudinal changes in observed means, variances, and covariances of the factors can be examined for developmental change.
Following the path diagram shown in Figure 14.2 (p. 349), the strict invariance model is specified as ….
strict.invar<-'#common factor loadings (with constraints)
eta1 =~ lambda_S*s_g3+
lambda_R*r_g3+
lambda_M*m_g3
eta2 =~ lambda_S*s_g5+
lambda_R*r_g5+
lambda_M*m_g5
eta3 =~ lambda_S*s_g8+
lambda_R*r_g8+
lambda_M*m_g8
#latent variable variances
eta1~~1*eta1
eta2~~eta2
eta3~~eta3
#latent variable covariances
eta1~~eta2
eta1~~eta3
eta2~~eta3
#unique variances (with constraints)
s_g3~~theta_S*s_g3
s_g5~~theta_S*s_g5
s_g8~~theta_S*s_g8
r_g3~~theta_R*r_g3
r_g5~~theta_R*r_g5
r_g8~~theta_R*r_g8
m_g3~~theta_M*m_g3
m_g5~~theta_M*m_g5
m_g8~~theta_M*m_g8
#unique covariances
s_g3~~s_g5
s_g3~~s_g8
s_g5~~s_g8
r_g3~~r_g5
r_g3~~r_g8
r_g5~~r_g8
m_g3~~m_g5
m_g3~~m_g8
m_g5~~m_g8
#latent variable intercepts
eta1~0*1
eta2~1
eta3~1
#observed variable intercepts (with constraints)
s_g3~tau_S*1
s_g5~tau_S*1
s_g8~tau_S*1
r_g3~tau_R*1
r_g5~tau_R*1
r_g8~tau_R*1
m_g3~tau_M*1
m_g5~tau_M*1
m_g8~tau_M*1'
lavaan.fit4<-lavaan(strict.invar,
data = eclsk,
meanstructure = TRUE,
estimator = "ML",
missing = "fiml",
fixed.x = FALSE)
summary(lavaan.fit4, fit.measures=TRUE)
## lavaan (0.5-22) converged normally after 260 iterations
##
## Used Total
## Number of observations 1478 2108
##
## Number of missing patterns 24
##
## Estimator ML
## Minimum Function Test Statistic 600.614
## Degrees of freedom 29
## P-value (Chi-square) 0.000
##
## Model test baseline model:
##
## Minimum Function Test Statistic 11669.413
## Degrees of freedom 36
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.951
## Tucker-Lewis Index (TLI) 0.939
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -42200.641
## Loglikelihood unrestricted model (H1) -41900.334
##
## Number of free parameters 25
## Akaike (AIC) 84451.282
## Bayesian (BIC) 84583.743
## Sample-size adjusted Bayesian (BIC) 84504.325
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.115
## 90 Percent Confidence Interval 0.108 0.124
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.078
##
## Parameter Estimates:
##
## Information Observed
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## eta1 =~
## s_g3 (lm_S) 15.176 0.319 47.527 0.000
## r_g3 (lm_R) 22.982 0.520 44.172 0.000
## m_g3 (lm_M) 21.236 0.473 44.918 0.000
## eta2 =~
## s_g5 (lm_S) 15.176 0.319 47.527 0.000
## r_g5 (lm_R) 22.982 0.520 44.172 0.000
## m_g5 (lm_M) 21.236 0.473 44.918 0.000
## eta3 =~
## s_g8 (lm_S) 15.176 0.319 47.527 0.000
## r_g8 (lm_R) 22.982 0.520 44.172 0.000
## m_g8 (lm_M) 21.236 0.473 44.918 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## eta1 ~~
## eta2 0.968 0.013 72.280 0.000
## eta3 0.904 0.019 48.240 0.000
## eta2 ~~
## eta3 0.949 0.027 35.094 0.000
## .s_g3 ~~
## .s_g5 28.853 3.118 9.254 0.000
## .s_g8 6.616 3.301 2.004 0.045
## .s_g5 ~~
## .s_g8 5.354 3.311 1.617 0.106
## .r_g3 ~~
## .r_g5 112.069 8.673 12.921 0.000
## .r_g8 69.529 9.345 7.440 0.000
## .r_g5 ~~
## .r_g8 79.188 9.640 8.214 0.000
## .m_g3 ~~
## .m_g5 112.589 7.308 15.407 0.000
## .m_g8 93.467 7.990 11.698 0.000
## .m_g5 ~~
## .m_g8 102.538 8.014 12.794 0.000
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## eta1 0.000
## eta2 1.020 0.025 41.017 0.000
## eta3 1.947 0.045 43.249 0.000
## .s_g3 (ta_S) 51.433 0.448 114.745 0.000
## .s_g5 (ta_S) 51.433 0.448 114.745 0.000
## .s_g8 (ta_S) 51.433 0.448 114.745 0.000
## .r_g3 (ta_R) 126.363 0.702 179.899 0.000
## .r_g5 (ta_R) 126.363 0.702 179.899 0.000
## .r_g8 (ta_R) 126.363 0.702 179.899 0.000
## .m_g3 (ta_M) 100.193 0.651 153.913 0.000
## .m_g5 (ta_M) 100.193 0.651 153.913 0.000
## .m_g8 (ta_M) 100.193 0.651 153.913 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## eta1 1.000
## eta2 1.011 0.026 38.332 0.000
## eta3 0.971 0.036 27.168 0.000
## .s_g3 (th_S) 60.274 2.939 20.506 0.000
## .s_g5 (th_S) 60.274 2.939 20.506 0.000
## .s_g8 (th_S) 60.274 2.939 20.506 0.000
## .r_g3 (th_R) 208.701 8.172 25.537 0.000
## .r_g5 (th_R) 208.701 8.172 25.537 0.000
## .r_g8 (th_R) 208.701 8.172 25.537 0.000
## .m_g3 (th_M) 173.321 7.236 23.952 0.000
## .m_g5 (th_M) 173.321 7.236 23.952 0.000
## .m_g8 (th_M) 173.321 7.236 23.952 0.000
These results match those given in Outputs 14.1 and 14.2 (pp. 359-360).
Fit statistics for the 4 invariance models are collected in the table below.
Model | Chi-sq (df) | AIC | BIC | RMSEA | CFI | TLI |
---|---|---|---|---|---|---|
Configural Invariance | 35.52 (15) | 83914.19 | 84120.83 | 0.030 | 0.998 | 0.996 |
Weak Invariance | 116.83 (19) | 83987.50 | 84172.94 | 0.059 | 0.992 | 0.984 |
Strong Invariance | 540.75 (23) | 84403.42 | 84567.67 | 0.123 | 0.955 | 0.930 |
Strict Invariance | 600.61 (29) | 84451.28 | 84583.74 | 0.115 | 0.951 | 0.939 |
Statsitical improvement/decrements in model fit that accompany each set of constraints can be examined using chi-square difference tests. With these data, all three pairs of tests show significant difference in chi-square, indicating that the more constrained models fit worse than the less-constrained models.
#Comparing configural and weak invariance
anova(lavaan.fit1,lavaan.fit2)
## Chi Square Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## lavaan.fit1 15 83914 84121 35.522
## lavaan.fit2 19 83987 84173 116.826 81.304 4 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Comparing weak and strong invariance
anova(lavaan.fit2,lavaan.fit3)
## Chi Square Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## lavaan.fit2 19 83987 84173 116.83
## lavaan.fit3 23 84403 84568 540.75 423.92 4 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Comparing strong and strict invariance
anova(lavaan.fit3,lavaan.fit4)
## Chi Square Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## lavaan.fit3 23 84403 84568 540.75
## lavaan.fit4 29 84451 84584 600.61 59.863 6 4.798e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Relatively worse fit of the more constrained models is also reflected in the AIC and BIC. However, the the CFI and TLI suggest that the strong and strict invariance models do have adequate model fit. Therefore, for the purposes of this example, where there is incentive to be able to examine change across time, the strict invariance model is accepted - and we move on to construct a second-order growth model for change in the invariant “academic achievement” factor.
Once strong or strict invariance is supported for the longitudinal factor model, changes in the latent factor scores can be examined using a second order growth model.
In this example, measurement invariance is imposed using a strict invariance longitudinal factor model, with growth in children’s latent academic achievement then modeled via a second-order latent basis growth model.
Following the path diagram shown in Figure 14.3 (p. 352), and fixing one factor loading to obtain an interpretable metric (see text on p. 351) the second-order growth model is specified as …
LG.strict.invar<-'#factor loadings (with constraints)
eta1 =~ 15.1749088*s_g3+
lambda_R*r_g3+
lambda_M*m_g3
eta2 =~ 15.1749088*s_g5+
lambda_R*r_g5+
lambda_M*m_g5
eta3 =~ 15.1749088*s_g8+
lambda_R*r_g8+
lambda_M*m_g8
#factor variances
eta1~~psi*eta1
eta2~~psi*eta2
eta3~~psi*eta3
#unique variances (with constraints)
s_g3~~theta_S*s_g3
s_g5~~theta_S*s_g5
s_g8~~theta_S*s_g8
r_g3~~theta_R*r_g3
r_g5~~theta_R*r_g5
r_g8~~theta_R*r_g8
m_g3~~theta_M*m_g3
m_g5~~theta_M*m_g5
m_g8~~theta_M*m_g8
#unique covariances
s_g3~~s_g5
s_g3~~s_g8
s_g5~~s_g8
r_g3~~r_g5
r_g3~~r_g8
r_g5~~r_g8
m_g3~~m_g5
m_g3~~m_g8
m_g5~~m_g8
#observed variable intercepts (with constraints)
s_g3~tau_S*1
s_g5~tau_S*1
s_g8~tau_S*1
r_g3~tau_R*1
r_g5~tau_R*1
r_g8~tau_R*1
m_g3~tau_M*1
m_g5~tau_M*1
m_g8~tau_M*1
#latent basis model
xi_1 =~ 1*eta1+1*eta2+1*eta3
xi_2 =~ 0*eta1+start(0.5)*eta2+1*eta3
xi_1~~start(.8)*xi_1
xi_2~~start(.5)*xi_2
xi_1~~start(0)*xi_2
xi_1~0*1
xi_2~1 '
lavaan.fit5<-lavaan(LG.strict.invar,
data = eclsk,
meanstructure = TRUE,
estimator = "ML",
missing = "fiml",
fixed.x = FALSE)
summary(lavaan.fit5, fit.measures=TRUE)
## lavaan (0.5-22) converged normally after 189 iterations
##
## Used Total
## Number of observations 1478 2108
##
## Number of missing patterns 24
##
## Estimator ML
## Minimum Function Test Statistic 606.985
## Degrees of freedom 31
## P-value (Chi-square) 0.000
##
## Model test baseline model:
##
## Minimum Function Test Statistic 11669.413
## Degrees of freedom 36
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.950
## Tucker-Lewis Index (TLI) 0.943
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -42203.827
## Loglikelihood unrestricted model (H1) -41900.334
##
## Number of free parameters 23
## Akaike (AIC) 84453.654
## Bayesian (BIC) 84575.518
## Sample-size adjusted Bayesian (BIC) 84502.454
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.112
## 90 Percent Confidence Interval 0.104 0.120
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.078
##
## Parameter Estimates:
##
## Information Observed
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## eta1 =~
## s_g3 15.175
## r_g3 (lm_R) 22.989 0.283 81.092 0.000
## m_g3 (lm_M) 21.235 0.248 85.618 0.000
## eta2 =~
## s_g5 15.175
## r_g5 (lm_R) 22.989 0.283 81.092 0.000
## m_g5 (lm_M) 21.235 0.248 85.618 0.000
## eta3 =~
## s_g8 15.175
## r_g8 (lm_R) 22.989 0.283 81.092 0.000
## m_g8 (lm_M) 21.235 0.248 85.618 0.000
## xi_1 =~
## eta1 1.000
## eta2 1.000
## eta3 1.000
## xi_2 =~
## eta1 0.000
## eta2 0.524 0.006 85.608 0.000
## eta3 1.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## .s_g3 ~~
## .s_g5 28.704 3.098 9.265 0.000
## .s_g8 6.807 3.311 2.056 0.040
## .s_g5 ~~
## .s_g8 5.742 3.289 1.746 0.081
## .r_g3 ~~
## .r_g5 113.232 8.622 13.134 0.000
## .r_g8 67.856 9.358 7.251 0.000
## .r_g5 ~~
## .r_g8 78.009 9.529 8.187 0.000
## .m_g3 ~~
## .m_g5 113.479 7.277 15.594 0.000
## .m_g8 92.852 8.000 11.607 0.000
## .m_g5 ~~
## .m_g8 100.693 7.932 12.694 0.000
## xi_1 ~~
## xi_2 -0.061 0.019 -3.114 0.002
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## .s_g3 (ta_S) 51.423 0.449 114.450 0.000
## .s_g5 (ta_S) 51.423 0.449 114.450 0.000
## .s_g8 (ta_S) 51.423 0.449 114.450 0.000
## .r_g3 (ta_R) 126.339 0.704 179.382 0.000
## .r_g5 (ta_R) 126.339 0.704 179.382 0.000
## .r_g8 (ta_R) 126.339 0.704 179.382 0.000
## .m_g3 (ta_M) 100.181 0.653 153.454 0.000
## .m_g5 (ta_M) 100.181 0.653 153.454 0.000
## .m_g8 (ta_M) 100.181 0.653 153.454 0.000
## xi_1 0.000
## xi_2 1.947 0.024 82.219 0.000
## eta1 0.000
## eta2 0.000
## eta3 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## eta1 (psi) 0.024 0.004 6.026 0.000
## eta2 (psi) 0.024 0.004 6.026 0.000
## eta3 (psi) 0.024 0.004 6.026 0.000
## .s_g3 (th_S) 60.273 2.929 20.580 0.000
## .s_g5 (th_S) 60.273 2.929 20.580 0.000
## .s_g8 (th_S) 60.273 2.929 20.580 0.000
## .r_g3 (th_R) 208.603 8.162 25.557 0.000
## .r_g5 (th_R) 208.603 8.162 25.557 0.000
## .r_g8 (th_R) 208.603 8.162 25.557 0.000
## .m_g3 (th_M) 173.812 7.217 24.083 0.000
## .m_g5 (th_M) 173.812 7.217 24.083 0.000
## .m_g8 (th_M) 173.812 7.217 24.083 0.000
## xi_1 0.983 0.042 23.184 0.000
## xi_2 0.110 0.018 6.277 0.000
These results match those given in Outputs 14.3 and 14.4.
We can visualize the predicted common factor trajectories by extracing the predicted factor estimates of the second order growth model, using these estimates to obtain growth curves of the estimated “educational achievement” score for each individual, and plotting them using ggplot.
#extract predicted factor scores
pred<-as.data.frame(lavPredict(lavaan.fit5))
#define function to fit factor trajectories
fit.line <- function(xi1,xi2,t) {
xi1 + xi2*t
}
time<-c(0,.524,1)
#fit trajectories using estimated factor scores
out=NULL; temp=NULL
for (i in 1:nrow(pred)){
temp<-fit.line(pred$xi_1[i],pred$xi_2[i],time)
out<-rbind(out,temp)
}
out<-as.data.frame(out)
colnames(out)<-c("pred3","pred5","pred8")
out$ID<-eclsk$id
#reshape estimated scores to long format (for plotting)
outlong<-reshape(out,
varying=c("pred3","pred5","pred8"),
timevar="grade",
idvar="ID",
direction="long",sep="")
outlong<-outlong[order(outlong$ID,outlong$grade),]
#fit average trajectory (for plotting)
avg<-as.data.frame(fit.line(0, 1.947,time))
avg$grade<-c(3,5,8)
avg$ID<-1
colnames(avg)<-c("fit","grade","ID")
#plot predicted common factor trajectories (average in red)
ggplot(data=outlong,aes(x=grade,y=pred,group=ID))+
ggtitle("Predicted Trajectories: Second Order Growth Model") +
xlab("Grade")+
ylab("Predicted 'Educational Achievement' Factor Score")+
geom_line() +
geom_line(data=avg,aes(x=grade,y=fit,group=ID),color="red",size=2)