1.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.2 Prelim - Reading in Repeated Measures Data

We use data from the NLSY-CYA (Center for Human Resource Research, 2009) that includes repeated measures of children’s math ability (math) and reading comprehension ability (rec) from the second through eighth grade.

Reading in the data

#set filepath
filepath <- "https://raw.githubusercontent.com/LRI-2/Data/main/GrowthModeling/nlsy_math_hyp_wide_R.dat"
#read in the text data file using the url() function
nlsy_data <- read.table(file=url(filepath),na.strings = ".")

#adding names for the columns of the data set
names(nlsy_data) <- c('id', 'female', 'lb_wght', 'anti_k1', 'math2', 'math3', 'math4', 'math5', 'math6', 'math7', 'math8', 'comp2', 'comp3', 'comp4', 'comp5', 'comp6', 'comp7', 'comp8', 'rec2', 'rec3', 'rec4', 'rec5', 'rec6', 'rec7', 'rec8', 'bpi2', 'bpi3', 'bpi4', 'bpi5', 'bpi6', 'bpi7', 'bpi8', 'asl2', 'asl3', 'asl4', 'asl5', 'asl6', 'asl7', 'asl8', 'ax2', 'ax3', 'ax4', 'ax5', 'ax6', 'ax7', 'ax8', 'hds2', 'hds3', 'hds4', 'hds5', 'hds6', 'hds7', 'hds8', 'hyp2', 'hyp3', 'hyp4', 'hyp5', 'hyp6', 'hyp7', 'hyp8', 'dpn2', 'dpn3', 'dpn4', 'dpn5', 'dpn6', 'dpn7', 'dpn8', 'wdn2', 'wdn3', 'wdn4', 'wdn5', 'wdn6', 'wdn7', 'wdn8', 'age2', 'age3', 'age4', 'age5', 'age6', 'age7', 'age8', 'men2', 'men3', 'men4', 'men5', 'men6', 'men7', 'men8', 'spring2', 'spring3', 'spring4', 'spring5', 'spring6', 'spring7', 'spring8', 'anti2', 'anti3', 'anti4', 'anti5', 'anti6', 'anti7', 'anti8')

#reduce data down to the id variable and the math and reading variables of interest
nlsy_data <- nlsy_data[ ,c('id', 'math2', 'math3', 'math4', 'math5', 'math6', 'math7', 'math8',
                           'rec2', 'rec3', 'rec4', 'rec5', 'rec6', 'rec7', 'rec8')]

describe(nlsy_data)

1.0.3 Prelim - Plotting the Repeated Measures Data

#reshaping wide to wide
data_long <- reshape(data=nlsy_data,
                    varying = c('math2', 'math3', 'math4', 'math5', 'math6', 'math7', 'math8',
                           'rec2', 'rec3', 'rec4', 'rec5', 'rec6', 'rec7', 'rec8'),
                    timevar=c("grade"), 
                    idvar=c("id"),
                    direction="long", sep="")

#sorting for easy viewing
#reorder by id and day
data_long <- data_long[order(data_long$id,data_long$grade), ]

#looking at the long data
head(data_long, 8)

#Plotting intraindividual change MATH
ggplot(data = data_long, aes(x = grade, y = math, group = id)) +
  geom_point(color="blue") + 
  geom_line(color="blue") +
  xlab("Grade") + 
  ylab("PIAT Mathematics") + 
  scale_x_continuous(limits=c(2,8), breaks=seq(2,8,by=1)) +
  scale_y_continuous(limits=c(0,90), breaks=seq(0,90,by=10))

## Warning: Removed 4310 rows containing missing values (geom_point).

## Warning: Removed 2787 row(s) containing missing values (geom_path).

#Plotting intraindividual change READING
ggplot(data = data_long, aes(x = grade, y = rec, group = id)) +
  geom_point(color="red") + 
  geom_line(color="red") +
  xlab("Grade") + 
  ylab("PIAT Reading Recognition") +
  scale_x_continuous(limits=c(2,8), breaks=seq(2,8,by=1)) +
  scale_y_continuous(limits=c(0,90), breaks=seq(0,90,by=10))

## Warning: Removed 4317 rows containing missing values (geom_point).

## Warning: Removed 2804 row(s) containing missing values (geom_path).

#plotting intraindividaul change in the bivariate space
#Plotting vector field with .25 unit time change 
ggplot(data = data_long, aes(x = rec, y = math, group = id)) +
  geom_point(alpha=.1) + 
  geom_line(alpha=.3, arrow = arrow(length = unit(0.1, "cm"))) +
  xlab("Reading Recognition") + 
  ylab("Mathematics") +
  scale_x_continuous(limits=c(10,90), breaks=seq(10,90,by=10)) +
  scale_y_continuous(limits=c(10,90), breaks=seq(10,90,by=10))

## Warning: Removed 4317 rows containing missing values (geom_point).

## Warning: Removed 4317 row(s) containing missing values (geom_path).

This last plot shows individual trajectories in the bivariate space. Our general research goal is to describe how movement in the horizontal direction is coupled with movement in the vertical direction.

2.1 Model Specification

We first specify the full coupling model. The model invokes the latent variables and builds the difference scores, incorporates the constant growth factor and the proportional change. Lines invoking each set of latent variables and arrows is indicated in the model specification.

The model diagram follows this kind of setup …

Model Specification

bdcm_lavaan <- ' #opening quote
#MATHEMATICS
  #latent true scores (loadings = 1)
    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 (initial free, others = 0)  
    lm1 ~ 1
    lm2 ~ 0*1
    lm3 ~ 0*1
    lm4 ~ 0*1
    lm5 ~ 0*1
    lm6 ~ 0*1
    lm7 ~ 0*1

  #latent true score variances (initial free, others = 0)
    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 0)
    math2 ~ 0*1
    math3 ~ 0*1
    math4 ~ 0*1
    math5 ~ 0*1
    math6 ~ 0*1
    math7 ~ 0*1
    math8 ~ 0*1

  #observed residual variances (constrained to equality)
    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 (fixed = 1)
    lm2 ~ 1*lm1
    lm3 ~ 1*lm2
    lm4 ~ 1*lm3
    lm5 ~ 1*lm4
    lm6 ~ 1*lm5
    lm7 ~ 1*lm6

  #latent change scores (fixed = 1)
    dm2 =~ 1*lm2
    dm3 =~ 1*lm3
    dm4 =~ 1*lm4
    dm5 =~ 1*lm5
    dm6 =~ 1*lm6
    dm7 =~ 1*lm7

  #latent change score means (constrained to 0)  
    dm2 ~ 0*1
    dm3 ~ 0*1
    dm4 ~ 0*1
    dm5 ~ 0*1
    dm6 ~ 0*1
    dm7 ~ 0*1

  #latent change score variances (constrained to 0)  
    dm2 ~~ 0*dm2
    dm3 ~~ 0*dm3
    dm4 ~~ 0*dm4
    dm5 ~~ 0*dm5
    dm6 ~~ 0*dm6
    dm7 ~~ 0*dm7

  #constant change factor (loadings = 1)
    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 (constrained equal)
    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 RECOGNITION
  #latent true scores (loadings = 1)
    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 (initial free, others = 0)    
    lr1 ~ 1
    lr2 ~ 0*1
    lr3 ~ 0*1
    lr4 ~ 0*1
    lr5 ~ 0*1
    lr6 ~ 0*1
    lr7 ~ 0*1

  #latent true score variances (initial free, others = 0)
    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 = 0)
    rec2 ~ 0*1
    rec3 ~ 0*1
    rec4 ~ 0*1
    rec5 ~ 0*1
    rec6 ~ 0*1
    rec7 ~ 0*1
    rec8 ~ 0*1

  #observed residual variances (constrained to eqaulity)
    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 (fixed = 1)
    lr2 ~ 1*lr1
    lr3 ~ 1*lr2
    lr4 ~ 1*lr3
    lr5 ~ 1*lr4
    lr6 ~ 1*lr5
    lr7 ~ 1*lr6

  #latent change scores (ficed = 1)
    dr2 =~ 1*lr2
    dr3 =~ 1*lr3
    dr4 =~ 1*lr4
    dr5 =~ 1*lr5
    dr6 =~ 1*lr6
    dr7 =~ 1*lr7

  #latent change score means (fixed = 0)  
    dr2 ~ 0*1
    dr3 ~ 0*1
    dr4 ~ 0*1
    dr5 ~ 0*1
    dr6 ~ 0*1
    dr7 ~ 0*1

  #latent change score variances (fixed = 0)
    dr2 ~~ 0*dr2
    dr3 ~~ 0*dr3
    dr4 ~~ 0*dr4
    dr5 ~~ 0*dr5
    dr6 ~~ 0*dr6
    dr7 ~~ 0*dr7

  #constant change factor (fixed = 1)
    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 (constrained to equality)
    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
    lm1 ~~ j2
    lr1 ~~ g2
    j2 ~~ 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 PARMETERS
  
  #math to changes in 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 changes in 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
' #closing quote

2.2 Model Estimation and Interpretation

We fit the model using lavaan.

#Model fitting
dcs_fullcoupling <- lavaan.fit <- lavaan(bdcm_lavaan,
                data = nlsy_data, #note that fitting uses wide 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, cluster = cluster, : lavaan WARNING: some cases are empty and will be ignored:
##   741

## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan
## WARNING: due to missing values, some pairwise combinations have less than 10%
## coverage

#Model summary
summary(dcs_fullcoupling, fit.measures=TRUE)

## lavaan 0.6-8 ended normally after 178 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        58
##   Number of equality constraints                    37
##                                                       
##                                                   Used       Total
##   Number of observations                           932         933
##   Number of missing patterns                        66            
##                                                                   
## Model Test User Model:
##                                                       
##   Test statistic                               166.098
##   Degrees of freedom                                98
##   P-value (Chi-square)                           0.000
## 
## Model Test Baseline Model:
## 
##   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
##                                                       
##   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 - lower         0.020
##   90 Percent confidence interval - upper         0.034
##   P-value RMSEA <= 0.05                          1.000
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.111
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Observed
##   Observed information based on                Hessian
## 
## 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
##     j2                6.163    3.103    1.986    0.047
##   g2 ~~                                               
##     lr1               9.791    3.083    3.176    0.001
##     j2                0.681    2.552    0.267    0.790
##  .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.843    0.000
##    .rec3  (sgm2_s)   33.217    1.763   18.843    0.000
##    .rec4  (sgm2_s)   33.217    1.763   18.843    0.000
##    .rec5  (sgm2_s)   33.217    1.763   18.843    0.000
##    .rec6  (sgm2_s)   33.217    1.763   18.843    0.000
##    .rec7  (sgm2_s)   33.217    1.763   18.843    0.000
##    .rec8  (sgm2_s)   33.217    1.763   18.843    0.000
##    .dr2               0.000                           
##    .dr3               0.000                           
##    .dr4               0.000                           
##    .dr5               0.000                           
##    .dr6               0.000                           
##    .dr7               0.000                           
##     j2               18.384    6.806    2.701    0.007

#Model Diagram 
# semPaths(dcs_fullcoupling, 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, better to map to the diagram from the book. 

The change equations based on the ouput could be written as

Change in Math:

\[\Delta_{math} = 15.09 - 0.293(math_{t-1}) + 0.053(rec_{t-1})\] Change in Reading Recognition:

\[\Delta_{rec} = 10.89 - 0.495(rec_{t-1}) + 0.391(math_{t-1})\] However, need to be careful in interpretation as not all of the parameters are significantly different than zero.

The different coupling parameters (math to change in reading, reading to change in math) can each be constrained to zero in order to test the relative fit of models with unidirectional coupling and no coupling in order to reject some models and retain others.

2.3 Predicted scores

#obtaining predicted scores
nlsy_predicted <- cbind(nlsy_data$id, as.data.frame(lavPredict(dcs_fullcoupling, type = "yhat")))
names(nlsy_predicted)[1] <- "id"
#looking at data
head(nlsy_predicted)

#reshaping wide to long
predicted_long <- reshape(data=nlsy_predicted,
                    varying = c('math2', 'math3', 'math4', 'math5', 'math6', 'math7', 'math8',
                           'rec2', 'rec3', 'rec4', 'rec5', 'rec6', 'rec7', 'rec8'),
                    timevar=c("grade"), 
                    idvar=c("id"),
                    direction="long", sep="")

#sorting for easy viewing
#reorder by id and day
predicted_long <- predicted_long[order(predicted_long$id,predicted_long$grade), ]

#looking at the long data
head(predicted_long, 14)

#Plotting intraindividual change MATH
ggplot(data = predicted_long, aes(x = grade, y = math, group = id)) +
  #geom_point(color="blue") + 
  geom_line(color="blue") +
  xlab("Grade") + 
  ylab("Predicted PIAT Mathematics") + 
  scale_x_continuous(limits=c(2,8), breaks=seq(2,8,by=1)) +
  scale_y_continuous(limits=c(0,90), breaks=seq(0,90,by=10))

## Warning: Removed 7 row(s) containing missing values (geom_path).

#Plotting intraindividual change MATH
ggplot(data = predicted_long, aes(x = grade, y = rec, group = id)) +
  #geom_point(color="red") + 
  geom_line(color="red") +
  xlab("Grade") + 
  ylab("Predicted PIAT Reading Recognition") +
  scale_x_continuous(limits=c(2,8), breaks=seq(2,8,by=1)) +
  scale_y_continuous(limits=c(0,90), breaks=seq(0,90,by=10))

## Warning: Removed 12 row(s) containing missing values (geom_path).

These figures (similar to Figure 17.2 in the book) show the nonlinear exponential-type shape that is captured by the dual change score model.

2.4 Vector field plot of the Bivariate Dynamics

It is informative to examine how the coupling manifests in the bivariate space. To portray that, we create a vector field plot that indicates the probable movement in the bivariate space from any given starting location. We take a shortcut to the final plot using the RAMpath library, which has special functions for fitting the bivariate coupling model.

#creating a grid of starting values 
df <- expand.grid(math=seq(10, 90, 5), rec=seq(10, 90, 5))

#calculating change scores for each starting value 
#changes in math based on output from model
df$dm <- with(df, 15.09 - 0.293*math + 0.053*rec)
#changes in reading based on output from model
df$dr <- with(df, 10.89 - 0.495*rec + 0.391*math)


#Plotting vector field with .25 unit time change 
ggplot(data = df, aes(x = rec, y = math)) +
  geom_point(data=data_long, aes(x=rec, y=math), alpha=.1) + 
  stat_ellipse(data=data_long, aes(x=rec, y=math)) +
  geom_segment(aes(x = rec, y = math, xend = rec+.25*dr, yend = math+.25*dm), arrow = arrow(length = unit(0.1, "cm"))) +
  xlab("Reading Recognition") + 
  ylab("Mathematics") +
  scale_x_continuous(limits=c(10,90), breaks=seq(10,90,by=10)) +
  scale_y_continuous(limits=c(10,90), breaks=seq(10,90,by=10))

## Warning: Removed 4317 rows containing non-finite values (stat_ellipse).

## Warning: Removed 4317 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_segment).

Now we’ve got a plot of the change dynamics!

The arrows in the vector field plot (similar to Figure 17.3 in the book) show the expected changes in the bivariate space for a .25 unit of time. Note that the data are located mostly within the eillipse (95%). Outside that area the directional vectors are interpolated based on what was observed and modeled in the actual data (i.e., within the ellipse). So, please be careful not to over-interpret what is happening outside the ellipse.