library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library (lavaan)
## This is lavaan 0.6-17
## lavaan is FREE software! Please report any bugs.
library (semPlot)
data <- read.csv("https://quantdev.ssri.psu.edu/sites/qdev/files/nlsy_data_subset.csv")
names(data)
## [1] "X" "id" "math2" "math3" "math4" "math5" "math6" "math7" "math8"
## [10] "rec2" "rec3" "rec4" "rec5" "rec6" "rec7" "rec8"
Class_Data <-(data[,-1])
names(Class_Data)
## [1] "id" "math2" "math3" "math4" "math5" "math6" "math7" "math8" "rec2"
## [10] "rec3" "rec4" "rec5" "rec6" "rec7" "rec8"
riclpm <- '
# Define intercept factors
i_math =~ 1*math2+1*math3+1*math4+1*math5
i_rec =~ 1*rec2+1*rec3+1*rec4+1*rec5
# Define single item latent variables
eta_math2 =~ 1*math2
eta_math3 =~ 1*math3
eta_math4 =~ 1*math4
eta_math5 =~ 1*math5
eta_rec2 =~ 1*rec2
eta_rec3 =~ 1*rec3
eta_rec4 =~ 1*rec4
eta_rec5 =~ 1*rec5
# Autoregressive effects
eta_math3 ~ a1*eta_math2
eta_math4 ~ a1*eta_math3
eta_math5 ~ a1*eta_math4
eta_rec3 ~ a2*eta_rec2
eta_rec4 ~ a2*eta_rec3
eta_rec5 ~ a2*eta_rec4
# Crosslagged effects
eta_rec3 ~ c1*eta_math2
eta_rec4 ~ c1*eta_math3
eta_rec5 ~ c1*eta_math4
eta_math3 ~ c2*eta_rec2
eta_math4 ~ c2*eta_rec3
eta_math5 ~ c2*eta_rec4
# Some further constraints on the variance structure
# 1. Set error variances of the observed variables to zero
math2 ~~ 0*math2
math3 ~~ 0*math3
math4 ~~ 0*math4
math5 ~~ 0*math5
rec2 ~~ 0*rec2
rec3 ~~ 0*rec3
rec4 ~~ 0*rec4
rec5 ~~ 0*rec5
# 2. Let lavaan estimate the variance of the latent variables
eta_math2 ~~ varmath2*eta_math2
eta_math3 ~~ varmath3*eta_math3
eta_math4 ~~ varmath4*eta_math4
eta_math5 ~~ varmath5*eta_math5
eta_rec2 ~~ varrec2*eta_rec2
eta_rec3 ~~ varrec3*eta_rec3
eta_rec4 ~~ varrec4*eta_rec4
eta_rec5 ~~ varrec5*eta_rec5
# 3. We also want estimates of the intercept factor variances and an
# estimate of their covariance
i_math ~~ varimath*i_math
i_rec ~~ varirec*i_rec
i_math ~~ covi*i_rec
# 4. We have to define that the covariance between the intercepts and
# the latents of the first time point are zero
eta_math2 ~~ 0*i_math
eta_rec2 ~~ 0*i_math
eta_math2 ~~ 0*i_rec
eta_rec2 ~~ 0*i_rec
# 5. Finally, we estimate the covariance between the latents of x and y
# of the first time point, the second time-point and so on. note that
# for the second to fourth time point the correlation is constrained to
# the same value
eta_math2 ~~ cov1*eta_rec2
eta_math3 ~~ e1*eta_rec3
eta_math4 ~~ e1*eta_rec4
eta_math5 ~~ e1*eta_rec5
# The model also contains a mean structure and we have to define some
# constraints for this part of the model. the assumption is that we
# only want estimates of the mean of the intercept factors. all other means
# are defined to be zero:
math2 ~ 0*1
math3 ~ 0*1
math4 ~ 0*1
math5 ~ 0*1
rec2 ~ 0*1
rec3 ~ 0*1
rec4 ~ 0*1
rec5 ~ 0*1
eta_math2 ~ 1
eta_math3 ~ 1
eta_math4 ~ 1
eta_math5 ~ 1
eta_rec2 ~ 1
eta_rec3 ~ 1
eta_rec4 ~ 1
eta_rec5 ~ 1
i_math ~ 1
i_rec ~ 1
## define correlations
cori := covi / (sqrt(varimath) * sqrt(varirec))
cor1 := cov1 / (sqrt(varmath2) * sqrt(varrec2))
cort2 := e1 / (sqrt(varmath3) * sqrt(varrec3))
cort3 := e1 / (sqrt(varmath4) * sqrt(varrec4))
cort4 := e1 / (sqrt(varmath5) * sqrt(varrec5))
'
riclpm_fit <- sem(riclpm, estimator= "MLR", data=Class_Data, mimic = "Mplus", missing = "FIML")
## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan WARNING: some cases are empty and will be ignored:
## 41 115 153 303 306 374 555 585 608 663 732 741 771 807 818 922
## 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; use lavInspect(fit, "coverage") to investigate.
## Warning in lav_mvnorm_missing_h1_estimate_moments(Y = X[[g]], wt = WT[[g]], : lavaan WARNING:
## Maximum number of iterations reached when computing the sample
## moments using EM; use the em.h1.iter.max= argument to increase the
## number of iterations
## Warning in lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats, : lavaan WARNING:
## The variance-covariance matrix of the estimated parameters (vcov)
## does not appear to be positive definite! The smallest eigenvalue
## (= 3.546990e-14) is close to zero. This may be a symptom that the
## model is not identified.
summary(riclpm_fit, fit.measures = TRUE, standardized=TRUE)
## lavaan 0.6.17 ended normally after 378 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 37
## Number of equality constraints 10
##
## Used Total
## Number of observations 917 933
## Number of missing patterns 17
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 48.768 46.534
## Degrees of freedom 17 17
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.048
## Yuan-Bentler correction (Mplus variant)
##
## Model Test Baseline Model:
##
## Test statistic 1381.532 1285.358
## Degrees of freedom 28 28
## P-value 0.000 0.000
## Scaling correction factor 1.075
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.977 0.977
## Tucker-Lewis Index (TLI) 0.961 0.961
##
## Robust Comparative Fit Index (CFI) 0.875
## Robust Tucker-Lewis Index (TLI) 0.795
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -10846.830 -10846.830
## Scaling correction factor 0.757
## for the MLR correction
## Loglikelihood unrestricted model (H1) -10822.446 -10822.446
## Scaling correction factor 1.042
## for the MLR correction
##
## Akaike (AIC) 21747.659 21747.659
## Bayesian (BIC) 21877.829 21877.829
## Sample-size adjusted Bayesian (SABIC) 21792.080 21792.080
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.045 0.044
## 90 Percent confidence interval - lower 0.031 0.029
## 90 Percent confidence interval - upper 0.060 0.058
## P-value H_0: RMSEA <= 0.050 0.682 0.748
## P-value H_0: RMSEA >= 0.080 0.000 0.000
##
## Robust RMSEA 0.194
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.342
## P-value H_0: Robust RMSEA <= 0.050 0.107
## P-value H_0: Robust RMSEA >= 0.080 0.862
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.069 0.069
##
## Parameter Estimates:
##
## Standard errors Sandwich
## Information bread Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## i_math =~
## math2 1.000 7.867 0.788
## math3 1.000 7.867 0.768
## math4 1.000 7.867 0.791
## math5 1.000 7.867 0.810
## i_rec =~
## rec2 1.000 8.628 0.852
## rec3 1.000 8.628 0.725
## rec4 1.000 8.628 0.697
## rec5 1.000 8.628 0.681
## eta_math2 =~
## math2 1.000 6.156 0.616
## eta_math3 =~
## math3 1.000 6.557 0.640
## eta_math4 =~
## math4 1.000 6.075 0.611
## eta_math5 =~
## math5 1.000 5.688 0.586
## eta_rec2 =~
## rec2 1.000 5.312 0.524
## eta_rec3 =~
## rec3 1.000 8.197 0.689
## eta_rec4 =~
## rec4 1.000 8.866 0.717
## eta_rec5 =~
## rec5 1.000 9.280 0.732
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## eta_math3 ~
## eta_math2 (a1) -0.213 0.123 -1.732 0.083 -0.200 -0.200
## eta_math4 ~
## eta_math3 (a1) -0.213 0.123 -1.732 0.083 -0.229 -0.229
## eta_math5 ~
## eta_math4 (a1) -0.213 0.123 -1.732 0.083 -0.227 -0.227
## eta_rec3 ~
## eta_rec2 (a2) 0.760 0.072 10.505 0.000 0.492 0.492
## eta_rec4 ~
## eta_rec3 (a2) 0.760 0.072 10.505 0.000 0.702 0.702
## eta_rec5 ~
## eta_rec4 (a2) 0.760 0.072 10.505 0.000 0.726 0.726
## eta_rec3 ~
## eta_math2 (c1) -0.017 0.122 -0.138 0.890 -0.013 -0.013
## eta_rec4 ~
## eta_math3 (c1) -0.017 0.122 -0.138 0.890 -0.012 -0.012
## eta_rec5 ~
## eta_math4 (c1) -0.017 0.122 -0.138 0.890 -0.011 -0.011
## eta_math3 ~
## eta_rec2 (c2) 0.073 0.116 0.626 0.531 0.059 0.059
## eta_math4 ~
## eta_rec3 (c2) 0.073 0.116 0.626 0.531 0.098 0.098
## eta_math5 ~
## eta_rec4 (c2) 0.073 0.116 0.626 0.531 0.113 0.113
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## i_math ~~
## i_rec (covi) 52.233 5.397 9.679 0.000 0.770 0.770
## et_mth2 0.000 0.000 0.000
## eta_rc2 0.000 0.000 0.000
## i_rec ~~
## et_mth2 0.000 0.000 0.000
## eta_rc2 0.000 0.000 0.000
## eta_math2 ~~
## eta_rc2 (cov1) 9.755 4.268 2.286 0.022 0.298 0.298
## .eta_math3 ~~
## .eta_rc3 (e1) 16.649 2.909 5.724 0.000 0.362 0.362
## .eta_math4 ~~
## .eta_rc4 (e1) 16.649 2.909 5.724 0.000 0.443 0.443
## .eta_math5 ~~
## .eta_rc5 (e1) 16.649 2.909 5.724 0.000 0.467 0.467
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .math2 0.000 0.000 0.000
## .math3 0.000 0.000 0.000
## .math4 0.000 0.000 0.000
## .math5 0.000 0.000 0.000
## .rec2 0.000 0.000 0.000
## .rec3 0.000 0.000 0.000
## .rec4 0.000 0.000 0.000
## .rec5 0.000 0.000 0.000
## eta_math2 -5.870 3.083 -1.904 0.057 -0.954 -0.954
## .eta_math3 -2.274 0.933 -2.437 0.015 -0.347 -0.347
## .eta_math4 4.589 1.825 2.514 0.012 0.755 0.755
## .eta_math5 9.295 2.848 3.264 0.001 1.634 1.634
## eta_rec2 39.767 1.076 36.957 0.000 7.486 7.486
## .eta_rec3 16.910 2.792 6.057 0.000 2.063 2.063
## .eta_rec4 17.350 3.135 5.534 0.000 1.957 1.957
## .eta_rec5 18.811 3.636 5.174 0.000 2.027 2.027
## i_math 38.171 3.150 12.117 0.000 4.852 4.852
## i_rec -5.773 1.117 -5.168 0.000 -0.669 -0.669
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .math2 0.000 0.000 0.000
## .math3 0.000 0.000 0.000
## .math4 0.000 0.000 0.000
## .math5 0.000 0.000 0.000
## .rec2 0.000 0.000 0.000
## .rec3 0.000 0.000 0.000
## .rec4 0.000 0.000 0.000
## .rec5 0.000 0.000 0.000
## et_mth2 (vrm2) 37.895 4.822 7.858 0.000 1.000 1.000
## .et_mth3 (vrm3) 41.433 4.684 8.845 0.000 0.964 0.964
## .et_mth4 (vrm4) 35.127 6.268 5.604 0.000 0.952 0.952
## .et_mth5 (vrm5) 30.826 4.332 7.115 0.000 0.953 0.953
## eta_rc2 (vrr2) 28.218 11.973 2.357 0.018 1.000 1.000
## .eta_rc3 (vrr3) 51.159 8.128 6.294 0.000 0.761 0.761
## .eta_rc4 (vrr4) 40.250 5.465 7.365 0.000 0.512 0.512
## .eta_rc5 (vrr5) 41.229 5.885 7.006 0.000 0.479 0.479
## i_math (vrmt) 61.885 4.073 15.195 0.000 1.000 1.000
## i_rec (vrrc) 74.439 10.943 6.803 0.000 1.000 1.000
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## cori 0.770 0.064 12.041 0.000 0.770 0.770
## cor1 0.298 0.119 2.503 0.012 0.298 0.298
## cort2 0.362 0.063 5.784 0.000 0.422 0.422
## cort3 0.443 0.080 5.564 0.000 0.518 0.518
## cort4 0.467 0.083 5.625 0.000 0.535 0.535
semPaths(riclpm_fit, "std", layout = "tree2", edge.label.cex = 0.8, curvePivot = TRUE)
