Overview

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.

Outline

The tutorial includes line-by-line code for the following models:

  1. Configurial Invariance model
  2. Weak Invariance model
  3. Strong Invariance model
  4. Strict Invariance model
  5. Second-order growth model of the invariant common factor

Preliminaries: Libraries and Data

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

Common factor model

Configural invariance

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

Weak invariance

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

Strong invariance

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

Strict invariance

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 summaries

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.

Second-order growth model

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.

Latent growth with strict invariance 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)