This tutorial walks through the fitting of a bivariate latent change score model in the structural equation modeling framework using the lavaan package. In this tutorial, we will be using a sample data set that includes repeated measures of children’s math and reading scores from the second through eighth grade from the NLSY-CYA (Center for Human Resource Research, 2009).
The example provided in this tutorial is from Chapter 17 of Grimm, Ram, and Estabrook (2016). Please refer to the chapter for further interpretations and insights about the analyses.
This tutorial provides line-by-line code to conduct a bivariate dual change model, which allows researchers to examine time-sequential associations between two variables/processes.
#set filepath
filepath <- "https://quantdev.ssri.psu.edu/sites/qdev/files/nlsy_data_subset.csv"
#read in the .csv file using the url() function
nlsy_data <- read.csv(file=url(filepath),header=TRUE)
head(nlsy_data)## X id math2 math3 math4 math5 math6 math7 math8 rec2 rec3 rec4 rec5
## 1 1 201 NA 38 NA 55 NA NA NA NA 35 NA 52
## 2 2 303 26 NA NA 33 NA NA NA 26 NA NA 35
## 3 3 2702 56 NA 58 NA NA NA 80 33 NA 50 NA
## 4 4 4303 NA 41 58 NA NA NA NA NA 34 43 NA
## 5 5 5002 NA NA 46 NA 54 NA 66 NA NA 46 NA
## 6 6 5005 35 NA 50 NA 60 NA 59 47 NA 67 NA
## rec6 rec7 rec8
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA 66
## 4 NA NA NA
## 5 70 NA 77
## 6 75 NA 78#calling the libraries we will need throughout this tutorial
library(lavaan)## This is lavaan 0.5-22## lavaan is BETA software! Please report any bugs.In this example, we are fitting a bivariate dual change score model to examine the coupling of mathematics and reading scores over time.
bdcm.lavaan <- '
#mathematics
#latent true scores
lm1 =~ 1*math2
lm2 =~ 1*math3
lm3 =~ 1*math4
lm4 =~ 1*math5
lm5 =~ 1*math6
lm6 =~ 1*math7
lm7 =~ 1*math8
#latent true score means
lm1 ~ 1
lm2 ~ 0*1
lm3 ~ 0*1
lm4 ~ 0*1
lm5 ~ 0*1
lm6 ~ 0*1
lm7 ~ 0*1
#latent true score variances
lm1 ~~ start(15)*lm1
lm2 ~~ 0*lm2
lm3 ~~ 0*lm3
lm4 ~~ 0*lm4
lm5 ~~ 0*lm5
lm6 ~~ 0*lm6
lm7 ~~ 0*lm7
#observed intercepts (fixed to zero)
math2 ~ 0*1
math3 ~ 0*1
math4 ~ 0*1
math5 ~ 0*1
math6 ~ 0*1
math7 ~ 0*1
math8 ~ 0*1
#observed residual variances
math2 ~~ sigma2_u*math2
math3 ~~ sigma2_u*math3
math4 ~~ sigma2_u*math4
math5 ~~ sigma2_u*math5
math6 ~~ sigma2_u*math6
math7 ~~ sigma2_u*math7
math8 ~~ sigma2_u*math8
#autoregressions
lm2 ~ 1*lm1
lm3 ~ 1*lm2
lm4 ~ 1*lm3
lm5 ~ 1*lm4
lm6 ~ 1*lm5
lm7 ~ 1*lm6
#latent change scores
dm2 =~ 1*lm2
dm3 =~ 1*lm3
dm4 =~ 1*lm4
dm5 =~ 1*lm5
dm6 =~ 1*lm6
dm7 =~ 1*lm7
#latent change score means
dm2 ~ 0*1
dm3 ~ 0*1
dm4 ~ 0*1
dm5 ~ 0*1
dm6 ~ 0*1
dm7 ~ 0*1
#latent change score variances
dm2 ~~ 0*dm2
dm3 ~~ 0*dm3
dm4 ~~ 0*dm4
dm5 ~~ 0*dm5
dm6 ~~ 0*dm6
dm7 ~~ 0*dm7
#constant change factor
g2 =~ 1*dm2 +
1*dm3 +
1*dm4 +
1*dm5 +
1*dm6 +
1*dm7
#constant change factor mean
g2 ~ start(15)*1
#constant change factor variance
g2 ~~ g2
#constant change factor covariance with the initial true score
g2 ~~ lm1
#proportional effects
dm2 ~ start(-.2)*pi_m * lm1
dm3 ~ start(-.2)*pi_m * lm2
dm4 ~ start(-.2)*pi_m * lm3
dm5 ~ start(-.2)*pi_m * lm4
dm6 ~ start(-.2)*pi_m * lm5
dm7 ~ start(-.2)*pi_m * lm6
#reading
#latent true scores
lr1 =~ 1*rec2
lr2 =~ 1*rec3
lr3 =~ 1*rec4
lr4 =~ 1*rec5
lr5 =~ 1*rec6
lr6 =~ 1*rec7
lr7 =~ 1*rec8
#latent true score means
lr1 ~ 1
lr2 ~ 0*1
lr3 ~ 0*1
lr4 ~ 0*1
lr5 ~ 0*1
lr6 ~ 0*1
lr7 ~ 0*1
#latent true score variances
lr1 ~~ start(15)*lr1
lr2 ~~ 0*lr2
lr3 ~~ 0*lr3
lr4 ~~ 0*lr4
lr5 ~~ 0*lr5
lr6 ~~ 0*lr6
lr7 ~~ 0*lr7
#observed intercept variances (fixed to zero)
rec2 ~ 0*1
rec3 ~ 0*1
rec4 ~ 0*1
rec5 ~ 0*1
rec6 ~ 0*1
rec7 ~ 0*1
rec8 ~ 0*1
#observed residual variances
rec2 ~~ sigma2_s*rec2
rec3 ~~ sigma2_s*rec3
rec4 ~~ sigma2_s*rec4
rec5 ~~ sigma2_s*rec5
rec6 ~~ sigma2_s*rec6
rec7 ~~ sigma2_s*rec7
rec8 ~~ sigma2_s*rec8
#autoregressions
lr2 ~ 1*lr1
lr3 ~ 1*lr2
lr4 ~ 1*lr3
lr5 ~ 1*lr4
lr6 ~ 1*lr5
lr7 ~ 1*lr6
#latent change scores
dr2 =~ 1*lr2
dr3 =~ 1*lr3
dr4 =~ 1*lr4
dr5 =~ 1*lr5
dr6 =~ 1*lr6
dr7 =~ 1*lr7
#latent change score means
dr2 ~ 0*1
dr3 ~ 0*1
dr4 ~ 0*1
dr5 ~ 0*1
dr6 ~ 0*1
dr7 ~ 0*1
#latent change score variances
dr2 ~~ 0*dr2
dr3 ~~ 0*dr3
dr4 ~~ 0*dr4
dr5 ~~ 0*dr5
dr6 ~~ 0*dr6
dr7 ~~ 0*dr7
#constant change factor
j2 =~ 1*dr2 +
1*dr3 +
1*dr4 +
1*dr5 +
1*dr6 +
1*dr7
#constant change factor mean
j2 ~ start(10)*1
#constant change factor variance
j2 ~~ j2
#constant change factor covariance with the initial true score
j2 ~~ lr1
#proportional effects
dr2 ~ start(-.2)*pi_r * lr1
dr3 ~ start(-.2)*pi_r * lr2
dr4 ~ start(-.2)*pi_r * lr3
dr5 ~ start(-.2)*pi_r * lr4
dr6 ~ start(-.2)*pi_r * lr5
dr7 ~ start(-.2)*pi_r * lr6
#bivariate information
#covariances between the latent growth factors
lm1 ~~ lr1
g2 ~~ j2
lm1 ~~ j2
lr1 ~~ g2
#residual covariances
math2 ~~ sigma_su*rec2
math3 ~~ sigma_su*rec3
math4 ~~ sigma_su*rec4
math5 ~~ sigma_su*rec5
math6 ~~ sigma_su*rec6
math7 ~~ sigma_su*rec7
#coupling parameters
#math to reading
dr2 ~ delta_r*lm1
dr3 ~ delta_r*lm2
dr4 ~ delta_r*lm3
dr5 ~ delta_r*lm4
dr6 ~ delta_r*lm5
dr7 ~ delta_r*lm6
#reading to math
dm2 ~ delta_m*lr1
dm3 ~ delta_m*lr2
dm4 ~ delta_m*lr3
dm5 ~ delta_m*lr4
dm6 ~ delta_m*lr5
dm7 ~ delta_m*lr6'
lavaan.fit <- lavaan(bdcm.lavaan,
data = nlsy_data,
meanstructure = TRUE,
estimator = "ML",
missing = "fiml",
fixed.x = FALSE,
mimic="mplus",
control=list(iter.max=500),
verbose=FALSE)## Warning in lav_data_full(data = data, group = group, group.label = group.label, : lavaan WARNING: some cases are empty and will be ignored:
## 741## Warning in lav_data_full(data = data, group = group, group.label =
## group.label, : lavaan WARNING: due to missing values, some pairwise
## combinations have less than 10% coveragesummary(lavaan.fit, fit.measures=TRUE)## lavaan (0.5-22) converged normally after 182 iterations
##
## Used Total
## Number of observations 932 933
##
## Number of missing patterns 66
##
## Estimator ML
## Minimum Function Test Statistic 166.098
## Degrees of freedom 98
## P-value (Chi-square) 0.000
##
## Model test baseline model:
##
## Minimum Function Test Statistic 2710.581
## Degrees of freedom 91
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.974
## Tucker-Lewis Index (TLI) 0.976
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -15680.484
## Loglikelihood unrestricted model (H1) -15597.435
##
## Number of free parameters 21
## Akaike (AIC) 31402.968
## Bayesian (BIC) 31504.552
## Sample-size adjusted Bayesian (BIC) 31437.857
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.027
## 90 Percent Confidence Interval 0.020 0.034
## P-value RMSEA <= 0.05 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.111
##
## Parameter Estimates:
##
## Information Observed
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## lm1 =~
## math2 1.000
## lm2 =~
## math3 1.000
## lm3 =~
## math4 1.000
## lm4 =~
## math5 1.000
## lm5 =~
## math6 1.000
## lm6 =~
## math7 1.000
## lm7 =~
## math8 1.000
## dm2 =~
## lm2 1.000
## dm3 =~
## lm3 1.000
## dm4 =~
## lm4 1.000
## dm5 =~
## lm5 1.000
## dm6 =~
## lm6 1.000
## dm7 =~
## lm7 1.000
## g2 =~
## dm2 1.000
## dm3 1.000
## dm4 1.000
## dm5 1.000
## dm6 1.000
## dm7 1.000
## lr1 =~
## rec2 1.000
## lr2 =~
## rec3 1.000
## lr3 =~
## rec4 1.000
## lr4 =~
## rec5 1.000
## lr5 =~
## rec6 1.000
## lr6 =~
## rec7 1.000
## lr7 =~
## rec8 1.000
## dr2 =~
## lr2 1.000
## dr3 =~
## lr3 1.000
## dr4 =~
## lr4 1.000
## dr5 =~
## lr5 1.000
## dr6 =~
## lr6 1.000
## dr7 =~
## lr7 1.000
## j2 =~
## dr2 1.000
## dr3 1.000
## dr4 1.000
## dr5 1.000
## dr6 1.000
## dr7 1.000
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## lm2 ~
## lm1 1.000
## lm3 ~
## lm2 1.000
## lm4 ~
## lm3 1.000
## lm5 ~
## lm4 1.000
## lm6 ~
## lm5 1.000
## lm7 ~
## lm6 1.000
## dm2 ~
## lm1 (pi_m) -0.293 0.115 -2.560 0.010
## dm3 ~
## lm2 (pi_m) -0.293 0.115 -2.560 0.010
## dm4 ~
## lm3 (pi_m) -0.293 0.115 -2.560 0.010
## dm5 ~
## lm4 (pi_m) -0.293 0.115 -2.560 0.010
## dm6 ~
## lm5 (pi_m) -0.293 0.115 -2.560 0.010
## dm7 ~
## lm6 (pi_m) -0.293 0.115 -2.560 0.010
## lr2 ~
## lr1 1.000
## lr3 ~
## lr2 1.000
## lr4 ~
## lr3 1.000
## lr5 ~
## lr4 1.000
## lr6 ~
## lr5 1.000
## lr7 ~
## lr6 1.000
## dr2 ~
## lr1 (pi_r) -0.495 0.115 -4.308 0.000
## dr3 ~
## lr2 (pi_r) -0.495 0.115 -4.308 0.000
## dr4 ~
## lr3 (pi_r) -0.495 0.115 -4.308 0.000
## dr5 ~
## lr4 (pi_r) -0.495 0.115 -4.308 0.000
## dr6 ~
## lr5 (pi_r) -0.495 0.115 -4.308 0.000
## dr7 ~
## lr6 (pi_r) -0.495 0.115 -4.308 0.000
## dr2 ~
## lm1 (dlt_r) 0.391 0.129 3.029 0.002
## dr3 ~
## lm2 (dlt_r) 0.391 0.129 3.029 0.002
## dr4 ~
## lm3 (dlt_r) 0.391 0.129 3.029 0.002
## dr5 ~
## lm4 (dlt_r) 0.391 0.129 3.029 0.002
## dr6 ~
## lm5 (dlt_r) 0.391 0.129 3.029 0.002
## dr7 ~
## lm6 (dlt_r) 0.391 0.129 3.029 0.002
## dm2 ~
## lr1 (dlt_m) 0.053 0.098 0.540 0.589
## dm3 ~
## lr2 (dlt_m) 0.053 0.098 0.540 0.589
## dm4 ~
## lr3 (dlt_m) 0.053 0.098 0.540 0.589
## dm5 ~
## lr4 (dlt_m) 0.053 0.098 0.540 0.589
## dm6 ~
## lr5 (dlt_m) 0.053 0.098 0.540 0.589
## dm7 ~
## lr6 (dlt_m) 0.053 0.098 0.540 0.589
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## lm1 ~~
## g2 14.165 2.164 6.545 0.000
## lr1 ~~
## j2 26.323 3.483 7.557 0.000
## lm1 ~~
## lr1 56.494 5.487 10.297 0.000
## g2 ~~
## j2 0.681 2.552 0.267 0.790
## lm1 ~~
## j2 6.163 3.103 1.986 0.047
## g2 ~~
## lr1 9.791 3.083 3.176 0.001
## .math2 ~~
## .rec2 (sgm_) 7.008 1.259 5.566 0.000
## .math3 ~~
## .rec3 (sgm_) 7.008 1.259 5.566 0.000
## .math4 ~~
## .rec4 (sgm_) 7.008 1.259 5.566 0.000
## .math5 ~~
## .rec5 (sgm_) 7.008 1.259 5.566 0.000
## .math6 ~~
## .rec6 (sgm_) 7.008 1.259 5.566 0.000
## .math7 ~~
## .rec7 (sgm_) 7.008 1.259 5.566 0.000
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## lm1 32.543 0.446 72.980 0.000
## .lm2 0.000
## .lm3 0.000
## .lm4 0.000
## .lm5 0.000
## .lm6 0.000
## .lm7 0.000
## .math2 0.000
## .math3 0.000
## .math4 0.000
## .math5 0.000
## .math6 0.000
## .math7 0.000
## .math8 0.000
## .dm2 0.000
## .dm3 0.000
## .dm4 0.000
## .dm5 0.000
## .dm6 0.000
## .dm7 0.000
## g2 15.092 0.926 16.291 0.000
## lr1 34.386 0.433 79.425 0.000
## .lr2 0.000
## .lr3 0.000
## .lr4 0.000
## .lr5 0.000
## .lr6 0.000
## .lr7 0.000
## .rec2 0.000
## .rec3 0.000
## .rec4 0.000
## .rec5 0.000
## .rec6 0.000
## .rec7 0.000
## .rec8 0.000
## .dr2 0.000
## .dr3 0.000
## .dr4 0.000
## .dr5 0.000
## .dr6 0.000
## .dr7 0.000
## j2 10.886 0.934 11.650 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## lm1 72.886 6.603 11.038 0.000
## .lm2 0.000
## .lm3 0.000
## .lm4 0.000
## .lm5 0.000
## .lm6 0.000
## .lm7 0.000
## .math2 (sigm2_) 31.306 1.609 19.460 0.000
## .math3 (sigm2_) 31.306 1.609 19.460 0.000
## .math4 (sigm2_) 31.306 1.609 19.460 0.000
## .math5 (sigm2_) 31.306 1.609 19.460 0.000
## .math6 (sigm2_) 31.306 1.609 19.460 0.000
## .math7 (sigm2_) 31.306 1.609 19.460 0.000
## .math8 (sigm2_) 31.306 1.609 19.460 0.000
## .dm2 0.000
## .dm3 0.000
## .dm4 0.000
## .dm5 0.000
## .dm6 0.000
## .dm7 0.000
## g2 5.743 1.697 3.384 0.001
## lr1 71.978 7.273 9.896 0.000
## .lr2 0.000
## .lr3 0.000
## .lr4 0.000
## .lr5 0.000
## .lr6 0.000
## .lr7 0.000
## .rec2 (sgm2_s) 33.217 1.763 18.844 0.000
## .rec3 (sgm2_s) 33.217 1.763 18.844 0.000
## .rec4 (sgm2_s) 33.217 1.763 18.844 0.000
## .rec5 (sgm2_s) 33.217 1.763 18.844 0.000
## .rec6 (sgm2_s) 33.217 1.763 18.844 0.000
## .rec7 (sgm2_s) 33.217 1.763 18.844 0.000
## .rec8 (sgm2_s) 33.217 1.763 18.844 0.000
## .dr2 0.000
## .dr3 0.000
## .dr4 0.000
## .dr5 0.000
## .dr6 0.000
## .dr7 0.000
## j2 18.384 6.805 2.701 0.007The warnings that occur when fitting this model are related to missing data.
We find that the mean trajectories for mathematics begins at 32.4 (intercept of lm1) and for reading begins at 34.39 (intercept of lr1). For mathematics, the scores change at a constant rate of 15.09 points a year (intercept of g2) and has a proportional change effect of -0.29 (pi_m). For reading, the scores change at a constant rate of 10.89 (intercept of j2) and has a proportional change effect of -0.50 (pi_r).
Additionally, the no coupling and two unidirectional coupling models can be conducted by modifying (i.e., removing parameters) from the above model.