## Overview

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.

## Outline

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.

1. Read in the data and call needed libraries.
2. Bivariate Dual Change Model.

## Step 0: Read in the data and call needed libraries.

#set filepath
filepath <- "https://quantdev.ssri.psu.edu/sites/qdev/files/nlsy_data_subset.csv"

#read in the .csv file using the url() function
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.

## Step 1: Bivariate Dual Change Model.

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

#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

dr2 ~ delta_r*lm1
dr3 ~ delta_r*lm2
dr4 ~ delta_r*lm3
dr5 ~ delta_r*lm4
dr6 ~ delta_r*lm5
dr7 ~ delta_r*lm6

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% coverage
summary(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.007

The 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.