library(pacman); p_load(psych, lavaan, dplyr)
FITM <- c("chisq", "df", "nPar", "cfi", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper", "aic", "bic")
psych::describe(data)Karwowski & Milerski (2021) should be applauded for analyzing and providing excellent data for investigating the effect of schooling on measured intelligence. They exploited a Polish educational reform that made education more intensive for younger students and then showed that it seemed they experienced larger gains than the kids one year older who were not exposed.
They did assess measurement invariance, but they only went up to scalar invariance. Scalar invariance is good enough when you’re trying to see if means are broadly comparable, but without invariant residual variances (i.e., strict factorial invariance), you’re in the dark on whether mean changes can be interpreted in common because the model will tend to conflate broad gains of all sorts that are in the residuals with latent gains, and this is all the more problematic if you use a simple model.
Using Karwowski & Milerski’s parcel-based model for g, I will conduct an MGCFA of their Time 1 and Time 2 data to assess whether there was evidence of baseline and endpoint differences, as well as differential gain for the younger, more intensively-schooled cohort. Differential gains will be assessed by comparing the middle schooler group at T1 to the primary schooler group at T2, since those samples will have similar ages and comparable academic experience, excluding supplementary schooling.
########################################################################
########## INITIAL TESTS - CFA FOR EACH TIME POINT SEPARATELY ############
########################################################################
IQparcT1 <- 'g =~ parc_anf1 + parc_wnf1 + parc_marf1 + parc_rotf1'
IQparcT2 <- 'g =~ parc_anf2 + parc_wnf2 + parc_marf2 + parc_rotf2'
###########################
###### SIMPLE MODEL #######
###########################
fit1 <- sem(IQparcT1, data=data, std.lv = T, missing = "FIML", meanstructure = TRUE)
fit2 <- sem(IQparcT2, data=data, std.lv = T, missing = "FIML", meanstructure = TRUE)
round(cbind(BASE1 = fitMeasures(fit1, FITM),
BASE2 = fitMeasures(fit2, FITM), 3))## BASE1 BASE2
## chisq 20 22 3
## df 2 2 3
## npar 12 12 3
## cfi 1 1 3
## rmsea 0 0 3
## rmsea.ci.lower 0 0 3
## rmsea.ci.upper 0 0 3
## aic 27499 27899 3
## bic 27566 27966 3#Their models
#############################################
######## CONFIGURAL INVARIANCE MODEL ########
#############################################
fitconf1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = T, group = "Cohort")
fitconf2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = T, group = "Cohort")
############################################
########## METRIC INVARIANCE MODEL #########
############################################
fitmetr1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings"))
fitmetr2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings"))
#########################################
####### SCALAR INVARIANCE MODEL #########
#########################################
fitsclr1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts"))
fitsclr2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts"))
#My models
#########################################
####### STRICT INVARIANCE MODEL #########
#########################################
fitstrct1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals"))
fitstrct2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals"))
#########################################
####### LATENT VARIANCE INVARIANCE MODEL #########
#########################################
fitlatv1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals", "lv.variances"))
fitlatv2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals", "lv.variances"))
#########################################
####### LATENT MEAN INVARIANCE MODEL #########
#########################################
fitlatm1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals", "lv.variances", "means"))
fitlatm2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals", "lv.variances", "means"))round(cbind(CONF1 = fitMeasures(fitconf1, FITM),
METR1 = fitMeasures(fitmetr1, FITM),
SCLR1 = fitMeasures(fitsclr1, FITM),
STRI1 = fitMeasures(fitstrct1, FITM),
LATV1 = fitMeasures(fitlatv1, FITM),
LATM1 = fitMeasures(fitlatm1, FITM)),3)## CONF1 METR1 SCLR1 STRI1 LATV1 LATM1
## chisq 20.471 27.987 32.357 46.516 52.854 59.940
## df 4.000 7.000 10.000 14.000 15.000 16.000
## npar 24.000 21.000 18.000 14.000 13.000 12.000
## cfi 0.992 0.990 0.989 0.984 0.982 0.979
## rmsea 0.065 0.055 0.048 0.049 0.051 0.053
## rmsea.ci.lower 0.039 0.035 0.030 0.034 0.036 0.039
## rmsea.ci.upper 0.094 0.077 0.067 0.065 0.066 0.068
## aic 27483.647 27485.163 27483.533 27489.691 27494.030 27499.115
## bic 27617.596 27602.369 27583.995 27567.828 27566.586 27566.090round(cbind(CONF2 = fitMeasures(fitconf2, FITM),
METR2 = fitMeasures(fitmetr2, FITM),
SCLR2 = fitMeasures(fitsclr2, FITM),
STRI2 = fitMeasures(fitstrct2, FITM),
LATV2 = fitMeasures(fitlatv2, FITM),
LATM2 = fitMeasures(fitlatm2, FITM)),3)## CONF2 METR2 SCLR2 STRI2 LATV2 LATM2
## chisq 22.585 25.971 33.995 48.085 49.253 49.276
## df 4.000 7.000 10.000 14.000 15.000 16.000
## npar 24.000 21.000 18.000 14.000 13.000 12.000
## cfi 0.994 0.994 0.992 0.989 0.989 0.989
## rmsea 0.069 0.053 0.049 0.050 0.048 0.046
## rmsea.ci.lower 0.043 0.032 0.032 0.035 0.034 0.032
## rmsea.ci.upper 0.098 0.075 0.068 0.066 0.064 0.061
## aic 27896.042 27893.429 27895.453 27901.543 27900.710 27898.734
## bic 28029.991 28010.634 27995.914 27979.680 27973.266 27965.708"Time 1"## [1] "Time 1"1-pchisq(fitMeasures(fitmetr1, "chisq")-fitMeasures(fitconf1, "chisq"), 3) #Metric## chisq
## 0.0571-pchisq(fitMeasures(fitsclr1, "chisq")-fitMeasures(fitmetr1, "chisq"), 3) #Scalar## chisq
## 0.2241-pchisq(fitMeasures(fitstrct1, "chisq")-fitMeasures(fitsclr1, "chisq"), 4) #Strict## chisq
## 0.0071-pchisq(fitMeasures(fitlatv1, "chisq")-fitMeasures(fitstrct1, "chisq"), 1) #Latent Variance## chisq
## 0.0121-pchisq(fitMeasures(fitlatm1, "chisq")-fitMeasures(fitlatv1, "chisq"), 1) #Latent Mean## chisq
## 0.008"Time 2"## [1] "Time 2"1-pchisq(fitMeasures(fitmetr2, "chisq")-fitMeasures(fitconf2, "chisq"), 3) #Metric## chisq
## 0.3361-pchisq(fitMeasures(fitsclr2, "chisq")-fitMeasures(fitmetr2, "chisq"), 3) #Scalar## chisq
## 0.0461-pchisq(fitMeasures(fitstrct2, "chisq")-fitMeasures(fitsclr2, "chisq"), 4) #Strict## chisq
## 0.0071-pchisq(fitMeasures(fitlatv2, "chisq")-fitMeasures(fitstrct2, "chisq"), 1) #Latent Variance## chisq
## 0.281-pchisq(fitMeasures(fitlatm2, "chisq")-fitMeasures(fitlatv2, "chisq"), 1) #Latent Mean## chisq
## 0.878#My models
#########################################
####### STRICT INVARIANCE MODEL #########
#########################################
pfitstrct1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals"),
group.partial = c("parc_rotf1 ~~ parc_rotf1"))
fitlatv1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals", "lv.variances"),
group.partial = c("parc_rotf1 ~~ parc_rotf1"))
fitlatm1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals", "means"),
group.partial = c("parc_rotf1 ~~ parc_rotf1"))
round(cbind(CONF1 = fitMeasures(fitconf1, FITM),
METR1 = fitMeasures(fitmetr1, FITM),
SCLR1 = fitMeasures(fitsclr1, FITM),
STRI1 = fitMeasures(fitstrct1, FITM),
PSTRI1 = fitMeasures(pfitstrct1, FITM),
LATV1 = fitMeasures(fitlatv1, FITM), #Drop Restraint
LATM1 = fitMeasures(fitlatm1, FITM)), 3) #Drop Restraint## CONF1 METR1 SCLR1 STRI1 PSTRI1 LATV1
## chisq 20.471 27.987 32.357 46.516 37.730 43.845
## df 4.000 7.000 10.000 14.000 13.000 14.000
## npar 24.000 21.000 18.000 14.000 15.000 14.000
## cfi 0.992 0.990 0.989 0.984 0.988 0.986
## rmsea 0.065 0.055 0.048 0.049 0.044 0.047
## rmsea.ci.lower 0.039 0.035 0.030 0.034 0.028 0.031
## rmsea.ci.upper 0.094 0.077 0.067 0.065 0.061 0.063
## aic 27483.647 27485.163 27483.533 27489.691 27482.905 27487.021
## bic 27617.596 27602.369 27583.995 27567.828 27566.624 27565.158
## LATM1
## chisq 44.819
## df 14.000
## npar 14.000
## cfi 0.985
## rmsea 0.047
## rmsea.ci.lower 0.032
## rmsea.ci.upper 0.063
## aic 27487.994
## bic 27566.1311-pchisq(fitMeasures(pfitstrct1, "chisq") - fitMeasures(fitsclr1, "chisq"), 3) #Partial Strict## chisq
## 0.1461-pchisq(fitMeasures(fitlatv1, "chisq") - fitMeasures(pfitstrct1, "chisq"), 1) #Latent Variances## chisq
## 0.0131-pchisq(fitMeasures(fitlatm1, "chisq") - fitMeasures(pfitstrct1, "chisq"), 1) #Means## chisq
## 0.008summary(pfitstrct1, fit = T, stand = T)## lavaan 0.6.14 ended normally after 31 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 25
## Number of equality constraints 10
##
## Number of observations per group:
## 1,00 985
## 2,00 976
## Number of missing patterns per group:
## 1,00 1
## 2,00 1
##
## Model Test User Model:
##
## Test statistic 37.730
## Degrees of freedom 13
## P-value (Chi-square) 0.000
## Test statistic for each group:
## 1,00 19.238
## 2,00 18.492
##
## Model Test Baseline Model:
##
## Test statistic 2075.599
## Degrees of freedom 12
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.988
## Tucker-Lewis Index (TLI) 0.989
##
## Robust Comparative Fit Index (CFI) 0.988
## Robust Tucker-Lewis Index (TLI) 0.989
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -13726.453
## Loglikelihood unrestricted model (H1) NA
##
## Akaike (AIC) 27482.905
## Bayesian (BIC) 27566.624
## Sample-size adjusted Bayesian (SABIC) 27518.968
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.044
## 90 Percent confidence interval - lower 0.028
## 90 Percent confidence interval - upper 0.061
## P-value H_0: RMSEA <= 0.050 0.700
## P-value H_0: RMSEA >= 0.080 0.000
##
## Robust RMSEA 0.044
## 90 Percent confidence interval - lower 0.028
## 90 Percent confidence interval - upper 0.061
## P-value H_0: Robust RMSEA <= 0.050 0.700
## P-value H_0: Robust RMSEA >= 0.080 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.029
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
##
## Group 1 [1,00]:
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## g =~
## prc_nf1 1.000 0.781 0.692
## prc_wn1 (.p2.) 2.111 0.077 27.436 0.000 1.648 0.787
## prc_mr1 (.p3.) 1.956 0.076 25.648 0.000 1.527 0.695
## prc_rt1 (.p4.) 0.626 0.035 17.782 0.000 0.489 0.457
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .prc_nf1 (.10.) 3.157 0.032 97.766 0.000 3.157 2.797
## .prc_wn1 (.11.) 4.470 0.063 70.860 0.000 4.470 2.135
## .prc_mr1 (.12.) 4.153 0.063 65.919 0.000 4.153 1.891
## .prc_rt1 (.13.) 2.592 0.028 92.832 0.000 2.592 2.422
## g 0.000 0.000 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .prc_nf1 (.p5.) 0.665 0.029 22.849 0.000 0.665 0.522
## .prc_wn1 (.p6.) 1.671 0.100 16.778 0.000 1.671 0.381
## .prc_mr1 (.p7.) 2.490 0.110 22.663 0.000 2.490 0.516
## .prc_rt1 0.906 0.044 20.548 0.000 0.906 0.791
## g 0.610 0.045 13.409 0.000 1.000 1.000
##
##
## Group 2 [2,00]:
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## g =~
## prc_nf1 1.000 0.861 0.726
## prc_wn1 (.p2.) 2.111 0.077 27.436 0.000 1.818 0.815
## prc_mr1 (.p3.) 1.956 0.076 25.648 0.000 1.685 0.730
## prc_rt1 (.p4.) 0.626 0.035 17.782 0.000 0.539 0.456
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .prc_nf1 (.10.) 3.157 0.032 97.766 0.000 3.157 2.662
## .prc_wn1 (.11.) 4.470 0.063 70.860 0.000 4.470 2.004
## .prc_mr1 (.12.) 4.153 0.063 65.919 0.000 4.153 1.799
## .prc_rt1 (.13.) 2.592 0.028 92.832 0.000 2.592 2.192
## g 0.110 0.042 2.659 0.008 0.128 0.128
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .prc_nf1 (.p5.) 0.665 0.029 22.849 0.000 0.665 0.472
## .prc_wn1 (.p6.) 1.671 0.100 16.778 0.000 1.671 0.336
## .prc_mr1 (.p7.) 2.490 0.110 22.663 0.000 2.490 0.467
## .prc_rt1 1.107 0.053 20.960 0.000 1.107 0.792
## g 0.742 0.054 13.641 0.000 1.000 1.000pfitstrct2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals"),
group.partial = c("parc_rotf2 ~~ parc_rotf2", "parc_wnf2 ~~ parc_wnf2"))
fitlatv2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals", "lv.variances"),
group.partial = c("parc_rotf2 ~~ parc_rotf2", "parc_wnf2 ~~ parc_wnf2"))
fitlatm2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals", "lv.variances", "means"),
group.partial = c("parc_rotf2 ~~ parc_rotf2", "parc_wnf2 ~~ parc_wnf2"))
round(cbind(CONF2 = fitMeasures(fitconf2, FITM),
METR2 = fitMeasures(fitmetr2, FITM),
SCLR2 = fitMeasures(fitsclr2, FITM),
STRI2 = fitMeasures(fitstrct2, FITM),
PSTRI2 = fitMeasures(pfitstrct2, FITM),
LATV2 = fitMeasures(fitlatv2, FITM),
LATM2 = fitMeasures(fitlatm2, FITM)),3)## CONF2 METR2 SCLR2 STRI2 PSTRI2 LATV2
## chisq 22.585 25.971 33.995 48.085 38.633 39.482
## df 4.000 7.000 10.000 14.000 12.000 13.000
## npar 24.000 21.000 18.000 14.000 16.000 15.000
## cfi 0.994 0.994 0.992 0.989 0.991 0.991
## rmsea 0.069 0.053 0.049 0.050 0.048 0.046
## rmsea.ci.lower 0.043 0.032 0.032 0.035 0.031 0.030
## rmsea.ci.upper 0.098 0.075 0.068 0.066 0.065 0.062
## aic 27896.042 27893.429 27895.453 27901.543 27896.091 27894.939
## bic 28029.991 28010.634 27995.914 27979.680 27985.390 27978.657
## LATM2
## chisq 39.506
## df 14.000
## npar 14.000
## cfi 0.992
## rmsea 0.043
## rmsea.ci.lower 0.028
## rmsea.ci.upper 0.059
## aic 27892.963
## bic 27971.1001-pchisq(fitMeasures(pfitstrct2, "chisq") - fitMeasures(fitsclr2, "chisq"), 2) #Partial Strict## chisq
## 0.0981-pchisq(fitMeasures(fitlatv2, "chisq") - fitMeasures(pfitstrct2, "chisq"), 1) #Latent Variances## chisq
## 0.3571-pchisq(fitMeasures(fitlatm2, "chisq") - fitMeasures(fitlatv2, "chisq"), 1) #Means## chisq
## 0.877df <- data[data$Cohort %in% c("1,00", "2,00"), ]
df$parc_anf <- ifelse(df$Cohort == "1,00", df$parc_anf2, df$parc_anf1)
df$parc_wnf <- ifelse(df$Cohort == "1,00", df$parc_wnf2, df$parc_wnf1)
df$parc_marf <- ifelse(df$Cohort == "1,00", df$parc_marf2, df$parc_marf1)
df$parc_rotf <- ifelse(df$Cohort == "1,00", df$parc_rotf2, df$parc_rotf1)
df <- df[, c("Cohort", "parc_anf", "parc_wnf", "parc_marf", "parc_rotf")]
IQparc <- 'g =~ parc_anf + parc_wnf + parc_marf + parc_rotf'
#############################################
######## CONFIGURAL INVARIANCE MODEL ########
#############################################
fitconf <- cfa(IQparc, data=df, missing = "FIML", meanstructure = TRUE, std.lv = T, group = "Cohort")
############################################
########## METRIC INVARIANCE MODEL #########
############################################
fitmetr <- cfa(IQparc, data=df, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings"))
pfitmetr <- cfa(IQparc, data=df, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings"),
group.partial = c("g =~ parc_rotf"))
#########################################
####### SCALAR INVARIANCE MODEL #########
#########################################
fitsclr <- cfa(IQparc, data=df, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts"),
group.partial = c("g =~ parc_rotf"))
#########################################
####### STRICT INVARIANCE MODEL #########
#########################################
fitstrct <- cfa(IQparc, data=df, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals"),
group.partial = c("g =~ parc_rotf"))
#########################################
####### LATENT VARIANCE INVARIANCE MODEL #########
#########################################
fitlatv <- cfa(IQparc, data=df, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals", "lv.variances"),
group.partial = c("g =~ parc_rotf"))
#########################################
####### LATENT MEAN INVARIANCE MODEL #########
#########################################
fitlatm <- cfa(IQparc, data=df, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals", "means"),
group.partial = c("g =~ parc_rotf"))
round(cbind(CONF = fitMeasures(fitconf, FITM),
METR = fitMeasures(fitmetr, FITM),
PMETR = fitMeasures(pfitmetr, FITM),
SCLR = fitMeasures(fitsclr, FITM),
STRI = fitMeasures(fitstrct, FITM),
LATV = fitMeasures(fitlatv, FITM), #Drop constraint
LATM = fitMeasures(fitlatm, FITM)), 3)## CONF METR PMETR SCLR STRI LATV
## chisq 21.428 29.944 22.401 24.828 30.109 36.915
## df 4.000 7.000 6.000 9.000 13.000 14.000
## npar 24.000 21.000 22.000 19.000 15.000 14.000
## cfi 0.993 0.991 0.994 0.994 0.993 0.991
## rmsea 0.067 0.058 0.053 0.042 0.037 0.041
## rmsea.ci.lower 0.041 0.037 0.031 0.023 0.019 0.025
## rmsea.ci.upper 0.096 0.080 0.077 0.063 0.054 0.057
## aic 27814.662 27817.178 27811.635 27808.062 27805.343 27810.149
## bic 27948.611 27934.383 27934.422 27914.105 27889.061 27888.286
## LATM
## chisq 31.403
## df 14.000
## npar 14.000
## cfi 0.993
## rmsea 0.036
## rmsea.ci.lower 0.019
## rmsea.ci.upper 0.052
## aic 27804.637
## bic 27882.7741-pchisq(fitMeasures(fitmetr, "chisq") - fitMeasures(fitconf, "chisq"), 3) #Metric## chisq
## 0.0361-pchisq(fitMeasures(pfitmetr, "chisq") - fitMeasures(fitconf, "chisq"), 2) #Metric## chisq
## 0.6151-pchisq(fitMeasures(fitsclr, "chisq") - fitMeasures(pfitmetr, "chisq"), 3) #Scalar## chisq
## 0.4891-pchisq(fitMeasures(fitstrct, "chisq") - fitMeasures(fitsclr, "chisq"), 4) #Strict## chisq
## 0.261-pchisq(fitMeasures(fitlatv, "chisq") - fitMeasures(fitstrct, "chisq"), 1) #Latent Variance## chisq
## 0.0091-pchisq(fitMeasures(fitlatm, "chisq") - fitMeasures(fitstrct, "chisq"), 1) #Latent Mean## chisq
## 0.255Oddly, mental rotation was influenced by different variables or had different reliability at both timepoints. I do not know why this was. At time two, there were different influences on reasoning too. I do not know why this happened, but partially invariant models had to be fit. There was meager evidence for a difference in intercepts at T2 (p = 0.046), and consistent with that, this had no effect on results, so I did not pursue a partially invariant scalar/strong model for this page. At T1, the latent variance of g differed, and again, I do not know why, but unfortunately, that needed to be released. At T2, the latent variances were equal.
The lack of longitudinal invariance for the exposed primary school cohort suggested that there was likely some effect, but it was expressed through changing the interpretation of indicator variables. Despite this, the levels of the items were still consistently interpreted.
The most surprising result came from the explicit test of whether the latent mean of g differed at T1 or T2. The answers differed: at T1, Group 2 was more intelligent, but at T2, they were equal. Because the differences dissolved, whatever impact intensive schooling might have had - and Karwowski & Milerski maintained it was “not systematic and quite unstable” - may have been general.
But the final analysis assessed whether this was really just the rate of developmental change being greater at younger ages by comparing Primary T2 and Middle T1. The problem with this is the retesting effect, but the assessment’s retest bonus should have at least faded a bit, if not totally, over the school year. With this said, the result was clear: changes in g as a result of more intensive education were not supported The primary school children’s mean for g was not distinguishable from the preexisting differences observed in the cohort that was not exposed to the more intensive primary schooling regime.
There was negative evidence that more intensive education increased the rate of development of g. As noted in the postscript below, verbal results were inconclusive because of poor model fit. Nonverbal results were consistent with the results for g, as the primary schooler T2 mean did not differ from the middle schooler T1 mean (p = 0.171).
Karwowski, M., & Milerski, B. (2021). Intensive schooling and cognitive ability: A case of Polish educational reform. Personality and Individual Differences, 183, 111121. https://doi.org/10.1016/j.paid.2021.111121
Data: https://osf.io/hm3df/
Note that because these are items, the estimator is wrong. However, DWLS did not seem to want to converge and WLSMV couldn’t handle the missingness well. Nonverbal first, since that was what they listed first, and what seemed to show catch-up gains in their study (see Fig. 2.).
########################################################################
########## INITIAL TESTS - CFA FOR EACH TIME POINT SEPARATELY ############
########################################################################
IQparcT1 <- 'g =~ mat1ok_f1 + mat2ok_f1 + mat3ok_f1 + mat4ok_f1 + mat5ok_f1 + mat6ok_f1 +
mat8ok_f1 + mat9ok_f1 + mat10ok_f1 + mat11ok_f1'
IQparcT2 <- 'g =~ mat1ok_f2 + mat2ok_f2 + mat3ok_f2 + mat4ok_f2 + mat5ok_f2 + mat6ok_f2 +
mat8ok_f2 + mat9ok_f2 + mat10ok_f2 + mat11ok_f2'
###########################
###### SIMPLE MODEL #######
###########################
fit1 <- sem(IQparcT1, data=data, std.lv = T, missing = "FIML", meanstructure = TRUE)
fit2 <- sem(IQparcT2, data=data, std.lv = T, missing = "FIML", meanstructure = TRUE)
round(cbind(BASE1 = fitMeasures(fit1, FITM),
BASE2 = fitMeasures(fit2, FITM), 3))## BASE1 BASE2
## chisq 122 115 3
## df 35 35 3
## npar 30 30 3
## cfi 1 1 3
## rmsea 0 0 3
## rmsea.ci.lower 0 0 3
## rmsea.ci.upper 0 0 3
## aic 19068 18789 3
## bic 19232 18952 3#Their models
#############################################
######## CONFIGURAL INVARIANCE MODEL ########
#############################################
fitconf1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = T, group = "Cohort")
fitconf2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = T, group = "Cohort")
############################################
########## METRIC INVARIANCE MODEL #########
############################################
fitmetr1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings"),
group.partial = c("g =~ mat4ok_f1"))
fitmetr2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings"))
#########################################
####### SCALAR INVARIANCE MODEL #########
#########################################
fitsclr1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts"),
group.partial = c("g =~ mat4ok_f1"))
fitsclr2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts"))
#My models
#########################################
####### STRICT INVARIANCE MODEL #########
#########################################
fitstrct1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals"),
group.partial = c("g =~ mat4ok_f1"))
fitstrct2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals"))
#########################################
####### LATENT VARIANCE INVARIANCE MODEL #########
#########################################
fitlatv1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals", "lv.variances"),
group.partial = c("g =~ mat4ok_f1"))
fitlatv2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals", "lv.variances"))
#########################################
####### LATENT MEAN INVARIANCE MODEL #########
#########################################
fitlatm1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals", "means"),
group.partial = c("g =~ mat4ok_f1"))
fitlatm2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals", "lv.variances", "means"))
round(cbind(CONF1 = fitMeasures(fitconf1, FITM),
METR1 = fitMeasures(fitmetr1, FITM),
SCLR1 = fitMeasures(fitsclr1, FITM),
STRI1 = fitMeasures(fitstrct1, FITM),
LATV1 = fitMeasures(fitlatv1, FITM), #Drop constraint
LATM1 = fitMeasures(fitlatm1, FITM)),3)## CONF1 METR1 SCLR1 STRI1 LATV1 LATM1
## chisq 165.004 176.453 183.180 190.882 193.339 206.423
## df 70.000 78.000 87.000 97.000 98.000 98.000
## npar 60.000 52.000 43.000 33.000 32.000 32.000
## cfi 0.915 0.912 0.914 0.916 0.915 0.904
## rmsea 0.039 0.038 0.035 0.033 0.033 0.035
## rmsea.ci.lower 0.032 0.030 0.028 0.026 0.026 0.029
## rmsea.ci.upper 0.047 0.045 0.043 0.040 0.040 0.042
## aic 19073.523 19068.973 19057.700 19045.402 19045.859 19058.943
## bic 19401.805 19353.484 19292.968 19225.957 19220.942 19234.026round(cbind(CONF2 = fitMeasures(fitconf2, FITM),
METR2 = fitMeasures(fitmetr2, FITM),
SCLR2 = fitMeasures(fitsclr2, FITM),
STRI2 = fitMeasures(fitstrct2, FITM),
LATV2 = fitMeasures(fitlatv2, FITM),
LATM2 = fitMeasures(fitlatm2, FITM)),3)## CONF2 METR2 SCLR2 STRI2 LATV2 LATM2
## chisq 172.186 176.635 186.943 193.128 194.316 194.316
## df 70.000 79.000 88.000 98.000 99.000 100.000
## npar 60.000 51.000 42.000 32.000 31.000 30.000
## cfi 0.930 0.933 0.932 0.935 0.934 0.935
## rmsea 0.042 0.038 0.037 0.034 0.034 0.034
## rmsea.ci.lower 0.034 0.031 0.029 0.027 0.027 0.026
## rmsea.ci.upper 0.050 0.046 0.044 0.041 0.041 0.041
## aic 18827.027 18813.476 18805.785 18791.969 18791.157 18789.158
## bic 19152.334 19089.987 19033.499 18965.466 18959.232 18951.811"Time 1"## [1] "Time 1"1-pchisq(fitMeasures(fitmetr1, "chisq")-fitMeasures(fitconf1, "chisq"), 8) #Metric## chisq
## 0.1781-pchisq(fitMeasures(fitsclr1, "chisq")-fitMeasures(fitmetr1, "chisq"), 9) #Scalar## chisq
## 0.6661-pchisq(fitMeasures(fitstrct1, "chisq")-fitMeasures(fitsclr1, "chisq"), 9) #Strict## chisq
## 0.5641-pchisq(fitMeasures(fitlatv1, "chisq")-fitMeasures(fitstrct1, "chisq"), 1) #Latent Variance## chisq
## 0.1171-pchisq(fitMeasures(fitlatm1, "chisq")-fitMeasures(fitstrct1, "chisq"), 1) #Latent Mean## chisq
## 0"Time 2"## [1] "Time 2"1-pchisq(fitMeasures(fitmetr2, "chisq")-fitMeasures(fitconf2, "chisq"), 9) #Metric## chisq
## 0.8791-pchisq(fitMeasures(fitsclr2, "chisq")-fitMeasures(fitmetr2, "chisq"), 9) #Scalar## chisq
## 0.3261-pchisq(fitMeasures(fitstrct2, "chisq")-fitMeasures(fitsclr2, "chisq"), 9) #Strict## chisq
## 0.7211-pchisq(fitMeasures(fitlatv2, "chisq")-fitMeasures(fitstrct2, "chisq"), 1) #Latent Variance## chisq
## 0.2761-pchisq(fitMeasures(fitlatm2, "chisq")-fitMeasures(fitlatv2, "chisq"), 1) #Latent Mean## chisq
## 0.981For the nonverbal results, the difference observed at T1 was real. The difference in the interpretation of the fourth matrices item was substantial, with Group 1 having a loading of 0.584 versus 0.892 for Group 2. This was not enough to explain the observed performance gaps and the latent mean difference ended up being 0.221 d in the strict model. The latent variance model could not be used, since the latent variances were unequal.
The T2 results were fully invariant, so there were real catch-up gains in nonverbal ability. I do not know why, because unfortunately, the way this was measured, “nonverbal intelligence” is full of variance from g. A more complex model is required, but there are not enough indicators in the dataset to do that. A more important issue also affected the verbal model, so I will discuss that later.
Before moving on, I should note that the authors mentioned in their code that all items modeled fit better if treated as categorical, and I think they should be scolded for this. The reason is because they changed the estimator and the interpretation of the items by using the wrong estimator. Do not treat variables that are not actually categorical as if they are categorical. You will produce a misfitted model because you used an inappropriate estimator. This is a much worse mistake than, say, using ML instead of DWLS or WLSMV when there are convergence issues.
On to verbal:
########################################################################
########## INITIAL TESTS - CFA FOR EACH TIME POINT SEPARATELY ############
########################################################################
IQparcT1 <- 'g =~ an1ok_f1 + an2ok_f1 + an4ok_f1 + wn4ok_f1 + wn8ok_f1 +
wn10ok_f1 + wn13ok_f1 + wn18ok_f1 + wn20ok_f1 + wn23ok_f1'
IQparcT2 <- 'g =~ an1ok_f2 + an2ok_f2 + an4ok_f2 + wn4ok_f2 + wn8ok_f2 +
wn10ok_f2 + wn13ok_f2 + wn18ok_f2 + wn20ok_f2 + wn23ok_f2'
###########################
###### SIMPLE MODEL #######
###########################
fit1 <- sem(IQparcT1, data=data, std.lv = T, missing = "FIML", meanstructure = TRUE)
fit2 <- sem(IQparcT2, data=data, std.lv = T, missing = "FIML", meanstructure = TRUE)
round(cbind(BASE1 = fitMeasures(fit1, FITM),
BASE2 = fitMeasures(fit2, FITM), 3))## BASE1 BASE2
## chisq 175 170 3
## df 35 35 3
## npar 30 30 3
## cfi 1 1 3
## rmsea 0 0 3
## rmsea.ci.lower 0 0 3
## rmsea.ci.upper 0 0 3
## aic 22603 21304 3
## bic 22768 21467 3#Their models
#############################################
######## CONFIGURAL INVARIANCE MODEL ########
#############################################
fitconf1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = T, group = "Cohort")
fitconf2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = T, group = "Cohort")
############################################
########## METRIC INVARIANCE MODEL #########
############################################
fitmetr1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings"))
fitmetr2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings"),
group.partial = c("g =~ wn23ok_f2"))
#########################################
####### SCALAR INVARIANCE MODEL #########
#########################################
fitsclr1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts"))
fitsclr2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts"),
group.partial = c("g =~ wn23ok_f2"))
#My models
#########################################
####### STRICT INVARIANCE MODEL #########
#########################################
fitstrct1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals"))
fitstrct2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals"),
group.partial = c("g =~ wn23ok_f2"))
#########################################
####### LATENT VARIANCE INVARIANCE MODEL #########
#########################################
fitlatv1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals", "lv.variances"))
fitlatv2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals", "lv.variances"),
group.partial = c("g =~ wn23ok_f2"))
#########################################
####### LATENT MEAN INVARIANCE MODEL #########
#########################################
fitlatm1 <- cfa(IQparcT1, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals", "lv.variances", "means"))
fitlatm2 <- cfa(IQparcT2, data=data, missing = "FIML", meanstructure = TRUE, std.lv = F, group = "Cohort",
group.equal = c("loadings", "intercepts", "residuals", "lv.variances", "means"),
group.partial = c("g =~ wn23ok_f2"))
round(cbind(CONF1 = fitMeasures(fitconf1, FITM),
METR1 = fitMeasures(fitmetr1, FITM),
SCLR1 = fitMeasures(fitsclr1, FITM),
STRI1 = fitMeasures(fitstrct1, FITM),
LATV1 = fitMeasures(fitlatv1, FITM),
LATM1 = fitMeasures(fitlatm1, FITM)),3)## CONF1 METR1 SCLR1 STRI1 LATV1 LATM1
## chisq 243.292 269.419 293.461 299.487 303.784 306.623
## df 70.000 79.000 88.000 98.000 99.000 100.000
## npar 60.000 51.000 42.000 32.000 31.000 30.000
## cfi 0.892 0.882 0.872 0.875 0.873 0.872
## rmsea 0.053 0.052 0.052 0.048 0.049 0.048
## rmsea.ci.lower 0.046 0.046 0.045 0.042 0.042 0.042
## rmsea.ci.upper 0.060 0.059 0.058 0.055 0.055 0.055
## aic 22600.088 22608.216 22614.258 22600.284 22602.581 22603.420
## bic 22928.370 22887.255 22844.055 22775.367 22772.193 22767.561round(cbind(CONF2 = fitMeasures(fitconf2, FITM),
METR2 = fitMeasures(fitmetr2, FITM),
SCLR2 = fitMeasures(fitsclr2, FITM),
STRI2 = fitMeasures(fitstrct2, FITM),
LATV2 = fitMeasures(fitlatv2, FITM),
LATM2 = fitMeasures(fitlatm2, FITM)),3)## CONF2 METR2 SCLR2 STRI2 LATV2 LATM2
## chisq 218.410 226.658 244.478 255.113 255.267 257.847
## df 70.000 78.000 87.000 97.000 98.000 99.000
## npar 60.000 52.000 43.000 33.000 32.000 31.000
## cfi 0.939 0.938 0.935 0.935 0.935 0.934
## rmsea 0.050 0.048 0.047 0.044 0.044 0.044
## rmsea.ci.lower 0.043 0.041 0.040 0.038 0.037 0.037
## rmsea.ci.upper 0.058 0.055 0.053 0.051 0.050 0.050
## aic 21310.845 21303.092 21302.912 21293.548 21291.702 21292.282
## bic 21636.151 21585.025 21536.049 21472.466 21465.199 21460.357"Time 1"## [1] "Time 1"1-pchisq(fitMeasures(fitmetr1, "chisq")-fitMeasures(fitconf1, "chisq"), 9) #Metric## chisq
## 0.0021-pchisq(fitMeasures(fitsclr1, "chisq")-fitMeasures(fitmetr1, "chisq"), 9) #Scalar## chisq
## 0.0041-pchisq(fitMeasures(fitstrct1, "chisq")-fitMeasures(fitsclr1, "chisq"), 9) #Strict## chisq
## 0.7371-pchisq(fitMeasures(fitlatv1, "chisq")-fitMeasures(fitstrct1, "chisq"), 1) #Latent Variance## chisq
## 0.0381-pchisq(fitMeasures(fitlatm1, "chisq")-fitMeasures(fitlatv1, "chisq"), 1) #Latent Mean## chisq
## 0.092"Time 2"## [1] "Time 2"1-pchisq(fitMeasures(fitmetr2, "chisq")-fitMeasures(fitconf2, "chisq"), 8) #Metric## chisq
## 0.411-pchisq(fitMeasures(fitsclr2, "chisq")-fitMeasures(fitmetr2, "chisq"), 9) #Scalar## chisq
## 0.0371-pchisq(fitMeasures(fitstrct2, "chisq")-fitMeasures(fitsclr2, "chisq"), 9) #Strict## chisq
## 0.3021-pchisq(fitMeasures(fitlatv2, "chisq")-fitMeasures(fitstrct2, "chisq"), 1) #Latent Variance## chisq
## 0.6941-pchisq(fitMeasures(fitlatm2, "chisq")-fitMeasures(fitlatv2, "chisq"), 1) #Latent Mean## chisq
## 0.108The T1 verbal results could not be analyzed credibly because the initial model fit was very poor and too many loadings had to be released to fit a working metric model. The similarity in means at T1 was uninterpretable. The T2 verbal results could be analyzed credibly because the model fit well enough and only a single loading had to be released to make all subsequent levels of measurement invariance fit just fine. Everything just worked.
The net result of these two findings is that the effect of intensive schooling on group differences in verbal intelligence - which, like nonverbal intelligence, is polluted by g variance - is unclear.
sessionInfo()## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dplyr_1.1.0 lavaan_0.6-14 psych_2.2.9 pacman_0.5.1
##
## loaded via a namespace (and not attached):
## [1] bslib_0.4.2 compiler_4.2.2 pillar_1.8.1 jquerylib_0.1.4
## [5] tools_4.2.2 digest_0.6.31 tibble_3.1.8 jsonlite_1.8.4
## [9] evaluate_0.20 lifecycle_1.0.3 nlme_3.1-160 lattice_0.20-45
## [13] pkgconfig_2.0.3 rlang_1.0.6 cli_3.6.0 rstudioapi_0.14
## [17] yaml_2.3.7 parallel_4.2.2 pbivnorm_0.6.0 xfun_0.37
## [21] fastmap_1.1.0 knitr_1.42 generics_0.1.3 sass_0.4.5
## [25] vctrs_0.5.2 tidyselect_1.2.0 stats4_4.2.2 grid_4.2.2
## [29] glue_1.6.2 R6_2.5.1 fansi_1.0.4 rmarkdown_2.20
## [33] magrittr_2.0.3 MASS_7.3-58.1 htmltools_0.5.4 mnormt_2.1.1
## [37] quadprog_1.5-8 utf8_1.2.3 cachem_1.0.6