Item Factor Analysis, Measurement Invariance + 2nd order Growth Model (ECLS-K)
Nilam Ram, Kevin Grimm et al.
1 Overview
In this tutorial we walk through example of testing measurement invariance and modeling change with latent variables measured by Ordinal indicators.
The example follows Chapter 15 of Growth Modeling: Structural Equation and Multilevel Modeling Approaches (Grimm, Ram & Estabrook, 2017).
Using 3-occasion data from the ECLS-K, we set up an item factor analysis model, test for factorial invariance, and then use a second-order growth model to describe change in the factor scores across time.
1.0.0.1 Prelim - Loading libraries used in this script.
library(psych)
library(ggplot2)
library(corrplot) #plotting correlation matrices
library(lavaan) #for fitting structural equation models
library(semPlot) #for automatically making diagrams 1.0.0.2 Prelim - Reading in Repeated Measures Data
For this example, we use an ECLS-K dataset that is in wideform. There are variables for parent’s and teachers’ ratings of children’s of children’s behavior in the fall and spring of kindergarten and spring of first grade on interpersonal skills, self-control, and externalizing behaviors. Ratings were made using a four-point (1-4) scale where higher values indicated more of the behavior.
These behavior ratings are thought to be indicators of a higher-order child behavior factor
We only use the teacher ratings here.
#set filepath for data file
filepath <- "https://raw.githubusercontent.com/LRI-2/Data/main/GrowthModeling/ECLS_Behavior.dat"
#read in the text data file using the url() function
dat <- read.table(file=url(filepath),na.strings = ".")
names(dat) <- c("id",
"p1_sc", "p1_soc", "p1_sad", "p1_imp",
"p2_sc", "p2_soc", "p2_sad", "p2_imp",
"p4_sc", "p4_soc", "p4_sad", "p4_imp",
"t1_sc", "t1_intp", "t1_ext", "t1_int",
"t2_sc", "t2_intp", "t2_ext", "t2_int",
"t4_sc", "t4_intp", "t4_ext", "t4_int")
#selecting only the teacher reports
dat <- dat[ ,c("id", "t1_sc", "t1_intp", "t1_ext", "t1_int",
"t2_sc", "t2_intp", "t2_ext", "t2_int",
"t4_sc", "t4_intp", "t4_ext", "t4_int")]1.0.0.3 Prelim - Descriptives
Lets have a quick look at the data file and the descriptives.
#look at data structure
head(dat,10)| id | t1_sc | t1_intp | t1_ext | t1_int | t2_sc | t2_intp | t2_ext | t2_int | t4_sc | t4_intp | t4_ext | t4_int |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 4 | 4 | 1 | 2 | NA | NA | NA | NA | 4 | 3 | 1 | 1 |
| 2 | 4 | 4 | 1 | 1 | 4 | 4 | 1 | 1 | 4 | 4 | 1 | 1 |
| 16 | 3 | 2 | 2 | 1 | 4 | 4 | 1 | 1 | 3 | 4 | 1 | 2 |
| 19 | 4 | 4 | 1 | 1 | 4 | 4 | 1 | 1 | 3 | 3 | 1 | 3 |
| 35 | 4 | 3 | 1 | 2 | 3 | 3 | 1 | 3 | 3 | 2 | 1 | 2 |
| 36 | 3 | 3 | 2 | 1 | 3 | 3 | 2 | 2 | 3 | 3 | 2 | 2 |
| 48 | 4 | 3 | 1 | 2 | 4 | 4 | 1 | 2 | 4 | 3 | 1 | 1 |
| 50 | 4 | 2 | 1 | 2 | 2 | 3 | 2 | 3 | 4 | 3 | 1 | 2 |
| 54 | 3 | 3 | 2 | 3 | 3 | 2 | 2 | 3 | 3 | 2 | 3 | 2 |
| 58 | 3 | NA | 3 | 2 | 3 | 2 | 4 | 3 | 3 | 3 | 3 | 2 |
#descriptives (means, sds)
describe(dat[,-1]) #-1 to remove the id column| vars | n | mean | sd | median | trimmed | mad | min | max | range | skew | kurtosis | se | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| t1_sc | 1 | 3198 | 3.215135 | 0.6949452 | 3 | 3.273828 | 1.4826 | 1 | 4 | 3 | -0.3897593 | -0.6238324 | 0.0122889 |
| t1_intp | 2 | 3162 | 3.024668 | 0.7159051 | 3 | 3.040712 | 0.0000 | 1 | 4 | 3 | -0.1654381 | -0.6702866 | 0.0127314 |
| t1_ext | 3 | 3262 | 1.586143 | 0.6981145 | 1 | 1.476628 | 0.0000 | 1 | 4 | 3 | 1.0780229 | 0.9799469 | 0.0122232 |
| t1_int | 4 | 3231 | 1.545342 | 0.6273211 | 1 | 1.475048 | 0.0000 | 1 | 4 | 3 | 0.8931608 | 0.6133514 | 0.0110362 |
| t2_sc | 5 | 3479 | 3.307560 | 0.7001419 | 3 | 3.389228 | 1.4826 | 1 | 4 | 3 | -0.5764546 | -0.5477944 | 0.0118702 |
| t2_intp | 6 | 3460 | 3.191330 | 0.7131807 | 3 | 3.244220 | 1.4826 | 1 | 4 | 3 | -0.3633800 | -0.7522430 | 0.0121244 |
| t2_ext | 7 | 3482 | 1.646467 | 0.7069649 | 2 | 1.536253 | 1.4826 | 1 | 4 | 3 | 0.9315256 | 0.6705202 | 0.0119807 |
| t2_int | 8 | 3466 | 1.587998 | 0.6323781 | 2 | 1.520908 | 1.4826 | 1 | 4 | 3 | 0.7793179 | 0.4288969 | 0.0107414 |
| t4_sc | 9 | 3619 | 3.284885 | 0.6969259 | 3 | 3.359337 | 1.4826 | 1 | 4 | 3 | -0.5016578 | -0.6670892 | 0.0115849 |
| t4_intp | 10 | 3610 | 3.125208 | 0.7246991 | 3 | 3.163435 | 1.4826 | 1 | 4 | 3 | -0.2816960 | -0.7838474 | 0.0120616 |
| t4_ext | 11 | 3621 | 1.631593 | 0.7122732 | 2 | 1.517777 | 1.4826 | 1 | 4 | 3 | 0.9590041 | 0.6168455 | 0.0118367 |
| t4_int | 12 | 3639 | 1.636988 | 0.6445657 | 2 | 1.567800 | 1.4826 | 1 | 4 | 3 | 0.7215088 | 0.4414765 | 0.0106850 |
#correlation matrix
round(cor(dat[,-1], use = "pairwise.complete"),2)## t1_sc t1_intp t1_ext t1_int t2_sc t2_intp t2_ext t2_int t4_sc t4_intp
## t1_sc 1.00 0.69 -0.59 -0.25 0.55 0.48 -0.48 -0.18 0.34 0.33
## t1_intp 0.69 1.00 -0.48 -0.31 0.47 0.53 -0.39 -0.23 0.29 0.32
## t1_ext -0.59 -0.48 1.00 0.20 -0.51 -0.45 0.64 0.14 -0.38 -0.34
## t1_int -0.25 -0.31 0.20 1.00 -0.18 -0.23 0.15 0.47 -0.10 -0.15
## t2_sc 0.55 0.47 -0.51 -0.18 1.00 0.70 -0.63 -0.26 0.39 0.36
## t2_intp 0.48 0.53 -0.45 -0.23 0.70 1.00 -0.54 -0.31 0.34 0.36
## t2_ext -0.48 -0.39 0.64 0.15 -0.63 -0.54 1.00 0.24 -0.41 -0.36
## t2_int -0.18 -0.23 0.14 0.47 -0.26 -0.31 0.24 1.00 -0.13 -0.16
## t4_sc 0.34 0.29 -0.38 -0.10 0.39 0.34 -0.41 -0.13 1.00 0.72
## t4_intp 0.33 0.32 -0.34 -0.15 0.36 0.36 -0.36 -0.16 0.72 1.00
## t4_ext -0.37 -0.30 0.46 0.09 -0.41 -0.36 0.49 0.11 -0.62 -0.55
## t4_int -0.15 -0.16 0.12 0.20 -0.14 -0.16 0.13 0.24 -0.28 -0.32
## t4_ext t4_int
## t1_sc -0.37 -0.15
## t1_intp -0.30 -0.16
## t1_ext 0.46 0.12
## t1_int 0.09 0.20
## t2_sc -0.41 -0.14
## t2_intp -0.36 -0.16
## t2_ext 0.49 0.13
## t2_int 0.11 0.24
## t4_sc -0.62 -0.28
## t4_intp -0.55 -0.32
## t4_ext 1.00 0.29
## t4_int 0.29 1.00#visusal correlation matrix
corrplot(cor(dat[,-1], use = "pairwise.complete"), order = "original", tl.col='black', tl.cex=.75) 
2 Longitudnal Item Factor Analysis Model
2.1 Configural Invariance model
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 “child behavior” factor for each of the 3 occasions, each time so that the factor is indicated by time-specific interpersonal skills, self-control, and externalizing behaviors.
While the common factor at each of the measurement occasions can be interpreted similarly (e.g., is named “academic achievement”), this model does not impose or assume that the time-specific factors measures the same construct , or that the factors are measured on the same scale.
The diagram would look like 
With some additional thresholds as we have multiple ordered response categories.
Model Specification
configural_invar <- ' #opening quote
#factor loadings
eta1 =~ 1*t1_sc + #for identification
t1_intp +
t1_ext
eta2 =~ 1*t2_sc + #for identification
t2_intp +
t2_ext
eta3 =~ 1*t4_sc + #for identification
t4_intp +
t4_ext
#latent variable variances
eta1~~eta1
eta2~~eta2
eta3~~eta3
#latent variable covariances
eta1~~eta2
eta1~~eta3
eta2~~eta3
#latent variable means
eta1~0*1 #for scaling
eta2~0*1 #for scaling
eta3~0*1 #for scaling
#propensity variances
t1_sc ~~ 1*t1_sc
t1_intp~~ 1*t1_intp
t1_ext ~~ 1*t1_ext
t2_sc ~~ 1*t2_sc
t2_intp~~ 1*t2_intp
t2_ext ~~ 1*t2_ext
t4_sc ~~ 1*t4_sc
t4_intp~~ 1*t4_intp
t4_ext ~~ 1*t4_ext
#unique covariances
#observed variable intercepts/thresholds (4 categories = 3 thresholds)
t1_sc |t1 + t2 + t3
t1_intp |t1 + t2 + t3
t1_ext |t1 + t2 + t3
t2_sc |t1 + t2 + t3
t2_intp |t1 + t2 + t3
t2_ext |t1 + t2 + t3
t4_sc |t1 + t2 + t3
t4_intp |t1 + t2 + t3
t4_ext |t1 + t2 + t3
' #closing quoteModel Estimation and Interpretation
#Model fitting
fit_configural <- cfa(configural_invar, data = dat, parameterization="theta", missing="pairwise",
ordered=c("t1_sc", "t1_intp", "t1_ext", "t1_int",
"t2_sc", "t2_intp", "t2_ext", "t2_int",
"t4_sc", "t4_intp", "t4_ext", "t4_int"))## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan WARNING: some cases are empty and will be ignored:
## 60summary(fit_configural, fit.measures = TRUE)## lavaan 0.6-8 ended normally after 126 iterations
##
## Estimator DWLS
## Optimization method NLMINB
## Number of model parameters 39
##
## Used Total
## Number of observations 3638 3639
## Number of missing patterns 48
##
## Model Test User Model:
## Standard Robust
## Test Statistic 814.553 1339.846
## Degrees of freedom 24 24
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 0.610
## Shift parameter 5.209
## simple second-order correction
##
## Model Test Baseline Model:
##
## Test statistic 84147.713 45237.710
## Degrees of freedom 36 36
## P-value 0.000 0.000
## Scaling correction factor 1.861
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.991 0.971
## Tucker-Lewis Index (TLI) 0.986 0.956
##
## Robust Comparative Fit Index (CFI) NA
## Robust Tucker-Lewis Index (TLI) NA
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.095 0.123
## 90 Percent confidence interval - lower 0.090 0.117
## 90 Percent confidence interval - upper 0.101 0.128
## P-value RMSEA <= 0.05 0.000 0.000
##
## Robust RMSEA NA
## 90 Percent confidence interval - lower NA
## 90 Percent confidence interval - upper NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.060 0.060
##
## Parameter Estimates:
##
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## eta1 =~
## t1_sc 1.000
## t1_intp 0.727 0.045 16.020 0.000
## t1_ext -0.825 0.053 -15.546 0.000
## eta2 =~
## t2_sc 1.000
## t2_intp 0.750 0.047 15.806 0.000
## t2_ext -0.786 0.048 -16.385 0.000
## eta3 =~
## t4_sc 1.000
## t4_intp 0.693 0.063 11.079 0.000
## t4_ext -0.596 0.046 -12.850 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## eta1 ~~
## eta2 3.865 0.262 14.749 0.000
## eta3 3.056 0.262 11.667 0.000
## eta2 ~~
## eta3 3.583 0.299 12.000 0.000
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## eta1 0.000
## eta2 0.000
## eta3 0.000
## .t1_sc 0.000
## .t1_intp 0.000
## .t1_ext 0.000
## .t2_sc 0.000
## .t2_intp 0.000
## .t2_ext 0.000
## .t4_sc 0.000
## .t4_intp 0.000
## .t4_ext 0.000
##
## Thresholds:
## Estimate Std.Err z-value P(>|z|)
## t1_sc|t1 -6.165 0.268 -23.019 0.000
## t1_sc|t2 -2.426 0.104 -23.399 0.000
## t1_sc|t3 0.785 0.056 14.034 0.000
## t1_intp|t1 -4.411 0.151 -29.143 0.000
## t1_intp|t2 -1.361 0.053 -25.848 0.000
## t1_intp|t3 1.171 0.048 24.262 0.000
## t1_ext|t1 0.094 0.044 2.119 0.034
## t1_ext|t2 2.718 0.103 26.349 0.000
## t1_ext|t3 4.224 0.159 26.499 0.000
## t2_sc|t1 -6.516 0.311 -20.983 0.000
## t2_sc|t2 -2.764 0.114 -24.229 0.000
## t2_sc|t3 0.358 0.052 6.827 0.000
## t2_intp|t1 -5.189 0.208 -24.905 0.000
## t2_intp|t2 -1.879 0.064 -29.301 0.000
## t2_intp|t3 0.681 0.044 15.513 0.000
## t2_ext|t1 -0.148 0.043 -3.434 0.001
## t2_ext|t2 2.608 0.092 28.229 0.000
## t2_ext|t3 4.251 0.151 28.122 0.000
## t4_sc|t1 -7.709 0.502 -15.365 0.000
## t4_sc|t2 -3.060 0.179 -17.078 0.000
## t4_sc|t3 0.539 0.063 8.572 0.000
## t4_intp|t1 -5.225 0.198 -26.423 0.000
## t4_intp|t2 -1.756 0.069 -25.566 0.000
## t4_intp|t3 0.920 0.050 18.439 0.000
## t4_ext|t1 -0.054 0.038 -1.416 0.157
## t4_ext|t2 2.334 0.074 31.375 0.000
## t4_ext|t3 3.886 0.120 32.441 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## eta1 4.425 0.395 11.197 0.000
## eta2 5.047 0.456 11.059 0.000
## eta3 6.719 0.857 7.837 0.000
## .t1_sc 1.000
## .t1_intp 1.000
## .t1_ext 1.000
## .t2_sc 1.000
## .t2_intp 1.000
## .t2_ext 1.000
## .t4_sc 1.000
## .t4_intp 1.000
## .t4_ext 1.000
##
## Scales y*:
## Estimate Std.Err z-value P(>|z|)
## t1_sc 0.429
## t1_intp 0.547
## t1_ext 0.499
## t2_sc 0.407
## t2_intp 0.510
## t2_ext 0.493
## t4_sc 0.360
## t4_intp 0.486
## t4_ext 0.543#Model Diagram
semPaths(fit_configural, what="est",
sizeLat = 7, sizeMan = 7, edge.label.cex = .75)
Note that semPaths does not quite know how to draw this kind of model. So, we check the output as well.
2.2 Weak Invariance Model
The weak (or metric) invariance model imposes the additional constraint that the factor loadings are identical at all measurement occasions. The equality constraints on the factor laodings impose that the covariance structures across time are proportional.
The diagram would look like 
Notice in the naming of the factor loadings that the \(\lambda\) are invariant across time.
Model Specification (additional constraints that \(\boldsymbol{\Lambda_t} = \boldsymbol{\Lambda}\))
weak_invar <- ' #opening quote
#factor loadings
eta1 =~ 1*t1_sc + #for identification
lambda_i*t1_intp +
lambda_e*t1_ext
eta2 =~ 1*t2_sc + #for identification
lambda_i*t2_intp +
lambda_e*t2_ext
eta3 =~ 1*t4_sc + #for identification
lambda_i*t4_intp +
lambda_e*t4_ext
#latent variable variances
eta1~~eta1
eta2~~eta2
eta3~~eta3
#latent variable covariances
eta1~~eta2
eta1~~eta3
eta2~~eta3
#latent variable means
eta1~0*1 #for scaling
eta2~0*1 #for scaling
eta3~0*1 #for scaling
# #propensity variances
# t1_sc ~~ 1*t1_sc
# t1_intp~~ 1*t1_intp
# t1_ext ~~ 1*t1_ext
# t2_sc ~~ 1*t2_sc
# t2_intp~~ 1*t2_intp
# t2_ext ~~ 1*t2_ext
# t4_sc ~~ 1*t4_sc
# t4_intp~~ 1*t4_intp
# t4_ext ~~ 1*t4_ext
#unique covariances
#observed variable intercepts
t1_sc |t1 + t2 + t3
t1_intp |t1 + t2 + t3
t1_ext |t1 + t2 + t3
t2_sc |t1 + t2 + t3
t2_intp |t1 + t2 + t3
t2_ext |t1 + t2 + t3
t4_sc |t1 + t2 + t3
t4_intp |t1 + t2 + t3
t4_ext |t1 + t2 + t3
' #closing quoteModel Estimation and Interpretation
#Model Estimation
fit_weak <- cfa(weak_invar, data = dat, parameterization="theta", missing="pairwise",
ordered=c("t1_sc", "t1_intp", "t1_ext", "t1_int",
"t2_sc", "t2_intp", "t2_ext", "t2_int",
"t4_sc", "t4_intp", "t4_ext", "t4_int"))## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan WARNING: some cases are empty and will be ignored:
## 60summary(fit_weak, fit.measures = TRUE)## lavaan 0.6-8 ended normally after 93 iterations
##
## Estimator DWLS
## Optimization method NLMINB
## Number of model parameters 39
## Number of equality constraints 4
##
## Used Total
## Number of observations 3638 3639
## Number of missing patterns 48
##
## Model Test User Model:
## Standard Robust
## Test Statistic 838.621 1198.931
## Degrees of freedom 28 28
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 0.703
## Shift parameter 5.928
## simple second-order correction
##
## Model Test Baseline Model:
##
## Test statistic 84147.713 45237.710
## Degrees of freedom 36 36
## P-value 0.000 0.000
## Scaling correction factor 1.861
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.990 0.974
## Tucker-Lewis Index (TLI) 0.988 0.967
##
## Robust Comparative Fit Index (CFI) NA
## Robust Tucker-Lewis Index (TLI) NA
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.089 0.107
## 90 Percent confidence interval - lower 0.084 0.102
## 90 Percent confidence interval - upper 0.094 0.112
## P-value RMSEA <= 0.05 0.000 0.000
##
## Robust RMSEA NA
## 90 Percent confidence interval - lower NA
## 90 Percent confidence interval - upper NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.060 0.060
##
## Parameter Estimates:
##
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## eta1 =~
## t1_ 1.000
## t1_ (lambda_i) 0.730 0.029 24.773 0.000
## t1_ (lambda_e) -0.744 0.029 -25.850 0.000
## eta2 =~
## t2_ 1.000
## t2_ (lambda_i) 0.730 0.029 24.773 0.000
## t2_ (lambda_e) -0.744 0.029 -25.850 0.000
## eta3 =~
## t4_ 1.000
## t4_ (lambda_i) 0.730 0.029 24.773 0.000
## t4_ (lambda_e) -0.744 0.029 -25.850 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## eta1 ~~
## eta2 4.094 0.260 15.757 0.000
## eta3 2.867 0.196 14.641 0.000
## eta2 ~~
## eta3 3.361 0.221 15.206 0.000
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## eta1 0.000
## eta2 0.000
## eta3 0.000
## .t1_sc 0.000
## .t1_intp 0.000
## .t1_ext 0.000
## .t2_sc 0.000
## .t2_intp 0.000
## .t2_ext 0.000
## .t4_sc 0.000
## .t4_intp 0.000
## .t4_ext 0.000
##
## Thresholds:
## Estimate Std.Err z-value P(>|z|)
## t1_sc|t1 -6.303 0.237 -26.542 0.000
## t1_sc|t2 -2.480 0.085 -29.086 0.000
## t1_sc|t3 0.803 0.055 14.583 0.000
## t1_intp|t1 -4.507 0.140 -32.271 0.000
## t1_intp|t2 -1.390 0.050 -27.530 0.000
## t1_intp|t3 1.196 0.047 25.575 0.000
## t1_ext|t1 0.089 0.042 2.129 0.033
## t1_ext|t2 2.571 0.072 35.577 0.000
## t1_ext|t3 3.995 0.116 34.478 0.000
## t2_sc|t1 -6.674 0.271 -24.591 0.000
## t2_sc|t2 -2.831 0.094 -30.239 0.000
## t2_sc|t3 0.367 0.053 6.863 0.000
## t2_intp|t1 -5.193 0.192 -27.068 0.000
## t2_intp|t2 -1.881 0.056 -33.375 0.000
## t2_intp|t3 0.681 0.043 15.930 0.000
## t2_ext|t1 -0.145 0.042 -3.433 0.001
## t2_ext|t2 2.558 0.072 35.594 0.000
## t2_ext|t3 4.169 0.121 34.519 0.000
## t4_sc|t1 -7.133 0.310 -23.019 0.000
## t4_sc|t2 -2.831 0.093 -30.459 0.000
## t4_sc|t3 0.499 0.054 9.190 0.000
## t4_intp|t1 -5.072 0.167 -30.368 0.000
## t4_intp|t2 -1.704 0.054 -31.587 0.000
## t4_intp|t3 0.893 0.044 20.278 0.000
## t4_ext|t1 -0.060 0.042 -1.415 0.157
## t4_ext|t2 2.569 0.073 35.087 0.000
## t4_ext|t3 4.278 0.122 34.996 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## eta1 4.672 0.309 15.127 0.000
## eta2 5.343 0.361 14.786 0.000
## eta3 5.608 0.363 15.435 0.000
## .t1_sc 1.000
## .t1_intp 1.000
## .t1_ext 1.000
## .t2_sc 1.000
## .t2_intp 1.000
## .t2_ext 1.000
## .t4_sc 1.000
## .t4_intp 1.000
## .t4_ext 1.000
##
## Scales y*:
## Estimate Std.Err z-value P(>|z|)
## t1_sc 0.420
## t1_intp 0.535
## t1_ext 0.528
## t2_sc 0.397
## t2_intp 0.510
## t2_ext 0.502
## t4_sc 0.389
## t4_intp 0.501
## t4_ext 0.493#Model Diagram
# semPaths(fit_weak, what="est",
# sizeLat = 7, sizeMan = 7, edge.label.cex = .75)2.3 Strong Invariance Model
The strong invariance model constrains observed variable intercepts (\(\tau\)), which in this case are the thresholds between ordered categories, to be equal across time, while means of the latent variables are allowed to differ across occasions. Since all mean-level change in the observed variables can be attributed to the factors, the scale of the latent variables is now equal across measurement occasions. Notice that in this set-up with the multiple response options modeled using thresholds, there are no additional resifual variances - so the Strong Invariance Model is also a Strict Factorial Invariance Model (or scalar invariance).
The diagram would look like 
Notice now that the thresholds \(\tau\) are invariant across time.
Model Specification (additional constraints on manifest intercepts, \(\boldsymbol{\tau_t} = \boldsymbol{\tau}\))
strong_invar <- ' #opening quote
#factor loadings
eta1 =~ 1*t1_sc + #for identification
lambda_i*t1_intp +
lambda_e*t1_ext
eta2 =~ 1*t2_sc + #for identification
lambda_i*t2_intp +
lambda_e*t2_ext
eta3 =~ 1*t4_sc + #for identification
lambda_i*t4_intp +
lambda_e*t4_ext
#latent variable variances
eta1~~eta1
eta2~~eta2
eta3~~eta3
#latent variable covariances
eta1~~eta2
eta1~~eta3
eta2~~eta3
#latent variable means
eta1~0*1 #for scaling
eta2~1 #free
eta3~1 #free
# #propensity variances
# t1_sc ~~ 1*t1_sc
# t1_intp~~ 1*t1_intp
# t1_ext ~~ 1*t1_ext
# t2_sc ~~ 1*t2_sc
# t2_intp~~ 1*t2_intp
# t2_ext ~~ 1*t2_ext
# t4_sc ~~ 1*t4_sc
# t4_intp~~ 1*t4_intp
# t4_ext ~~ 1*t4_ext
#unique covariances
#observed variable intercepts (now with equality constraints)
t1_sc |tau_s1*t1 + tau_s2*t2 + tau_s3*t3
t1_intp |tau_i1*t1 + tau_i2*t2 + tau_i3*t3
t1_ext |tau_e1*t1 + tau_e2*t2 + tau_e3*t3
t2_sc |tau_s1*t1 + tau_s2*t2 + tau_s3*t3
t2_intp |tau_i1*t1 + tau_i2*t2 + tau_i3*t3
t2_ext |tau_e1*t1 + tau_e2*t2 + tau_e3*t3
t4_sc |tau_s1*t1 + tau_s2*t2 + tau_s3*t3
t4_intp |tau_i1*t1 + tau_i2*t2 + tau_i3*t3
t4_ext |tau_e1*t1 + tau_e2*t2 + tau_e3*t3
' #closing quoteModel Estimation & Interpretation
fit_strong <- cfa(strong_invar, data = dat, parameterization="theta", missing="pairwise",
ordered=c("t1_sc", "t1_intp", "t1_ext", "t1_int",
"t2_sc", "t2_intp", "t2_ext", "t2_int",
"t4_sc", "t4_intp", "t4_ext", "t4_int"))## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan WARNING: some cases are empty and will be ignored:
## 60summary(fit_strong, fit.measures = TRUE)## lavaan 0.6-8 ended normally after 85 iterations
##
## Estimator DWLS
## Optimization method NLMINB
## Number of model parameters 41
## Number of equality constraints 22
##
## Used Total
## Number of observations 3638 3639
## Number of missing patterns 48
##
## Model Test User Model:
## Standard Robust
## Test Statistic 962.008 1407.889
## Degrees of freedom 44 44
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 0.687
## Shift parameter 7.608
## simple second-order correction
##
## Model Test Baseline Model:
##
## Test statistic 84147.713 45237.710
## Degrees of freedom 36 36
## P-value 0.000 0.000
## Scaling correction factor 1.861
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.989 0.970
## Tucker-Lewis Index (TLI) 0.991 0.975
##
## Robust Comparative Fit Index (CFI) NA
## Robust Tucker-Lewis Index (TLI) NA
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.076 0.092
## 90 Percent confidence interval - lower 0.072 0.088
## 90 Percent confidence interval - upper 0.080 0.096
## P-value RMSEA <= 0.05 0.000 0.000
##
## Robust RMSEA NA
## 90 Percent confidence interval - lower NA
## 90 Percent confidence interval - upper NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.060 0.060
##
## Parameter Estimates:
##
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## eta1 =~
## t1_ 1.000
## t1_ (lambda_i) 0.730 0.030 24.644 0.000
## t1_ (lambda_e) -0.740 0.029 -25.838 0.000
## eta2 =~
## t2_ 1.000
## t2_ (lambda_i) 0.730 0.030 24.644 0.000
## t2_ (lambda_e) -0.740 0.029 -25.838 0.000
## eta3 =~
## t4_ 1.000
## t4_ (lambda_i) 0.730 0.030 24.644 0.000
## t4_ (lambda_e) -0.740 0.029 -25.838 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## eta1 ~~
## eta2 4.156 0.262 15.837 0.000
## eta3 2.857 0.194 14.713 0.000
## eta2 ~~
## eta3 3.363 0.219 15.369 0.000
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## eta1 0.000
## eta2 0.322 0.040 7.995 0.000
## eta3 0.222 0.047 4.706 0.000
## .t1_sc 0.000
## .t1_intp 0.000
## .t1_ext 0.000
## .t2_sc 0.000
## .t2_intp 0.000
## .t2_ext 0.000
## .t4_sc 0.000
## .t4_intp 0.000
## .t4_ext 0.000
##
## Thresholds:
## Estimate Std.Err z-value P(>|z|)
## t1_|1 (t_s1) -6.503 0.193 -33.781 0.000
## t1_|2 (t_s2) -2.537 0.075 -33.638 0.000
## t1_|3 (t_s3) 0.738 0.048 15.360 0.000
## t1_|1 (tau_i1) -4.760 0.109 -43.520 0.000
## t1_|2 (tau_i2) -1.523 0.043 -35.274 0.000
## t1_|3 (tau_i3) 1.053 0.038 27.665 0.000
## t1_|1 (tau_e1) -0.169 0.037 -4.639 0.000
## t1_|2 (tau_e2) 2.421 0.057 42.739 0.000
## t1_|3 (tau_e3) 4.000 0.086 46.341 0.000
## t2_|1 (t_s1) -6.503 0.193 -33.781 0.000
## t2_|2 (t_s2) -2.537 0.075 -33.638 0.000
## t2_|3 (t_s3) 0.738 0.048 15.360 0.000
## t2_|1 (tau_i1) -4.760 0.109 -43.520 0.000
## t2_|2 (tau_i2) -1.523 0.043 -35.274 0.000
## t2_|3 (tau_i3) 1.053 0.038 27.665 0.000
## t2_|1 (tau_e1) -0.169 0.037 -4.639 0.000
## t2_|2 (tau_e2) 2.421 0.057 42.739 0.000
## t2_|3 (tau_e3) 4.000 0.086 46.341 0.000
## t4_|1 (t_s1) -6.503 0.193 -33.781 0.000
## t4_|2 (t_s2) -2.537 0.075 -33.638 0.000
## t4_|3 (t_s3) 0.738 0.048 15.360 0.000
## t4_|1 (tau_i1) -4.760 0.109 -43.520 0.000
## t4_|2 (tau_i2) -1.523 0.043 -35.274 0.000
## t4_|3 (tau_i3) 1.053 0.038 27.665 0.000
## t4_|1 (tau_e1) -0.169 0.037 -4.639 0.000
## t4_|2 (tau_e2) 2.421 0.057 42.739 0.000
## t4_|3 (tau_e3) 4.000 0.086 46.341 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## eta1 4.736 0.302 15.656 0.000
## eta2 5.466 0.348 15.704 0.000
## eta3 5.478 0.336 16.284 0.000
## .t1_sc 1.000
## .t1_intp 1.000
## .t1_ext 1.000
## .t2_sc 1.000
## .t2_intp 1.000
## .t2_ext 1.000
## .t4_sc 1.000
## .t4_intp 1.000
## .t4_ext 1.000
##
## Scales y*:
## Estimate Std.Err z-value P(>|z|)
## t1_sc 0.418
## t1_intp 0.533
## t1_ext 0.528
## t2_sc 0.393
## t2_intp 0.505
## t2_ext 0.501
## t4_sc 0.393
## t4_intp 0.505
## t4_ext 0.500#Model Diagram
# semPaths(fit_strong, what="est",
# sizeLat = 7, sizeMan = 7, edge.label.cex = .75)2.4 Item Response Characteristic Curves
We can use the parameters to construct and look at Item Response Characteristic Curves
This is done by pulling parameters estimates from lavaan (facgtor loading and three thresholds), and then calculating the response characteristic curve. We deomstrate with one item.
Item Characteristic Curve for one item, t1_sc
#pulling parameters from the output
lambda_s <- parameterEstimates(fit_strong)[1,5]
tau_s1 <- parameterEstimates(fit_strong)[19,5]
tau_s2 <- parameterEstimates(fit_strong)[20,5]
tau_s3 <- parameterEstimates(fit_strong)[21,5]
#creating a sequence of ability scores
behavior <- seq(-3,3,.01)
#calculating probability of response curves
t1_scStar = lambda_s*behavior
p123 = pnorm(t1_scStar - tau_s1)
p23 = pnorm(t1_scStar - tau_s2)
p3 = pnorm(t1_scStar - tau_s3)
p0 = 1 - p123
p1 = p123 - p23
p2 = p23 - p3
#putting int oa data frame
temp = data.frame(behavior,p0,p1,p2,p3)
#Plotting Response Characteristic Curve (Item t1_sc)
ggplot(temp, aes(x=behavior, y=p0)) + geom_line(size=2) +
geom_line(data=temp, aes(x=behavior, y=p1), size=2, color="blue") +
geom_line(data=temp, aes(x=behavior, y=p2), size=2, color="red") +
geom_line(data=temp, aes(x=behavior, y=p3), size=2, color="green") +
scale_x_continuous(limits=c(-3,3), breaks=c(-3,-2,-1,0,1,2,3),
name='Child Behavior Latent Variable') +
scale_y_continuous(limits=c(0,1), breaks=c(0,.2,.4,.6,.8,1), name='P(y_i = k)')
We see how the probability of response changes for each category across ability levels. This plot is similar to Figure 15.6 in the book.
2.5 Evaluation of Invariance Models
We can do formal evaluations of how the various forms of invaraince compare.
Model comparison…
Looking at fit statistics
#Compare model fit statistics
round(cbind(configural=inspect(fit_configural, 'fit.measures'),
weak=inspect(fit_weak, 'fit.measures'),
strong=inspect(fit_strong, 'fit.measures')),3)## configural weak strong
## npar 39.000 35.000 19.000
## fmin 0.112 0.115 0.132
## chisq 814.553 838.621 962.008
## df 24.000 28.000 44.000
## pvalue 0.000 0.000 0.000
## chisq.scaled 1339.846 1198.931 1407.889
## df.scaled 24.000 28.000 44.000
## pvalue.scaled 0.000 0.000 0.000
## chisq.scaling.factor 0.610 0.703 0.687
## baseline.chisq 84147.713 84147.713 84147.713
## baseline.df 36.000 36.000 36.000
## baseline.pvalue 0.000 0.000 0.000
## baseline.chisq.scaled 45237.710 45237.710 45237.710
## baseline.df.scaled 36.000 36.000 36.000
## baseline.pvalue.scaled 0.000 0.000 0.000
## baseline.chisq.scaling.factor 1.861 1.861 1.861
## cfi 0.991 0.990 0.989
## tli 0.986 0.988 0.991
## nnfi 0.986 0.988 0.991
## rfi 0.985 0.987 NA
## nfi 0.990 0.990 NA
## pnfi 0.660 0.770 1.208
## ifi 0.991 0.990 0.989
## rni 0.991 0.990 0.989
## cfi.scaled 0.971 0.974 0.970
## tli.scaled 0.956 0.967 0.975
## cfi.robust NA NA NA
## tli.robust NA NA NA
## nnfi.scaled 0.956 0.967 0.975
## nnfi.robust NA NA NA
## rfi.scaled 0.956 0.966 NA
## nfi.scaled 0.970 0.973 NA
## ifi.scaled 0.971 0.974 0.970
## rni.scaled 0.971 0.974 0.970
## rni.robust NA NA NA
## rmsea 0.095 0.089 0.076
## rmsea.ci.lower 0.090 0.084 0.072
## rmsea.ci.upper 0.101 0.094 0.080
## rmsea.pvalue 0.000 0.000 0.000
## rmsea.scaled 0.123 0.107 0.092
## rmsea.ci.lower.scaled 0.117 0.102 0.088
## rmsea.ci.upper.scaled 0.128 0.112 0.096
## rmsea.pvalue.scaled 0.000 0.000 0.000
## rmsea.robust NA NA NA
## rmsea.ci.lower.robust NA NA NA
## rmsea.ci.upper.robust NA NA NA
## rmsea.pvalue.robust NA NA NA
## rmr 0.055 0.055 0.094
## rmr_nomean 0.060 0.060 0.060
## srmr 0.060 0.060 0.060
## srmr_bentler 0.055 0.055 0.094
## srmr_bentler_nomean 0.060 0.060 0.060
## crmr 0.060 0.060 0.103
## crmr_nomean 0.068 0.067 0.067
## srmr_mplus NA NA NA
## srmr_mplus_nomean NA NA NA
## cn_05 163.594 180.274 229.656
## cn_01 192.906 210.377 260.765
## gfi 0.993 0.992 0.991
## agfi 0.981 0.983 0.988
## pgfi 0.378 0.441 0.692
## mfi 0.897 0.895 0.881Relatively worse fit of the more constrained models is 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. With strong incentive to study change, we may conclude that the Strong Invariance model fits the data reasonably well.
Formal statistical improvements/decrements in model fit that accompany each set of constraints can be examined using chi-square difference tests.
#Chi-square difference test for nested models
anova(fit_configural, fit_weak, fit_strong)| Df | AIC | BIC | Chisq | Chisq diff | Df diff | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|
| fit_configural | 24 | NA | NA | 814.5535 | NA | NA | NA |
| fit_weak | 28 | NA | NA | 838.6211 | 23.68773 | 4 | 9.23e-05 |
| fit_strong | 44 | NA | NA | 962.0083 | 194.38052 | 16 | 0.00e+00 |
#lavTestLRT(fit_configural, fit_weak, fit_strong) #does the same as anovaWith these data, all three pairs of tests show significant difference in chi-square. In all cases, we reject the null hypothesis that the models fit the same. The more constrained models fit worse than the less-constrained models. That suggests that there is NOT strong factorial invariance.
However, remember that we often have incentives to go on towards the examination of change. Given our goals, we may accept the evidence of strong invariance as sufficient (it is still a good fitting model) and carry on with the analysis of factor-level change.
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. That is, we impose invariance and move forward.
3 Second-order Growth Model
Once strong or strict invariance is supported or imposed for the longitudinal factor model, we can examine changes in the latent factor scores.
Working from the strong/strict factorial invariance model, we expand to model change in children’s latent academic achievement scores using a second-order latent basis growth model. Note that a few other adjustments have to be made for identification and scaling. Important to make sure the model matches one’s diagram as intended.
The diagram would look like 
There are both 1st order factors that establish the construct of interest and 2nd order factors that model change in that 1st order factor. Very cool!
Model specification
growth_invar <- ' #opening quote
#factor loadings
eta1 =~ 1*t1_sc + #for identification
lambda_i*t1_intp +
lambda_e*t1_ext
eta2 =~ 1*t2_sc + #for identification
lambda_i*t2_intp +
lambda_e*t2_ext
eta3 =~ 1*t4_sc + #for identification
lambda_i*t4_intp +
lambda_e*t4_ext
#latent variable variances
eta1~~psi*eta1 #equality constraints
eta2~~psi*eta2
eta3~~psi*eta3
#latent variable covariances
eta1~~0*eta2 #fixed to zero
eta1~~0*eta3
eta2~~0*eta3
#latent variable means
eta1~0*1 #fixed to zero
eta2~0*1
eta3~0*1
# #propensity variances
# t1_sc ~~ 1*t1_sc
# t1_intp~~ 1*t1_intp
# t1_ext ~~ 1*t1_ext
# t2_sc ~~ 1*t2_sc
# t2_intp~~ 1*t2_intp
# t2_ext ~~ 1*t2_ext
# t4_sc ~~ 1*t4_sc
# t4_intp~~ 1*t4_intp
# t4_ext ~~ 1*t4_ext
#unique covariances
#observed variable intercepts
t1_sc |tau_s1*t1 + tau_s2*t2 + tau_s3*t3
t1_intp |tau_i1*t1 + tau_i2*t2 + tau_i3*t3
t1_ext |tau_e1*t1 + tau_e2*t2 + tau_e3*t3
t2_sc |tau_s1*t1 + tau_s2*t2 + tau_s3*t3
t2_intp |tau_i1*t1 + tau_i2*t2 + tau_i3*t3
t2_ext |tau_e1*t1 + tau_e2*t2 + tau_e3*t3
t4_sc |tau_s1*t1 + tau_s2*t2 + tau_s3*t3
t4_intp |tau_i1*t1 + tau_i2*t2 + tau_i3*t3
t4_ext |tau_e1*t1 + tau_e2*t2 + tau_e3*t3
#second-order latent basis growth
#growth factors
xi_1 =~ 1*eta1+ #intercept factor
1*eta2+
1*eta3
xi_2 =~ 0*eta1+ #latent basis slope factor
1*eta2+
3*eta3
#factor variances & covariance
xi_1~~xi_1
xi_2~~xi_2
xi_1~~xi_2
#factor intercepts
xi_1~0*1 #for scaling
xi_2~1
' #closing quoteModel Estimation & Interpretation
#Model estimation
fit_growth <- cfa(growth_invar, data = dat, parameterization="theta",
ordered=c("t1_sc", "t1_intp", "t1_ext", "t1_int",
"t2_sc", "t2_intp", "t2_ext", "t2_int",
"t4_sc", "t4_intp", "t4_ext", "t4_int"))
summary(fit_growth, fit.measures = TRUE)## lavaan 0.6-8 ended normally after 61 iterations
##
## Estimator DWLS
## Optimization method NLMINB
## Number of model parameters 40
## Number of equality constraints 24
##
## Used Total
## Number of observations 2922 3639
##
## Model Test User Model:
## Standard Robust
## Test Statistic 967.785 1270.049
## Degrees of freedom 47 47
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 0.768
## Shift parameter 9.356
## simple second-order correction
##
## Model Test Baseline Model:
##
## Test statistic 74458.518 38959.199
## Degrees of freedom 36 36
## P-value 0.000 0.000
## Scaling correction factor 1.912
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.988 0.969
## Tucker-Lewis Index (TLI) 0.991 0.976
##
## Robust Comparative Fit Index (CFI) NA
## Robust Tucker-Lewis Index (TLI) NA
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.082 0.094
## 90 Percent confidence interval - lower 0.077 0.090
## 90 Percent confidence interval - upper 0.086 0.099
## P-value RMSEA <= 0.05 0.000 0.000
##
## Robust RMSEA NA
## 90 Percent confidence interval - lower NA
## 90 Percent confidence interval - upper NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.062 0.062
##
## Parameter Estimates:
##
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## eta1 =~
## t1_ 1.000
## t1_ (lambda_i) 0.729 0.031 23.251 0.000
## t1_ (lambda_e) -0.757 0.032 -23.672 0.000
## eta2 =~
## t2_ 1.000
## t2_ (lambda_i) 0.729 0.031 23.251 0.000
## t2_ (lambda_e) -0.757 0.032 -23.672 0.000
## eta3 =~
## t4_ 1.000
## t4_ (lambda_i) 0.729 0.031 23.251 0.000
## t4_ (lambda_e) -0.757 0.032 -23.672 0.000
## xi_1 =~
## et1 1.000
## et2 1.000
## et3 1.000
## xi_2 =~
## et1 0.000
## et2 1.000
## et3 3.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## .eta1 ~~
## .eta2 0.000
## .eta3 0.000
## .eta2 ~~
## .eta3 0.000
## xi_1 ~~
## xi_2 -0.606 0.066 -9.121 0.000
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## .eta1 0.000
## .eta2 0.000
## .eta3 0.000
## xi_1 0.000
## xi_2 0.041 0.017 2.382 0.017
## .t1_sc 0.000
## .t1_intp 0.000
## .t1_ext 0.000
## .t2_sc 0.000
## .t2_intp 0.000
## .t2_ext 0.000
## .t4_sc 0.000
## .t4_intp 0.000
## .t4_ext 0.000
##
## Thresholds:
## Estimate Std.Err z-value P(>|z|)
## t1_|1 (t_s1) -6.658 0.215 -30.916 0.000
## t1_|2 (t_s2) -2.678 0.085 -31.503 0.000
## t1_|3 (t_s3) 0.613 0.052 11.896 0.000
## t1_|1 (tau_i1) -4.918 0.124 -39.657 0.000
## t1_|2 (tau_i2) -1.629 0.048 -33.931 0.000
## t1_|3 (tau_i3) 0.939 0.041 23.101 0.000
## t1_|1 (tau_e1) -0.107 0.041 -2.594 0.009
## t1_|2 (tau_e2) 2.542 0.066 38.813 0.000
## t1_|3 (tau_e3) 4.153 0.099 42.113 0.000
## t2_|1 (t_s1) -6.658 0.215 -30.916 0.000
## t2_|2 (t_s2) -2.678 0.085 -31.503 0.000
## t2_|3 (t_s3) 0.613 0.052 11.896 0.000
## t2_|1 (tau_i1) -4.918 0.124 -39.657 0.000
## t2_|2 (tau_i2) -1.629 0.048 -33.931 0.000
## t2_|3 (tau_i3) 0.939 0.041 23.101 0.000
## t2_|1 (tau_e1) -0.107 0.041 -2.594 0.009
## t2_|2 (tau_e2) 2.542 0.066 38.813 0.000
## t2_|3 (tau_e3) 4.153 0.099 42.113 0.000
## t4_|1 (t_s1) -6.658 0.215 -30.916 0.000
## t4_|2 (t_s2) -2.678 0.085 -31.503 0.000
## t4_|3 (t_s3) 0.613 0.052 11.896 0.000
## t4_|1 (tau_i1) -4.918 0.124 -39.657 0.000
## t4_|2 (tau_i2) -1.629 0.048 -33.931 0.000
## t4_|3 (tau_i3) 0.939 0.041 23.101 0.000
## t4_|1 (tau_e1) -0.107 0.041 -2.594 0.009
## t4_|2 (tau_e2) 2.542 0.066 38.813 0.000
## t4_|3 (tau_e3) 4.153 0.099 42.113 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .eta1 (psi) 0.922 0.077 12.042 0.000
## .eta2 (psi) 0.922 0.077 12.042 0.000
## .eta3 (psi) 0.922 0.077 12.042 0.000
## xi_1 4.828 0.338 14.282 0.000
## xi_2 0.357 0.031 11.600 0.000
## .t1_sc 1.000
## .t1_intp 1.000
## .t1_ext 1.000
## .t2_sc 1.000
## .t2_intp 1.000
## .t2_ext 1.000
## .t4_sc 1.000
## .t4_intp 1.000
## .t4_ext 1.000
##
## Scales y*:
## Estimate Std.Err z-value P(>|z|)
## t1_sc 0.385
## t1_intp 0.496
## t1_ext 0.483
## t2_sc 0.412
## t2_intp 0.527
## t2_ext 0.513
## t4_sc 0.398
## t4_intp 0.511
## t4_ext 0.497#Model Diagram
semPaths(fit_growth, what="est",
sizeLat = 7, sizeMan = 7, edge.label.cex = .75)
Note that the diagram is not quite accurate - but does give us a general picture of the model we are working with. Helpful to re-draw by hand and map the output.
3.1 Plotting Factor Score Trajectories
Finally, we can visualize the predicted factor score trajectories by extracting the predicted factor estimates of the second order growth model and plotting.
Plotting factor trajectories.
#extract predicted factor scores
pred<-as.data.frame(lavPredict(fit_growth))
#define function to fit factor trajectories
fit.line <- function(xi1,xi2,t) {
xi1 + xi2*t
}
time<-c(0,1,3)
#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("pred0","pred1","pred3")
out$ID<-1:nrow(out)
#reshape estimated scores to long format (for plotting)
outlong<-reshape(out,
varying=c("pred0","pred1","pred3"),
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, 0.041,time))
avg$grade<-c(0,1,3)
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 Child Behavior Factor Score")+
geom_line() +
geom_line(data=avg,aes(x=grade,y=fit,group=ID),color="red",size=2) +
geom_hline(aes(yintercept=0), color="red")
That was fun! The ability to incorporate item-level factor analysis within the growth modeling framework provides new opportunites for precision and the kinds of data that can be used for investigation and study of interindividual differences in intraindivudal change. Very exciting!
As always, be careful! Thanks for playing!