Very much a work-in-progress.

Setup

Packages

library(pacman)
p_load(kirkegaard, dplyr, psych, lavaan, semPlot, umx)

Functions

#From Gunn, Grimm & Edwards (2019)

SDI2.UDI2 = function(p, loading1, loading2, intercept1, intercept2, stdev2, fmean2, fsd2, nodewidth = 0.01, lowerLV = -5, upperLV = 5) {
   LV = seq(lowerLV,upperLV,nodewidth)
   DiffExpItem12 = matrix(NA,length(LV),p)
   pdfLV2 = matrix(NA,length(LV),1)
   SDI2numerator = matrix(NA,length(LV),p)
   UDI2numerator = matrix(NA,length(LV),p)
   SDI2 = matrix(NA,p,1)
   UDI2 = matrix(NA,p,1)
   for(j in 1:p){
    for(k in 1:length(LV)){
      DiffExpItem12[k,j] <- (intercept1[j]-intercept2[j]) + (loading1[j]-loading2[j])*LV[k]
      pdfLV2[k] = dnorm(LV[k], mean=fmean2, sd=fsd2)
      SDI2numerator[k,j] = DiffExpItem12[k,j]*pdfLV2[k]*nodewidth
      UDI2numerator[k,j] = abs(SDI2numerator[k,j])
    }
    SDI2[j] <- sum(SDI2numerator[,j])/stdev2[j]
    UDI2[j] <- sum(UDI2numerator[,j])/stdev2[j]
   }
  colnames(SDI2) = "SDI2"
  colnames(UDI2) = "UDI2"
  return(list(SDI2=round(SDI2,4), UDI2=round(UDI2,4)))}

#Based on Naaman (2016)

NZ <- function(N, S) {
  NZ = qnorm(1-(N^-(6/5))/S) #'1-' included instead of 'abs()'; change it if you don't like it
  return(NZ)}

#Based on work from Millsap & Olivera-Aguilar (2012); also works for intercepts

MOES <- function(T1, T2, S1, S2) {
  ES = (T1 - T2) / (S1^2 - S2^2)
  return(ES)}

Data

Too large for DT :( - you can get the data at the osf.

#VES <- read.csv(".../VES_dataset.csv")

VES %<>% rename(
  WBD = WAIS_BD, #WAIS Block Design
  WGI = WAIS_GI, #WAIS General Information
  GPTR = GPT_right, #Grooved Pegboard Task - Right Hand
  GPTL = GPT_left, #Grooved Pegboard Task - Left Hand
  CD = copy_direct, #Rey-Osterrieth Complex Figure Drawing - Direct, i.e., in sight
  CI = copy_immediate, #Complex Figure Drawing - Immediate(ly after removal from sight)
  CY = copy_delayed, #Complex Figure Drawing - Delayed, i.e., asked to draw the figure after time has passed from viewing
  ACVE = ACB_verbal_early, #Army Classification Battery Verbal subtest administered at introduction
  ACVL = ACB_verbal_later, #Verbal subtest administered at follow-up
  ACAE = ACB_arithmetic_early, #Arithmetic subtest administered at introduction
  ACAL = ACB_arithmetic_later) #Arithmetic subtest administered at follow-up

VES$CC = (standardize(VES$CI) + standardize(VES$CY))/2 #Very similar, which unnecessarily harms path model fit.

early_tests = c("ACVE", "ACAE", "PA", "GIT", "AFQT") #Tests taken at introduction
later_tests = c("WRAT", "CVLT", "WCST", "WBD", "WGI", "GPTR", "GPTL", "PASAT", "CD", "CC", "WLGT", "ACVL", "ACAL") #Tests taken at follow-up. Unlisted above are the Wisconsin Card Sorting Test (WCST), California Verbal Learning Test (CVLT), Wide-Range Achievement Test (WRAT), Word List Generation Test (WLGT), and Paced Auditory Serial Addition Test (PASAT)
tests = c(early_tests, later_tests) #All tests

#Model fit measures

FITM <- c("chisq", "df", "nPar", "cfi", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper", "aic", "bic")

Rationale

Replication of Ritchie, Bates & Deary (2015) including measurement invariance assessment (MI) via (1) MGCFA; (2) LSEM; (3) MNLFA. Current upload (24/5/2020) is just the headline result and replication, not the additional analyses, which include testing MI, gains to the sumscore, comparisons of apparent and real effects on \(g\), assessment of the (lack of) relationship of education with sensory measures, computation of whether residual variances change for more educated groups as expected under an effect model, and computation of “economic gains” (increased income) due to education and their (lack of) mediation by cognitive improvement. More analyses may be performed if I think of them or people care. This dataset will also be used to compare factor and network models and at that point network effects of education will be assessed (see: Fried; note: add cites).

Analysis

#describe(VES) <- too long to look good with chosen theme.

Path Model

Five tests were administered at entrance and an additional 15 at follow-up. Two of these - the immediated and delayed Rey-Osterrieth Complex Figure Drawing tasks - are highly intercorrelated and specific, so they’ve been combined (see ‘Data’ above), leaving us with 18 items where \(\mu = 100, \sigma = 15\). For this initial analysis and upload, I’ll output the three models Ritchie, Bates & Deary (2015) assessed. Briefly, these are (A) education effects only on \(g\); (B) education effects on both \(g\) and individual subtests, and; (C) education effects solely relegated to individual subtests.

The variable in the path models below, ‘\(gl\)’, corresponds to the general factor of intelligence measured at follow-up, and ‘\(ge\)’ is the same variable measured at induction into the study. \(ge\) is regressed on \(gl\) rather than being covaried with it because it is the temporally antecedent form of the same construct. The rest should be common sense but anyone legitimately struggling with them can send me an email or read the VES documentation.

The first model assesses merely how correlated the \(g\) factors from the different timepoints are. Initial testing occurred at approximately high school graduation (i.e., ~18 years-of-age) and follow-up occurred around age 35. The same tests administered at both time-points (from the ACB) are residually covaried for goodness-of-fit and the substantive reason that they’re the same tests, and as such expected to share method/test/item variance. Because group factors are not modeled and the battery is large, the fit is poor, so a handful of residual covariances (selected if modification indices met or exceeded 100) are allowed, although the full analysis will include group factors. It is debatable whether other residual covariances should be allowed between early and later subtests; in the model below, four are allowed; they are substantively inconsequential.

GCOR.model<-'

ge =~ ACVE + ACAE + AFQT + GIT + PA

gl =~ WRAT + CVLT + WCST + WBD + WGI + CD + CC + PASAT + GPTL + GPTR + WLGT + ACVL + ACAL

ge ~ gl

ACVE ~~ ACVL
ACAE ~~ ACAL
CD ~~ CC
GPTL ~~ GPTR

ACVE ~~ WRAT 
AFQT ~~ PA + WBD
PA ~~ WBD + CC
WRAT ~~ WGI + WLGT + ACVL
CC ~~ CVLT + WBD
WGI ~~ PASAT'

GCOR.fit <- sem(GCOR.model, data = VES, std.lv=T)

parameterEstimates(GCOR.fit, stand = T) %>% 
  filter(op == "~") %>% 
  dplyr::select(Early = lhs, Later = rhs, "Unstandardized Regression" = est, SE = se, Z = z, 'p-value' = pvalue, Beta = std.all) %>% 
  kable(digits = 3, format = "pandoc")

Early	Later	Unstandardized Regression	SE	Z	p-value	Beta
ge	gl	2.88	0.086	33.7	0	0.945

The stability of the \(g\) factor over time is considerable (\(r = 0.945\)).

A.model <- '
ge =~ ACVE + ACAE + AFQT + GIT + PA

gl =~ WRAT + CVLT + WCST + WBD + WGI + CD + CC + PASAT + GPTL + GPTR + WLGT + ACVL + ACAL

ge ~ gl + education

gl ~ education

ACVE ~~ ACVL
ACAE ~~ ACAL
CD ~~ CC
GPTL ~~ GPTR

ACVE ~~ WRAT 
AFQT ~~ PA + WBD
PA ~~ WBD + CC
WRAT ~~ WGI + WLGT + ACVL
CC ~~ CVLT + WBD
WGI ~~ PASAT'

B.model <- '
ge =~ ACVE + ACAE + AFQT + GIT + PA

gl =~ WRAT + CVLT + WCST + WBD + WGI + CD + CC + PASAT + GPTL + GPTR + WLGT + ACVL + ACAL

gl + education ~ ge

gl + WRAT + CVLT + WCST + WBD + WGI + CD + CC + PASAT + GPTL + GPTR + WLGT + ACVL + ACAL ~ education

ACVE ~~ ACVL
ACAE ~~ ACAL
CD ~~ CC
GPTL ~~ GPTR

ACVE ~~ WRAT 
AFQT ~~ PA + WBD
PA ~~ WBD + CC
WRAT ~~ WGI + WLGT + ACVL
CC ~~ CVLT + WBD
WGI ~~ PASAT'

C.model <- '
ge =~ ACVE + ACAE + AFQT + GIT + PA

gl =~ WRAT + CVLT + WCST + WBD + WGI + CD + CC + PASAT + GPTL + GPTR + WLGT + ACVL + ACAL

gl + education ~ ge

WRAT + CVLT + WCST + WBD + WGI + CD + CC + PASAT + GPTL + GPTR + WLGT + ACVL + ACAL ~ education

ACVE ~~ ACVL
ACAE ~~ ACAL
CD ~~ CC
GPTL ~~ GPTR

ACVE ~~ WRAT 
AFQT ~~ PA + WBD
PA ~~ WBD + CC
WRAT ~~ WGI + WLGT + ACVL
CC ~~ CVLT + WBD
WGI ~~ PASAT'

A.fit <- sem(A.model, data = VES, std.lv = T)
B.fit <- sem(B.model, data = VES, std.lv = T)
C.fit <- sem(C.model, data = VES, std.lv = T)

round(cbind(JUSTG = fitMeasures(A.fit, FITM),
            GPLUS = fitMeasures(B.fit, FITM),
            JUSTS = fitMeasures(C.fit, FITM)), 3)

                    JUSTG      GPLUS      JUSTS
chisq            2456.636   2154.666   2154.667
df                135.000    122.000    123.000
npar               54.000     68.000     67.000
cfi                 0.952      0.958      0.958
rmsea               0.063      0.062      0.062
rmsea.ci.lower      0.061      0.060      0.060
rmsea.ci.upper      0.065      0.064      0.064
aic            595165.045 614311.103 614309.103
bic            595509.054 614744.300 614735.930

It appears that model C - effects on specific subtests - fits best. At this point, it’s appropriate to look at A, B, and C, to see if they show any signs of misfit like Heywood cases.

parameterEstimates(A.fit, stand = T) %>% 
  filter(op == "~") %>% 
  dplyr::select(Dependent = lhs, Independent = rhs, "Unstandardized Regression" = est, SE = se, Z = z, 'p-value' = pvalue, Beta = std.all) %>% 
  kable(digits = 3, format = "pandoc")

Dependent	Independent	Unstandardized Regression	SE	Z	p-value	Beta
ge	gl	2.305	0.071	32.37	0.000	0.934
ge	education	0.026	0.013	2.05	0.041	0.019
gl	education	0.322	0.008	39.29	0.000	0.594

parameterEstimates(B.fit, stand = T) %>% 
  filter(op == "~") %>% 
  dplyr::select(Dependent = lhs, Independent = rhs, "Unstandardized Regression" = est, SE = se, Z = z, 'p-value' = pvalue, Beta = std.all) %>% 
  kable(digits = 3, format = "pandoc")

Dependent	Independent	Unstandardized Regression	SE	Z	p-value	Beta
gl	ge	2.79	NA	NA	NA	0.431
education	ge	1.30	NA	NA	NA	0.566
gl	education	1.91	NA	NA	NA	0.678
WRAT	education	-5.40	NA	NA	NA	-0.833
CVLT	education	-4.31	NA	NA	NA	-0.659
WCST	education	-4.33	NA	NA	NA	-0.663
WBD	education	-6.43	NA	NA	NA	-0.997
WGI	education	-5.33	NA	NA	NA	-0.815
CD	education	-4.35	NA	NA	NA	-0.679
CC	education	-4.97	NA	NA	NA	-0.771
PASAT	education	-6.00	NA	NA	NA	-0.919
GPTL	education	-2.90	NA	NA	NA	-0.450
GPTR	education	-3.23	NA	NA	NA	-0.495
WLGT	education	-3.54	NA	NA	NA	-0.542
ACVL	education	-7.41	NA	NA	NA	-1.140
ACAL	education	-7.62	NA	NA	NA	-1.169

parameterEstimates(C.fit, stand = T) %>% 
  filter(op == "~") %>% 
  dplyr::select(Dependent = lhs, Independent = rhs, "Unstandardized Regression" = est, SE = se, Z = z, 'p-value' = pvalue, Beta = std.all) %>% 
  kable(digits = 3, format = "pandoc")

Dependent	Independent	Unstandardized Regression	SE	Z	p-value	Beta
gl	ge	2.789	0.086	32.550	0.000	0.941
education	ge	1.297	0.033	38.773	0.000	0.566
WRAT	education	0.911	0.080	11.442	0.000	0.140
CVLT	education	-0.137	0.112	-1.228	0.219	-0.021
WCST	education	0.026	0.110	0.233	0.816	0.004
WBD	education	-0.269	0.088	-3.041	0.002	-0.042
WGI	education	1.284	0.082	15.716	0.000	0.196
CD	education	-0.142	0.109	-1.302	0.193	-0.022
CC	education	-0.247	0.106	-2.320	0.020	-0.038
PASAT	education	-0.120	0.103	-1.169	0.242	-0.018
GPTL	education	-0.012	0.114	-0.109	0.913	-0.002
GPTR	education	-0.189	0.115	-1.636	0.102	-0.029
WLGT	education	0.725	0.106	6.827	0.000	0.111
ACVL	education	0.343	0.067	5.086	0.000	0.053
ACAL	education	0.208	0.072	2.869	0.004	0.032

B is unidentified, which is unsurprising given that \(ge\) is a latent variable. The residualized models section (which includes models with \(ge\) as an observed variable/everything residualized for \(ge\), etc.) includes assessments of what this does. It does not change the qualitative result. It appears that B generates two negative ultra-Heywood cases with ACVL and ACAL, so I’ll prune these. Moreover, I’ll prune the insignificant paths in C. There is nothing to change about the A result (yet).

It appears that model B is misspecified thanks to the appearance of an ultra-Heywood case between education and \(gl\), where it’s extremely negatively associated with it (\(\Beta = -1.239\)). The only model making any sense at present is model A, since all of its paths remained significant. However, at this sample size, \(p = 0.04\) is actually evidence against a significant effect of \(ge\) on educational attainment due to Lindley’s paradox. With these caveats in mind, I will prune the insignificant paths to subtests for B and C and retain the path to \(gl\) from education in B.

B2.model <- '
ge =~ ACVE + ACAE + AFQT + GIT + PA

gl =~ WRAT + CVLT + WCST + WBD + WGI + CD + CC + PASAT + GPTL + GPTR + WLGT + ACVL + ACAL

gl + education ~ ge

gl + WRAT + CVLT + WCST + WBD + WGI + CD + CC + PASAT + GPTL + GPTR + WLGT ~ education

ACVE ~~ ACVL
ACAE ~~ ACAL
CD ~~ CC
GPTL ~~ GPTR

ACVE ~~ WRAT 
AFQT ~~ PA + WBD
PA ~~ WBD + CC
WRAT ~~ WGI + WLGT + ACVL
CC ~~ CVLT + WBD
WGI ~~ PASAT'

C2.model <- '
ge =~ ACVE + ACAE + AFQT + GIT + PA

gl =~ WRAT + CVLT + WCST + WBD + WGI + CD + CC + PASAT + GPTL + GPTR + WLGT + ACVL + ACAL

gl + education ~ ge

WRAT + WBD + WGI + CC + WLGT + ACVL + ACAL ~ education

ACVE ~~ ACVL
ACAE ~~ ACAL
CD ~~ CC
GPTL ~~ GPTR

ACVE ~~ WRAT 
AFQT ~~ PA + WBD
PA ~~ WBD + CC
WRAT ~~ WGI + WLGT + ACVL
CC ~~ CVLT + WBD
WGI ~~ PASAT'

B2.fit <- sem(B2.model, data = VES, std.lv = T)
C2.fit <- sem(C2.model, data = VES, std.lv = T)

round(cbind(JUSTG = fitMeasures(A.fit, FITM),
            GPLUS2 = fitMeasures(B2.fit, FITM),
            JUSTS2 = fitMeasures(C2.fit, FITM)), 3)

                    JUSTG     GPLUS2     JUSTS2
chisq            2456.636   2157.257   2162.577
df                135.000    124.000    129.000
npar               54.000     66.000     61.000
cfi                 0.952      0.958      0.958
rmsea               0.063      0.062      0.060
rmsea.ci.lower      0.061      0.059      0.058
rmsea.ci.upper      0.065      0.064      0.063
aic            595165.045 614309.694 614305.013
bic            595509.054 614730.150 614693.617

Model C still fits best between it and B, but more may need pruned.

parameterEstimates(B2.fit, stand = T) %>% 
  filter(op == "~") %>% 
  dplyr::select(Dependent = lhs, Independent = rhs, "Unstandardized Regression" = est, SE = se, Z = z, 'p-value' = pvalue, Beta = std.all) %>% 
  kable(digits = 3, format = "pandoc")

Dependent	Independent	Unstandardized Regression	SE	Z	p-value	Beta
gl	ge	2.789	0.086	32.56	0.000	0.914
education	ge	1.296	0.033	38.75	0.000	0.565
gl	education	0.069	0.014	5.03	0.000	0.052
WRAT	education	0.661	0.077	8.60	0.000	0.102
CVLT	education	-0.287	0.115	-2.50	0.013	-0.044
WCST	education	-0.131	0.113	-1.16	0.248	-0.020
WBD	education	-0.491	0.093	-5.30	0.000	-0.076
WGI	education	1.045	0.084	12.45	0.000	0.160
CD	education	-0.293	0.112	-2.61	0.009	-0.046
CC	education	-0.416	0.110	-3.79	0.000	-0.065
PASAT	education	-0.330	0.106	-3.12	0.002	-0.051
GPTL	education	-0.116	0.117	-0.99	0.322	-0.018
GPTR	education	-0.298	0.119	-2.50	0.012	-0.046
WLGT	education	0.570	0.109	5.21	0.000	0.087

parameterEstimates(C2.fit, stand = T) %>% 
  filter(op == "~") %>% 
  dplyr::select(Dependent = lhs, Independent = rhs, "Unstandardized Regression" = est, SE = se, Z = z, 'p-value' = pvalue, Beta = std.all) %>% 
  kable(digits = 3, format = "pandoc")

Dependent	Independent	Unstandardized Regression	SE	Z	p-value	Beta
gl	ge	2.768	0.084	33.13	0.000	0.940
education	ge	1.291	0.033	38.74	0.000	0.563
WRAT	education	0.933	0.078	11.91	0.000	0.144
WBD	education	-0.245	0.087	-2.81	0.005	-0.038
WGI	education	1.303	0.080	16.31	0.000	0.199
CC	education	-0.169	0.097	-1.75	0.081	-0.026
WLGT	education	0.747	0.105	7.11	0.000	0.114
ACVL	education	0.372	0.065	5.70	0.000	0.057
ACAL	education	0.240	0.070	3.43	0.001	0.037

Some parameters need to be pruned from models B2 and C2.

Model B still generates a negative ultra-Heywood case with \(gl\), but this time it can’t be reasonably attributed to its other indicators. Model C has generated yet another insignificant effect of education on a subtest, so I’ll prune that too.

B3.model <- '
ge =~ ACVE + ACAE + AFQT + GIT + PA

gl =~ WRAT + CVLT + WCST + WBD + WGI + CD + CC + PASAT + GPTL + GPTR + WLGT + ACVL + ACAL

gl + education ~ ge

gl + WRAT + CVLT + WBD + WGI + CD + CC + PASAT + GPTR + WLGT ~ education

ACVE ~~ ACVL
ACAE ~~ ACAL
CD ~~ CC
GPTL ~~ GPTR

ACVE ~~ WRAT 
AFQT ~~ PA + WBD
PA ~~ WBD + CC
WRAT ~~ WGI + WLGT + ACVL
CC ~~ CVLT + WBD
WGI ~~ PASAT'

C3.model <- '
ge =~ ACVE + ACAE + AFQT + GIT + PA

gl =~ WRAT + CVLT + WCST + WBD + WGI + CD + CC + PASAT + GPTL + GPTR + WLGT + ACVL + ACAL

gl + education ~ ge

WRAT + WBD + WGI + WLGT + ACVL + ACAL ~ education

ACVE ~~ ACVL
ACAE ~~ ACAL
CD ~~ CC
GPTL ~~ GPTR

ACVE ~~ WRAT 
AFQT ~~ PA + WBD
PA ~~ WBD + CC
WRAT ~~ WGI + WLGT + ACVL
CC ~~ CVLT + WBD
WGI ~~ PASAT'

B3.fit <- sem(B3.model, data = VES, std.lv = T)
C3.fit <- sem(C3.model, data = VES, std.lv = T)

round(cbind(JUSTG = fitMeasures(A.fit, FITM),
            GPLUS2 = fitMeasures(B3.fit, FITM),
            JUSTS3 = fitMeasures(C3.fit, FITM)), 3)

                    JUSTG     GPLUS2     JUSTS3
chisq            2456.636   2159.496   2165.639
df                135.000    126.000    130.000
npar               54.000     64.000     60.000
cfi                 0.952      0.958      0.958
rmsea               0.063      0.061      0.060
rmsea.ci.lower      0.061      0.059      0.058
rmsea.ci.upper      0.065      0.063      0.062
aic            595165.045 614307.933 614306.075
bic            595509.054 614715.648 614688.308

parameterEstimates(B3.fit, stand = T) %>% 
  filter(op == "~") %>% 
  dplyr::select(Dependent = lhs, Independent = rhs, "Unstandardized Regression" = est, SE = se, Z = z, 'p-value' = pvalue, Beta = std.all) %>% 
  kable(digits = 3, format = "pandoc")

Dependent	Independent	Unstandardized Regression	SE	Z	p-value	Beta
gl	ge	2.790	0.086	32.56	0.000	0.916
education	ge	1.295	0.033	38.71	0.000	0.565
gl	education	0.066	0.014	4.87	0.000	0.050
WRAT	education	0.668	0.077	8.72	0.000	0.103
CVLT	education	-0.277	0.114	-2.42	0.015	-0.042
WBD	education	-0.479	0.092	-5.20	0.000	-0.074
WGI	education	1.059	0.083	12.75	0.000	0.162
CD	education	-0.282	0.112	-2.52	0.012	-0.044
CC	education	-0.405	0.109	-3.71	0.000	-0.063
PASAT	education	-0.316	0.105	-3.01	0.003	-0.048
GPTR	education	-0.224	0.094	-2.38	0.018	-0.034
WLGT	education	0.580	0.109	5.32	0.000	0.089

parameterEstimates(C3.fit, stand = T) %>% 
  filter(op == "~") %>% 
  dplyr::select(Dependent = lhs, Independent = rhs, "Unstandardized Regression" = est, SE = se, Z = z, 'p-value' = pvalue, Beta = std.all) %>% 
  kable(digits = 3, format = "pandoc")

Dependent	Independent	Unstandardized Regression	SE	Z	p-value	Beta
gl	ge	2.760	0.083	33.21	0.000	0.940
education	ge	1.290	0.033	38.71	0.000	0.563
WRAT	education	0.939	0.078	11.99	0.000	0.145
WBD	education	-0.213	0.085	-2.49	0.013	-0.033
WGI	education	1.312	0.080	16.47	0.000	0.201
WLGT	education	0.752	0.105	7.16	0.000	0.115
ACVL	education	0.379	0.065	5.83	0.000	0.058
ACAL	education	0.248	0.070	3.56	0.000	0.038

No more paths require pruning by conventional standards, but, Lindley’s paradox is likely to affect the results since the sample is large. As such, it’s appropriate to prune based on a scaled critical value. For this, I use Naaman’s (2016) “almost-sure” scaling. Now, parameters with a \(Z\)-score below 4.09 are pruned. This means that for model A, the path from \(ge\) to educational attainment is removed.

NZ(4318, 2)

[1] 4.09

A2.model <- '
ge =~ ACVE + ACAE + AFQT + GIT + PA

gl =~ WRAT + CVLT + WCST + WBD + WGI + CD + CC + PASAT + GPTL + GPTR + WLGT + ACVL + ACAL

ge ~ gl

gl ~ education

ACVE ~~ ACVL
ACAE ~~ ACAL
CD ~~ CC
GPTL ~~ GPTR

ACVE ~~ WRAT 
AFQT ~~ PA + WBD
PA ~~ WBD + CC
WRAT ~~ WGI + WLGT + ACVL
CC ~~ CVLT + WBD
WGI ~~ PASAT'

B4.model <- '
ge =~ ACVE + ACAE + AFQT + GIT + PA

gl =~ WRAT + CVLT + WCST + WBD + WGI + CD + CC + PASAT + GPTL + GPTR + WLGT + ACVL + ACAL

gl + education ~ ge

gl + WRAT + WBD + WGI + WLGT ~ education

ACVE ~~ ACVL
ACAE ~~ ACAL
CD ~~ CC
GPTL ~~ GPTR

ACVE ~~ WRAT 
AFQT ~~ PA + WBD
PA ~~ WBD + CC
WRAT ~~ WGI + WLGT + ACVL
CC ~~ CVLT + WBD
WGI ~~ PASAT'

C4.model <- '
ge =~ ACVE + ACAE + AFQT + GIT + PA

gl =~ WRAT + CVLT + WCST + WBD + WGI + CD + CC + PASAT + GPTL + GPTR + WLGT + ACVL + ACAL

gl + education ~ ge

WRAT + WGI + WLGT + ACVL ~ education

ACVE ~~ ACVL
ACAE ~~ ACAL
CD ~~ CC
GPTL ~~ GPTR

ACVE ~~ WRAT 
AFQT ~~ PA + WBD
PA ~~ WBD + CC
WRAT ~~ WGI + WLGT + ACVL
CC ~~ CVLT + WBD
WGI ~~ PASAT'

A2.fit <- sem(A2.model, data = VES, std.lv = T)
B4.fit <- sem(B4.model, data = VES, std.lv = T)
C4.fit <- sem(C4.model, data = VES, std.lv = T)

round(cbind(JUSTG2 = fitMeasures(A2.fit, FITM),
            GPLUS4 = fitMeasures(B4.fit, FITM),
            JUSTS4 = fitMeasures(C4.fit, FITM)), 3)

                   JUSTG2     GPLUS4     JUSTS4
chisq            2460.669   2190.069   2186.998
df                136.000    131.000    132.000
npar               53.000     59.000     58.000
cfi                 0.952      0.957      0.957
rmsea               0.063      0.060      0.060
rmsea.ci.lower      0.061      0.058      0.058
rmsea.ci.upper      0.065      0.063      0.062
aic            595167.078 614328.505 614323.434
bic            595504.717 614704.367 614692.926

Models B and C are roughly tied in terms of absolute fit, but C is more parsimonious with not-significantly-worse exact fit.

parameterEstimates(A2.fit, stand = T) %>% 
  filter(op == "~") %>% 
  dplyr::select(Dependent = lhs, Independent = rhs, "Unstandardized Regression" = est, SE = se, Z = z, 'p-value' = pvalue, Beta = std.all) %>% 
  kable(digits = 3, format = "pandoc")

Dependent	Independent	Unstandardized Regression	SE	Z	p-value	Beta
ge	gl	2.344	0.070	33.4	0	0.946
gl	education	0.325	0.008	40.1	0	0.597

parameterEstimates(B4.fit, stand = T) %>% 
  filter(op == "~") %>% 
  dplyr::select(Dependent = lhs, Independent = rhs, "Unstandardized Regression" = est, SE = se, Z = z, 'p-value' = pvalue, Beta = std.all) %>% 
  kable(digits = 3, format = "pandoc")

Dependent	Independent	Unstandardized Regression	SE	Z	Beta
gl	ge	2.795	0.086	32.55	0.923
education	ge	1.287	0.033	38.48	0.561
gl	education	0.049	0.013	3.75	0.037
WRAT	education	0.705	0.075	9.35	0.109
WBD	education	-0.377	0.088	-4.26	-0.058
WGI	education	1.099	0.081	13.55	0.168
WLGT	education	0.634	0.107	5.91	0.097

parameterEstimates(C4.fit, stand = T) %>% 
  filter(op == "~") %>% 
  dplyr::select(Dependent = lhs, Independent = rhs, "Unstandardized Regression" = est, SE = se, Z = z, 'p-value' = pvalue, Beta = std.all) %>% 
  kable(digits = 3, format = "pandoc")

Dependent	Independent	Unstandardized Regression	SE	Z	Beta
gl	ge	2.793	0.083	33.77	0.941
education	ge	1.296	0.033	38.97	0.565
WRAT	education	0.908	0.077	11.76	0.140
WGI	education	1.264	0.078	16.29	0.194
WLGT	education	0.727	0.105	6.93	0.111
ACVL	education	0.339	0.063	5.40	0.052

It appears that model B is not viable for its theoretical backing, since it requires dropping the path from education to \(gl\).

B5.model <- '
ge =~ ACVE + ACAE + AFQT + GIT + PA

gl =~ WRAT + CVLT + WCST + WBD + WGI + CD + CC + PASAT + GPTL + GPTR + WLGT + ACVL + ACAL

gl + education ~ ge

WRAT + WBD + WGI + WLGT ~ education

ACVE ~~ ACVL
ACAE ~~ ACAL
CD ~~ CC
GPTL ~~ GPTR

ACVE ~~ WRAT 
AFQT ~~ PA + WBD
PA ~~ WBD + CC
WRAT ~~ WGI + WLGT + ACVL
CC ~~ CVLT + WBD
WGI ~~ PASAT'

B5.fit <- sem(B5.model, data = VES, std.lv = T)

round(cbind(JUSTG2 = fitMeasures(A2.fit, FITM),
            GPLUS5 = fitMeasures(B5.fit, FITM),
            JUSTS4 = fitMeasures(C4.fit, FITM)), 3)

                   JUSTG2     GPLUS5     JUSTS4
chisq            2460.669   2203.909   2186.998
df                136.000    132.000    132.000
npar               53.000     58.000     58.000
cfi                 0.952      0.957      0.957
rmsea               0.063      0.060      0.060
rmsea.ci.lower      0.061      0.058      0.058
rmsea.ci.upper      0.065      0.063      0.062
aic            595167.078 614340.345 614323.434
bic            595504.717 614709.837 614692.926

parameterEstimates(B5.fit, stand = T) %>% 
  filter(op == "~") %>% 
  dplyr::select(Dependent = lhs, Independent = rhs, "Unstandardized Regression" = est, SE = se, Z = z, 'p-value' = pvalue, Beta = std.all) %>% 
  kable(digits = 3, format = "pandoc")

Dependent	Independent	Unstandardized Regression	SE	Z	Beta
gl	ge	2.877	0.086	33.52	0.945
education	ge	1.300	0.033	39.13	0.567
WRAT	education	0.756	0.073	10.39	0.117
WBD	education	-0.302	0.084	-3.57	-0.047
WGI	education	1.175	0.077	15.29	0.180
WLGT	education	0.677	0.105	6.46	0.104

Now, to perform an ANOVA on the final models. I will include B5 in the first ANOVA even though it’s a theoretically bankrupt model. In the second ANOVA, I’ll include B4 instead of B5, since it was the last theoretically admissible model, even though the path it offered was not significant. The important test here is the \(\chi^2\) exact-fit test. The same scaling will be used here.

anova(A2.fit, B5.fit, C4.fit)

Chi-Squared Difference Test

        Df    AIC    BIC Chisq Chisq diff Df diff          Pr(>Chisq)    
B5.fit 132 614340 614710  2204                                           
C4.fit 132 614323 614693  2187      -16.9       0                        
A2.fit 136 595167 595505  2461      273.7       4 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(A2.fit, B4.fit, C4.fit)

Chi-Squared Difference Test

        Df    AIC    BIC Chisq Chisq diff Df diff          Pr(>Chisq)    
B4.fit 131 614329 614704  2190                                           
C4.fit 132 614323 614693  2187       -3.1       1                   1    
A2.fit 136 595167 595505  2461      273.7       4 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A shift of 273.7 \(\chi^2\) is extremely significant, making model A2 less parsimonious than C4, which is the most parsimonious, leaving C4 the best fitting model. Since these models went through multiple stages of pruning, the use of AIC/BIC may not be clear. Anyone who disagrees should email me why. As I see it, these models were not theoretically preselected and the model with extremely high \(\chi^2\) relative to the rest - thanks to an implausible assumption, perhaps - should have the best AIC/BIC, despite failing to recover the data well. I would certainly like to see more data on this if anyone has it. Either way, the models to be contrasted are models of specific effects and effects directly on \(gl\). One way to make ground here is to conduct a power analysis, which I do below, but this may be superfluous since the rejection of the A2 model was so incredibly severe and my criteria were already adjusted to be quite stringent, allowing it a great deal of leeway.

In the accepted model, C, the effect seems to be relegated to the WRAT, WGI, WLGT, and ACVL subtests. It is clear based on the direction of effect for WBD, that the pruning might be helpful insofar as it removes an implausible effect (i.e., education worsening spatial ability). However, it is possible that there is attrition from education by level of spatial ability (to-do: cite the study on spatial ability and lower academic achievement James likes so much, and the articles in NBER/AEA on association of spatial ability and market returns; or, explain that this could be that spatial is just bad for achievement) or that education really does have a negative effect. As I show even further below, the effect of education on these subtests is in no small part to generate measurement variance.

EDUVESLATS <- list(
  ge = early_tests,
  gl = later_tests,
  education = c("WRAT", "WGI", "WLGT", "ACVL"))

semPaths(C4.fit, "model", "std", title = F, residuals = F, groups = "EDUVESLATS", pastel = T, mar = c(4, 1, 3, 1), intercepts = F, layout = "tree", bifactor = c("gl", "ge"), exoCov = F, curvePivot = T)

Power Analysis

Do Later

Residualization Path Model

Do Later

Ethnic Confounding

Do Later

MGCFA

Do Later

Indicator Means

Do Later

Intercepts

Do Later

Residual Variances

Do Later

LSEM

Do Later

MNLFA

Do Later

IQ Gains

Do Later

Naive and Proper g Effects?

Do Later

Income Gains

Do Later

Mediation by IQ?

Do Later

Mediation by Intelligence?

Do Later

Did it Affect Nerve Conduction Velocity? (Or Any Other Sensory Measure?)

Do Later

To-Do

Make a ‘To-Do’ list
Note Emil’s earlier investigation of the same dataset included both iterations of the ACB subtests in the earlier tests rather than in their proper locations, explaining his result
Goldberger and using gated/ungated \(p-values\)
Show lower initial approximate fit threshold models
Show triviality of overfitting a model for information theoretic fits

References

To-Do

Gunn, H. J., Grimm, K. J., & Edwards, M. C. (2019). Evaluation of Six Effect Size Measures of Measurement Non-Invariance for Continuous Outcomes. Structural Equation Modeling: A Multidisciplinary Journal, 0(0), 1-12. https://doi.org/10.1080/10705511.2019.1689507

Naaman, M. (2016). Almost sure hypothesis testing and a resolution of the Jeffreys-Lindley paradox. Electronic Journal of Statistics, 10(1), 1526-1550. https://doi.org/10.1214/16-EJS1146

Millsap, R. E., & Olivera-Aguilar, M. (2012). Investigating measurement invariance using confirmatory factor analysis. In R. H. Hoyle, Handbook of structural equation modeling (pp. 380-392). New York, NY: Guildford Press

The Effect of Education On Intelligence in the Vietnam Experience Study