A better cross-lagged panel model, from Hamaker et al. (2015)

October 20, 2017

This walk-through explains, briefly, why and how to run a RI-CLPM in R.

Critique of cross-lagged pannel models

This post summarizes critiques of the traditional cross-lagged panel model (CLPM), and an improved model by Hamaker, Kuiper, and Grasman (2015).

The primary point Hamaker and colleagues make regarding the CLPM is that it assumes that there are “no trait-like individual differences that endure.” That is, looking at the structure of a CLPM it is clear that individual-level stability must be accounted for entirely by the auto-regressive path between waves. As they put it, it imposes an assumption that there is no between-subject variance of time-invariant, trait-like stability, but only temporal stability, wave to wave, of subjects around the mean score for any particular wave.

RI-CLPM

A key insight of the paper is that “we need to separate the within-person level from the between-person level” (p. 104). The model they propose, the Random Intercept CLPM (RI-CLPM) separates each person’s score on a variable at each wave into the group mean for that wave ($\mu_{t}, \pi_{t}$), an individual’s stable score over all waves (the random intercept; $\kappa_{i}, \omega_{i}$) and then an individual level deviation at each wave from the score expected by adding the group-wave-mean and individual trait ($p_{it}, q_{it}$).

The model looks like this:

RI-CLPM Diagram

Effectively, now, the paths $\alpha_{t}$ (or $\delta_{t}$) between $p_{it}$ (or $q_{it}$) and $p_{i(t+1)}$ (or $q_{i(t+1)}$) no longer capture rank-order stability of individuals, but rather a within-person carry-over effect.

If it is positive, it implies that occasions on which a person scored above his or her expected score are likely to be followed by occasions on which he or she still scores above the expected score again, and conversely. (p. 104)

(Hamaker, Kuiper, & Grasman, 2015)

More importantly, since $\kappa$ and $\omega$ separate out individual-level stability, the cross-lagged paths $\beta_{t}$ and $\gamma_{t}$ are now straightforward to interpret as the within person effect of one variable on the subsequent measurement of a second variable. This interpretive boost is allowed now because, for example, $\beta_{t}$ is the estimate of the additional explanatory power of deviations from trait-stable levels on variable $y_{t}$ on the deviations of the observed variable $x_{t+1}$ from the group mean and individual trait ($\mu_{t+1} + \kappa_{i}$) after accounting for the expected within-person carry-over effect, $\alpha_{t}$.

See the paper (Figure 2) for a demonstration of how terribly traditional CLPM performs when you have a data generating process that matches the RI-CLPM – that is, when you have stable individual differences.

Implemmenting the RI-CLPM in R

First, make sure we can load the necessary packages. If you haven’t installed them, I’ve included (and commented out) the lines that will allow you to do that. The corresponding git repository is here.

#if you need to install anything, uncomment the below install lines for now
#install.packages('lavaan')
#install.packages('tidyverse')
require(lavaan)
require(tidyverse)

Some data

Now, we need some data. I’m using a data set presented on at a methods symposium at SRCD in 1997. Supporting documentation can be found in this pdf. Data and code for importing it was helpfully provided by Sanjay Srivastava.

The variables we’re considering are a measure of antisocial behavior (anti) and reading recognition (read). See the docs for descriptions of the other variables. And for the purpose of the model fitting below, x <- anit and y <- read. Following are some descriptions of the raw data:

antiread <- read.table("srcddata.dat",
                       na.strings = c("999.00"),
                       col.names = c("anti1", "anti2", "anti3", "anti4", 
                                     "read1", "read2", "read3", "read4",
                                     "gen", "momage", "kidage", "homecog", 
                                     "homeemo", "id")
) %>%
  rename(x1 = anti1, x2 = anti2, x3 = anti3, x4 = anti4,
         y1 = read1, y2 = read2, y3 = read3, y4 = read4) %>%
  select(matches('[xy][1-4]'))

knitr::kable(summary(antiread), format = 'markdown')
  x1 x2 x3 x4 y1 y2 y3 y4
  Min. :0.000 Min. : 0.000 Min. : 0.000 Min. : 0.000 Min. :0.100 Min. :1.600 Min. :2.200 Min. :2.500
  1st Qu.:0.000 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.:1.800 1st Qu.:3.300 1st Qu.:4.200 1st Qu.:4.925
  Median :1.000 Median : 1.500 Median : 1.000 Median : 1.500 Median :2.300 Median :4.100 Median :5.000 Median :5.800
  Mean :1.662 Mean : 2.027 Mean : 1.828 Mean : 2.061 Mean :2.523 Mean :4.076 Mean :5.005 Mean :5.774
  3rd Qu.:3.000 3rd Qu.: 3.000 3rd Qu.: 3.000 3rd Qu.: 3.000 3rd Qu.:3.000 3rd Qu.:4.900 3rd Qu.:5.800 3rd Qu.:6.675
  Max. :9.000 Max. :10.000 Max. :10.000 Max. :10.000 Max. :7.200 Max. :8.200 Max. :8.400 Max. :8.300
  NA NA’s :31 NA’s :108 NA’s :111 NA NA’s :30 NA’s :130 NA’s :135
antiread %>%
  select(-x4,-y4) %>%
  mutate(pid = 1:n()) %>%
  gather(key, value, -pid) %>%
  extract(col = key, into = c('var', 'wave'), regex = '(\\w)(\\d)') %>%
  ggplot(aes(x = value)) +
  geom_density(alpha = 1) + 
  facet_grid(wave~var, scales = 'free') + 
  theme_classic()
## Warning: Removed 299 rows containing non-finite values (stat_density).

center

antireadLong <- antiread %>%
  select(-x4,-y4) %>%
  mutate(pid = 1:n()) %>%
  gather(key, value, -pid) %>%
  extract(col = key, into = c('var', 'wave'), regex = '(\\w)(\\d)')

antireadLong %>%
  ggplot(aes(x = wave, y = value, color = var, group = var)) +
  geom_point(position = position_jitter(w = .2), alpha = .1) +
  geom_line(stat = 'identity', aes(group = interaction(var, pid)), alpha = .04) + 
  geom_line(stat = 'smooth', method = 'lm', size = 1) + 
  theme_classic()
## Warning: Removed 299 rows containing non-finite values (stat_smooth).
## Warning: Removed 299 rows containing missing values (geom_point).
## Warning: Removed 271 rows containing missing values (geom_path).

center

Fitting a RI-CLPM

In the below lavaan code, I’ll be using the notation from the diagram. I am explicitly specifying everything in the diagram, which is why in the call to lavaan I set a bunch of auto options to false. This is because often lavaan will try to automatically estimate things that you don’t usually write out but often want estimated, like residuals. Because this model is unorthodox, I want to be as explicit as possible.

The lavaan code below uses syntax that can be found in their help docs for the basic stuff as well as the more advanced labeling and constraining.

riclpmModel <- 
'
#Note, the data contain x1-3 and y1-3
#Latent mean Structure with intercepts

kappa =~ 1*x1 + 1*x2 + 1*x3
omega =~ 1*y1 + 1*y2 + 1*y3

x1 ~ mu1*1 #intercepts
x2 ~ mu2*1
x3 ~ mu3*1
y1 ~ pi1*1
y2 ~ pi2*1
y3 ~ pi3*1

kappa ~~ kappa #variance
omega ~~ omega #variance
kappa ~~ omega #covariance

#laten vars for AR and cross-lagged effects
p1 =~ 1*x1 #each factor loading set to 1
p2 =~ 1*x2
p3 =~ 1*x3
q1 =~ 1*y1
q2 =~ 1*y2
q3 =~ 1*y3

#Later, we may constrain autoregression and cross-lagged
#effects to be the same across both lags.
p3 ~ alpha3*p2 + beta3*q2
p2 ~ alpha2*p1 + beta2*q1

q3 ~ delta3*q2 + gamma3*p2
q2 ~ delta2*q1 + gamma2*p1

p1 ~~ p1 #variance
p2 ~~ u2*p2
p3 ~~ u3*p3
q1 ~~ q1 #variance
q2 ~~ v2*q2
q3 ~~ v3*q3

p1 ~~ q1 #p1 and q1 covariance
p2 ~~ q2 #p2 and q2 covariance
p3 ~~ q3 #p2 and q2 covariance'

fit <- lavaan(riclpmModel, data = antiread,
              missing = 'ML', #for the missing data!
              int.ov.free = F,
              int.lv.free = F,
              auto.fix.first = F,
              auto.fix.single = F,
              auto.cov.lv.x = F,
              auto.cov.y = F,
              auto.var = F)
summary(fit, standardized = T)
## lavaan (0.5-23.1097) converged normally after  83 iterations
## 
##   Number of observations                           405
## 
##   Number of missing patterns                        14
## 
##   Estimator                                         ML
##   Minimum Function Test Statistic                3.213
##   Degrees of freedom                                 1
##   P-value (Chi-square)                           0.073
## 
## Parameter Estimates:
## 
##   Information                                 Observed
##   Standard Errors                             Standard
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   kappa =~                                                              
##     x1                1.000                               1.110    0.671
##     x2                1.000                               1.110    0.548
##     x3                1.000                               1.110    0.571
##   omega =~                                                              
##     y1                1.000                               0.734    0.794
##     y2                1.000                               0.734    0.679
##     y3                1.000                               0.734    0.622
##   p1 =~                                                                 
##     x1                1.000                               1.226    0.741
##   p2 =~                                                                 
##     x2                1.000                               1.692    0.836
##   p3 =~                                                                 
##     x3                1.000                               1.594    0.821
##   q1 =~                                                                 
##     y1                1.000                               0.562    0.608
##   q2 =~                                                                 
##     y2                1.000                               0.794    0.734
##   q3 =~                                                                 
##     y3                1.000                               0.924    0.783
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   p3 ~                                                                  
##     p2      (alp3)    0.353    0.078    4.549    0.000    0.374    0.374
##     q2      (bet3)   -0.274    0.154   -1.788    0.074   -0.137   -0.137
##   p2 ~                                                                  
##     p1      (alp2)    0.162    0.169    0.961    0.336    0.118    0.118
##     q1      (bet2)   -0.090    0.457   -0.196    0.845   -0.030   -0.030
##   q3 ~                                                                  
##     q2      (dlt3)    0.738    0.073   10.048    0.000    0.634    0.634
##     p2      (gmm3)   -0.008    0.034   -0.241    0.810   -0.015   -0.015
##   q2 ~                                                                  
##     q1      (dlt2)    0.374    0.350    1.068    0.286    0.265    0.265
##     p1      (gmm2)   -0.058    0.079   -0.733    0.464   -0.089   -0.089
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   kappa ~~                                                              
##     omega            -0.118    0.157   -0.754    0.451   -0.145   -0.145
##   p1 ~~                                                                 
##     q1                0.017    0.160    0.107    0.915    0.025    0.025
##  .p2 ~~                                                                 
##    .q2               -0.117    0.114   -1.025    0.305   -0.091   -0.091
##  .p3 ~~                                                                 
##    .q3               -0.115    0.071   -1.624    0.104   -0.111   -0.111
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .x1       (mu1)    1.662    0.082   20.217    0.000    1.662    1.005
##    .x2       (mu2)    1.985    0.103   19.189    0.000    1.985    0.981
##    .x3       (mu3)    1.898    0.107   17.658    0.000    1.898    0.977
##    .y1       (pi1)    2.523    0.046   54.925    0.000    2.523    2.729
##    .y2       (pi2)    4.066    0.055   74.267    0.000    4.066    3.760
##    .y3       (pi3)    5.023    0.064   78.328    0.000    5.023    4.256
##     kappa             0.000                               0.000    0.000
##     omega             0.000                               0.000    0.000
##     p1                0.000                               0.000    0.000
##    .p2                0.000                               0.000    0.000
##    .p3                0.000                               0.000    0.000
##     q1                0.000                               0.000    0.000
##    .q2                0.000                               0.000    0.000
##    .q3                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     kappa             1.232    0.271    4.550    0.000    1.000    1.000
##     omega             0.539    0.177    3.042    0.002    1.000    1.000
##     p1                1.504    0.281    5.361    0.000    1.000    1.000
##    .p2        (u2)    2.821    0.316    8.929    0.000    0.985    0.985
##    .p3        (u3)    2.110    0.204   10.341    0.000    0.830    0.830
##     q1                0.316    0.172    1.832    0.067    1.000    1.000
##    .q2        (v2)    0.582    0.086    6.811    0.000    0.923    0.923
##    .q3        (v3)    0.509    0.046   11.125    0.000    0.596    0.596
##    .x1                0.000                               0.000    0.000
##    .x2                0.000                               0.000    0.000
##    .x3                0.000                               0.000    0.000
##    .y1                0.000                               0.000    0.000
##    .y2                0.000                               0.000    0.000
##    .y3                0.000                               0.000    0.000

Comparing fits

RI-CLPM v CLPM

Because the traditional CLPM is nested in the RI-CLPM, we can compare model fit. The correct reference distribution for this comparison is $\chi^2$ but, as Hamaker and colleagues say

However, we can make use of the fact that the regular chi-square difference test is conservative in this context, meaning that, if it is significant, we are certain that the correct (i.e., chi-bar-square difference) test will be significant too. (p. 105)

(Hamaker, Kuiper, & Grasman, 2015)

We estimate the traditional CLPM by setting the variance and covariance of $\kappa$ and $\omega$ to 0.

clpmModel <- #yes, "Model" is redundant
'
#Note, the data contain x1-3 and y1-3
#Latent mean Structure with intercepts

kappa =~ 1*x1 + 1*x2 + 1*x3
omega =~ 1*y1 + 1*y2 + 1*y3

x1 ~ mu1*1 #intercepts
x2 ~ mu2*1
x3 ~ mu3*1
y1 ~ pi1*1
y2 ~ pi2*1
y3 ~ pi3*1

kappa ~~ 0*kappa #variance nope
omega ~~ 0*omega #variance nope
kappa ~~ 0*omega #covariance not even

#laten vars for AR and cross-lagged effects
p1 =~ 1*x1 #each factor loading set to 1
p2 =~ 1*x2
p3 =~ 1*x3
q1 =~ 1*y1
q2 =~ 1*y2
q3 =~ 1*y3

p3 ~ alpha3*p2 + beta3*q2
p2 ~ alpha2*p1 + beta2*q1

q3 ~ delta3*q2 + gamma3*p2
q2 ~ delta2*q1 + gamma2*p1

p1 ~~ p1 #variance
p2 ~~ u2*p2
p3 ~~ u3*p3
q1 ~~ q1 #variance
q2 ~~ v2*q2
q3 ~~ v3*q3

p1 ~~ q1 #p1 and q1 covariance
p2 ~~ q2 #p2 and q2 covariance
p3 ~~ q3 #p2 and q2 covariance'

fitCLPM <- lavaan(clpmModel, data = antiread,
              missing = 'ML', #for the missing data!
              int.ov.free = F,
              int.lv.free = F,
              auto.fix.first = F,
              auto.fix.single = F,
              auto.cov.lv.x = F,
              auto.cov.y = F,
              auto.var = F)

anova(fit, fitCLPM)
## Chi Square Difference Test
## 
##         Df    AIC    BIC   Chisq Chisq diff Df diff Pr(>Chisq)    
## fit      1 6814.1 6918.2  3.2127                                  
## fitCLPM  4 6825.6 6917.7 20.6779     17.465       3  0.0005669 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The CLPM fits much worse, with no model comparison statistic favoring the CLPM (the BIC advantage of 1 point is negligible). We can print out the standardized estimates to compare to the unconstrained RI-CLPM above.

summary(fitCLPM, standardize = T)
## lavaan (0.5-23.1097) converged normally after  47 iterations
## 
##   Number of observations                           405
## 
##   Number of missing patterns                        14
## 
##   Estimator                                         ML
##   Minimum Function Test Statistic               20.678
##   Degrees of freedom                                 4
##   P-value (Chi-square)                           0.000
## 
## Parameter Estimates:
## 
##   Information                                 Observed
##   Standard Errors                             Standard
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   kappa =~                                                              
##     x1                1.000                               0.000    0.000
##     x2                1.000                               0.000    0.000
##     x3                1.000                               0.000    0.000
##   omega =~                                                              
##     y1                1.000                               0.000    0.000
##     y2                1.000                               0.000    0.000
##     y3                1.000                               0.000    0.000
##   p1 =~                                                                 
##     x1                1.000                               1.656    1.000
##   p2 =~                                                                 
##     x2                1.000                               2.024    1.000
##   p3 =~                                                                 
##     x3                1.000                               1.955    1.000
##   q1 =~                                                                 
##     y1                1.000                               0.924    1.000
##   q2 =~                                                                 
##     y2                1.000                               1.082    1.000
##   q3 =~                                                                 
##     y3                1.000                               1.174    1.000
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   p3 ~                                                                  
##     p2      (alp3)    0.559    0.050   11.150    0.000    0.579    0.579
##     q2      (bet3)   -0.145    0.085   -1.698    0.090   -0.080   -0.080
##   p2 ~                                                                  
##     p1      (alp2)    0.538    0.056    9.679    0.000    0.440    0.440
##     q1      (bet2)   -0.093    0.101   -0.920    0.357   -0.042   -0.042
##   q3 ~                                                                  
##     q2      (dlt3)    0.849    0.042   20.402    0.000    0.783    0.783
##     p2      (gmm3)   -0.016    0.024   -0.677    0.498   -0.028   -0.028
##   q2 ~                                                                  
##     q1      (dlt2)    0.763    0.045   16.903    0.000    0.651    0.651
##     p1      (gmm2)   -0.047    0.025   -1.871    0.061   -0.072   -0.072
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   kappa ~~                                                              
##     omega             0.000                                 NaN      NaN
##   p1 ~~                                                                 
##     q1               -0.107    0.076   -1.403    0.161   -0.070   -0.070
##  .p2 ~~                                                                 
##    .q2               -0.085    0.077   -1.098    0.272   -0.058   -0.058
##  .p3 ~~                                                                 
##    .q3               -0.121    0.071   -1.710    0.087   -0.106   -0.106
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .x1       (mu1)    1.662    0.082   20.194    0.000    1.662    1.003
##    .x2       (mu2)    1.983    0.103   19.189    0.000    1.983    0.980
##    .x3       (mu3)    1.902    0.109   17.503    0.000    1.902    0.973
##    .y1       (pi1)    2.523    0.046   54.956    0.000    2.523    2.731
##    .y2       (pi2)    4.066    0.055   74.277    0.000    4.066    3.759
##    .y3       (pi3)    5.023    0.064   78.548    0.000    5.023    4.279
##     kappa             0.000                                 NaN      NaN
##     omega             0.000                                 NaN      NaN
##     p1                0.000                               0.000    0.000
##    .p2                0.000                               0.000    0.000
##    .p3                0.000                               0.000    0.000
##     q1                0.000                               0.000    0.000
##    .q2                0.000                               0.000    0.000
##    .q3                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     kappa             0.000                                 NaN      NaN
##     omega             0.000                                 NaN      NaN
##     p1                2.742    0.193   14.230    0.000    1.000    1.000
##    .p2        (u2)    3.284    0.239   13.726    0.000    0.802    0.802
##    .p3        (u3)    2.474    0.206   12.035    0.000    0.648    0.648
##     q1                0.853    0.060   14.230    0.000    1.000    1.000
##    .q2        (v2)    0.660    0.048   13.698    0.000    0.564    0.564
##    .q3        (v3)    0.525    0.046   11.509    0.000    0.381    0.381
##    .x1                0.000                               0.000    0.000
##    .x2                0.000                               0.000    0.000
##    .x3                0.000                               0.000    0.000
##    .y1                0.000                               0.000    0.000
##    .y2                0.000                               0.000    0.000
##    .y3                0.000                               0.000    0.000

Adding constraints to RI-CLPM

For parsimony, I usually try to constraint my autoregressive and cross-lagged paths to be the same across intervals. Oh, and residuals too. I’ll do this in the following code and then check the fit against the unconstrained model. To do this, all I have to do is make sure the paths have the same name, like alpha instead of alpha2 and alpha3.

riclpmModelConstrainedARCL <- 
'
#Note, the data contain x1-3 and y1-3
#Latent mean Structure with intercepts

kappa =~ 1*x1 + 1*x2 + 1*x3
omega =~ 1*y1 + 1*y2 + 1*y3

x1 ~ mu1*1 #intercepts
x2 ~ mu2*1
x3 ~ mu3*1
y1 ~ pi1*1
y2 ~ pi2*1
y3 ~ pi3*1

kappa ~~ kappa #variance
omega ~~ omega #variance
kappa ~~ omega #covariance

#laten vars for AR and cross-lagged effects
p1 =~ 1*x1 #each factor loading set to 1
p2 =~ 1*x2
p3 =~ 1*x3
q1 =~ 1*y1
q2 =~ 1*y2
q3 =~ 1*y3

#constrain autoregression and cross lagged effects to be the same across both lags.
p3 ~ alpha*p2 + beta*q2
p2 ~ alpha*p1 + beta*q1

q3 ~ delta*q2 + gamma*p2
q2 ~ delta*q1 + gamma*p1

p1 ~~ p1 #variance
p2 ~~ u*p2
p3 ~~ u*p3
q1 ~~ q1 #variance
q2 ~~ v*q2
q3 ~~ v*q3

p1 ~~ q1 #p1 and q1 covariance
p2 ~~ uv*q2 #p2 and q2 covariance should also be constrained to be the same as
p3 ~~ uv*q3 #p3 and q3 covariance'

fitConstrainedARCL <- lavaan(riclpmModelConstrainedARCL, data = antiread,
              missing = 'ML', #for the missing data!
              int.ov.free = F,
              int.lv.free = F,
              auto.fix.first = F,
              auto.fix.single = F,
              auto.cov.lv.x = F,
              auto.cov.y = F,
              auto.var = F)

anova(fit, fitConstrainedARCL)
## Chi Square Difference Test
## 
##                    Df    AIC    BIC   Chisq Chisq diff Df diff Pr(>Chisq)
## fit                 1 6814.1 6918.2  3.2127                              
## fitConstrainedARCL  8 6820.2 6896.2 23.2640     20.051       7    0.00546
##                      
## fit                  
## fitConstrainedARCL **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Well, according to AIC and the $\chi^2$ test, the constrained model fits worse. But BIC loves the constrained model because it hates parameters. Interpretive ease hates parameters too (most of the time), so let’s look at the summary for our simplified model.

summary(fitConstrainedARCL, standardized = T)
## lavaan (0.5-23.1097) converged normally after  78 iterations
## 
##   Number of observations                           405
## 
##   Number of missing patterns                        14
## 
##   Estimator                                         ML
##   Minimum Function Test Statistic               23.264
##   Degrees of freedom                                 8
##   P-value (Chi-square)                           0.003
## 
## Parameter Estimates:
## 
##   Information                                 Observed
##   Standard Errors                             Standard
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   kappa =~                                                              
##     x1                1.000                               1.025    0.616
##     x2                1.000                               1.025    0.525
##     x3                1.000                               1.025    0.516
##   omega =~                                                              
##     y1                1.000                               0.566    0.612
##     y2                1.000                               0.566    0.519
##     y3                1.000                               0.566    0.484
##   p1 =~                                                                 
##     x1                1.000                               1.311    0.788
##   p2 =~                                                                 
##     x2                1.000                               1.662    0.851
##   p3 =~                                                                 
##     x3                1.000                               1.702    0.857
##   q1 =~                                                                 
##     y1                1.000                               0.732    0.791
##   q2 =~                                                                 
##     y2                1.000                               0.932    0.855
##   q3 =~                                                                 
##     y3                1.000                               1.024    0.875
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   p3 ~                                                                  
##     p2      (alph)    0.306    0.092    3.331    0.001    0.299    0.299
##     q2      (beta)   -0.212    0.148   -1.433    0.152   -0.116   -0.116
##   p2 ~                                                                  
##     p1      (alph)    0.306    0.092    3.331    0.001    0.241    0.241
##     q1      (beta)   -0.212    0.148   -1.433    0.152   -0.093   -0.093
##   q3 ~                                                                  
##     q2      (delt)    0.713    0.076    9.391    0.000    0.648    0.648
##     p2      (gamm)   -0.039    0.031   -1.281    0.200   -0.063   -0.063
##   q2 ~                                                                  
##     q1      (delt)    0.713    0.076    9.391    0.000    0.560    0.560
##     p1      (gamm)   -0.039    0.031   -1.281    0.200   -0.055   -0.055
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   kappa ~~                                                              
##     omega            -0.046    0.143   -0.321    0.748   -0.079   -0.079
##   p1 ~~                                                                 
##     q1               -0.056    0.146   -0.387    0.699   -0.059   -0.059
##  .p2 ~~                                                                 
##    .q2        (uv)   -0.116    0.060   -1.949    0.051   -0.094   -0.094
##  .p3 ~~                                                                 
##    .q3        (uv)   -0.116    0.060   -1.949    0.051   -0.094   -0.094
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .x1       (mu1)    1.662    0.083   20.095    0.000    1.662    0.999
##    .x2       (mu2)    1.990    0.100   19.944    0.000    1.990    1.019
##    .x3       (mu3)    1.890    0.111   16.976    0.000    1.890    0.951
##    .y1       (pi1)    2.523    0.046   54.891    0.000    2.523    2.728
##    .y2       (pi2)    4.067    0.055   73.785    0.000    4.067    3.731
##    .y3       (pi3)    5.018    0.064   78.033    0.000    5.018    4.290
##     kappa             0.000                               0.000    0.000
##     omega             0.000                               0.000    0.000
##     p1                0.000                               0.000    0.000
##    .p2                0.000                               0.000    0.000
##    .p3                0.000                               0.000    0.000
##     q1                0.000                               0.000    0.000
##    .q2                0.000                               0.000    0.000
##    .q3                0.000                               0.000    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     kappa             1.052    0.253    4.161    0.000    1.000    1.000
##     omega             0.320    0.167    1.915    0.055    1.000    1.000
##     p1                1.718    0.254    6.764    0.000    1.000    1.000
##    .p2         (u)    2.569    0.199   12.895    0.000    0.930    0.930
##    .p3         (u)    2.569    0.199   12.895    0.000    0.887    0.887
##     q1                0.535    0.166    3.225    0.001    1.000    1.000
##    .q2         (v)    0.590    0.036   16.503    0.000    0.680    0.680
##    .q3         (v)    0.590    0.036   16.503    0.000    0.563    0.563
##    .x1                0.000                               0.000    0.000
##    .x2                0.000                               0.000    0.000
##    .x3                0.000                               0.000    0.000
##    .y1                0.000                               0.000    0.000
##    .y2                0.000                               0.000    0.000
##    .y3                0.000                               0.000    0.000

Plotting model fit

Now, we can plot the model-fitted values. To examine how the model relates to the data, we’ll follow the principle of the model which is to partition the within versus between subject variance. You can reinforce the corresponding intuition by looking back at the path diagram: keep in mind that every observed value will be exactly equal to the wave mean, the individual’s latent intercept, and the per-wave latent residual (that is, p and q). So first, we’ll plot the individual variation around the wave-wise means (the stable, between subject individual differences captured by $\kappa$ and $\omega$), along with the observed values. You can see that there is a lot of distance between the lines (that is, the expected values based on the random intercept and wave-wise means) and the observed values. It is that deviation that the within-subject portion of the model (the cross-lagged part) is attempting to explain.

#get the model-expected means
means <- fitted(fitConstrainedARCL)$mean
meansDF <- data.frame(mean = means, key = names(means)) %>%
  extract(col = key, into = c('var', 'wave'), regex = '(\\w)(\\d)')

#plot the model-expected random intercepts
predict(fitConstrainedARCL) %>%
  as.data.frame %>%
  mutate(pid = 1:n()) %>%
  gather(key, latentvalue, -pid, -kappa, -omega) %>%
  extract(col = key, into = c('latentvar', 'wave'), regex = '(\\w)(\\d)') %>%
  mutate(var = c(p = 'x', q = 'y')[latentvar]) %>%
  left_join(meansDF) %>% #those means from above
  left_join(antireadLong, by = c('pid', 'wave', 'var')) %>% #the raw data
  mutate(expectedLine = ifelse(var == 'x', kappa, omega) + mean,
         wave = as.numeric(wave)) %>%
  rowwise() %>%
  ggplot(aes(x = wave, y = expectedLine, color = var, group = var)) +
  geom_point(aes(x = wave, y = value, group = interaction(var, pid)), alpha = .1, position = position_jitter(w = .2, h = 0)) +
  geom_line(aes(y = expectedLine, group = interaction(var, pid)), stat = 'identity', alpha = .1) + 
  geom_line(aes(y = mean), stat = 'identity', alpha = 1, size = 1, color = 'black') + 
  facet_wrap(~var, ncol = 2) + 
  theme_classic()

center

We can also look at the correlations between the latent residuals, which are essentially the parts of the observed scores that are not accounted for by the stable individual differences in the above graph.

That is, the autoregressive parameters $\alpha^{*}_{t}$ and $\delta^{*}_{t}$ do not represent
the stability of the rank order of individuals from one occasion to
the next, but rather the amount of within-person carry-over effect
(cf., Hamaker, 2012; Kuppens, Allen, & Sheeber, 2010; Suls,
Green, & Hillis, 1998): If it is positive, it implies that occasions on
which a person scored above his or her expected score are likely to
be followed by occasions on which he or she still scores above the
expected score again, and vice versa. (p. 104-105)

(Hamaker, Kuiper, & Grasman, 2015)

So you may interpret the raw correlations in the graph below as the basis for the constrained model estimates of the path weights above. In the interest of checking against over-interpreting these relations, though, here is the authors’ footnote to the above statement:

One could also say these autoregressive parameters indicate the stability of the rank-order of individual deviations, but this is less appealing
from a substantive viewpoint. (p. 105)

(Hamaker, Kuiper, & Grasman, 2015)
library(GGally)
predict(fitConstrainedARCL) %>%
  as.data.frame %>%
  select(-kappa, -omega) %>%
  ggpairs(lower = list(continuous = wrap(ggally_smooth, alpha = .5))) + 
  theme_classic()

center

Note how different, and one might say diminished, the relations in the above graph are versus the relations in the scatter-plot matrix of the raw data, below. The strength of this model seems to lie in its ability to keep one from being fooled into a within-person explanation of what are largely between-person relations.

antiread %>% 
  select(-x4, -y4) %>%
  ggpairs(lower = list(continuous = wrap(ggally_smooth, alpha = .5))) + 
  theme_classic()

center

  1. Hamaker, E. L., Kuiper, R. M., & Grasman, R. P. P. P. (2015). A critique of the cross-lagged panel model. Psychological Methods, 20(1), 102–116. https://doi.org/10.1037/a0038889
A better cross-lagged panel model, from Hamaker et al. (2015) - October 20, 2017 - John C. Flournoy