# ============================================================
# FFCWS: Admixture Regression Analysis
# African Ancestry Sample (Both-Black + BW Parents)
# City-Weighted, HC1 Robust SEs
# ============================================================
library(haven)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(sandwich)
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(mice)
##
## Attaching package: 'mice'
## The following object is masked from 'package:stats':
##
## filter
## The following objects are masked from 'package:base':
##
## cbind, rbind
library(sirt)
## Warning: package 'sirt' was built under R version 4.5.3
## - sirt 4.2-133 (2025-09-27 12:57:51)
library(lavaan)
## This is lavaan 0.6-19
## lavaan is FREE software! Please report any bugs.
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.2
library(mediation)
## Warning: package 'mediation' was built under R version 4.5.2
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
## Loading required package: Matrix
## Loading required package: mvtnorm
## mediation: Causal Mediation Analysis
## Version: 4.5.1
library(mgcv)
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
## This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.
##
## Attaching package: 'mgcv'
## The following object is masked from 'package:sirt':
##
## rmvn
# ==============================================================
# 1. Load Data
# ==============================================================
d <- read_sav("C:/Users/mh198/Documents/Data/FFCWS/DS0001/31622-0001-Data.sav")
PGS <- read_sav("C:/Users/mh198/Documents/Data/FFCWS//DS0007/31622-0007-Data.sav")
d <- merge(d, PGS[, c("IDNUM", "K5ADMIXA", "K5PC1A", "K5PC2A",
"K5PC3A", "K5PC4A", "K5PC5A")],
by = "IDNUM", all.x = TRUE)
# ==============================================================
# 2. African Ancestry PGS Sample
# ==============================================================
d <- d[d$K5ADMIXA %in% c(0, 1), ]
# ==============================================================
# 3. Black Parent Identification
# ==============================================================
d$father_black_nh <- d$CF1ETHRACE == 2
d$mother_black_nh <- d$CM1ETHRACE == 2
d$father_black_hisp <- (d$CF1ETHRACE == 3) & (d$F1H3 == 2)
d$mother_black_hisp <- (d$CM1ETHRACE == 3) & (d$M1H3 == 2)
d$father_black_hisp[is.na(d$father_black_hisp)] <- FALSE
d$mother_black_hisp[is.na(d$mother_black_hisp)] <- FALSE
d$father_any_black <- d$father_black_nh | d$father_black_hisp
d$mother_any_black <- d$mother_black_nh | d$mother_black_hisp
d$any_black_parent <- d$father_any_black | d$mother_any_black
# ==============================================================
# 4. Continental African Exclusion
# ==============================================================
d$mother_african_born <- (d$M1H2 == 2) & (d$M1H2B == 1)
d$father_african_born <- (d$F1H2 == 2) & (d$F1H2B == 1)
d$mother_african_born[is.na(d$mother_african_born)] <- FALSE
d$father_african_born[is.na(d$father_african_born)] <- FALSE
d$gp_african <- (d$F3H1A == 101) | (d$F3H1B == 101) |
(d$M3H1A == 101) | (d$M3H1B == 101)
d$gp_african[is.na(d$gp_african)] <- FALSE
d$african_origin <- d$mother_african_born | d$father_african_born | d$gp_african
# ==============================================================
# 5. Generation & SIRE Dummies
# ==============================================================
d$parent_fb <- (d$M1H2 == 2) | (d$F1H2 == 2)
d$parent_fb[is.na(d$parent_fb)] <- FALSE
d$gp_fb <- (d$F3H1A >= 101) | (d$F3H1B >= 101) |
(d$M3H1A >= 101) | (d$M3H1B >= 101)
d$gp_fb[is.na(d$gp_fb)] <- FALSE
d$gen_2nd <- d$parent_fb & !d$african_origin
d$gen_3rd <- (!d$parent_fb) & d$gp_fb & !d$african_origin
d$sire <- NA_character_
d$sire[d$father_any_black & d$mother_any_black] <- "BB"
d$sire[(d$father_any_black & d$CM1ETHRACE == 1) |
(d$CF1ETHRACE == 1 & d$mother_any_black)] <- "BW"
d$sire[(d$father_any_black & d$CM1ETHRACE == 3 & !d$mother_black_hisp) |
(d$CF1ETHRACE == 3 & !d$father_black_hisp & d$mother_any_black)] <- "BH"
d$sire[d$any_black_parent & is.na(d$sire)] <- "BO"
d$sire[!d$any_black_parent] <- "None"
d$sire <- factor(d$sire, levels = c("BB", "BW", "BH", "BO", "None"))
# ==============================================================
# 6. Child TVIP Exclusion Flag
# ==============================================================
d$ch_tvip3_score <- ifelse(d$CH3TVIPSTD >= 0, d$CH3TVIPSTD, NA)
d$took_tvip3 <- as.numeric(!is.na(d$ch_tvip3_score))
cat("Child TVIP at age 3:", sum(d$took_tvip3), "cases flagged for exclusion\n")
## Child TVIP at age 3: 1 cases flagged for exclusion
# ==============================================================
# 7. Demographics & SES (3-component)
# ==============================================================
d$sex <- ifelse(d$CM1BSEX %in% c(1, 2), d$CM1BSEX - 1, NA)
d$age <- ifelse(d$CH5AGEM >= 0, d$CH5AGEM, NA)
d$single <- ifelse(d$CM1MARF == 0 & d$CM1COHF == 0, 1,
ifelse(d$CM1MARF %in% c(0, 1) & d$CM1COHF %in% c(0, 1), 0, NA))
d$M_edu <- NA
for (col in c("CM5EDU", "CM4EDU", "CM3EDU", "CM2EDU", "CM1EDU")) {
d$M_edu <- ifelse(is.na(d$M_edu) & d[[col]] >= 1, d[[col]], d$M_edu)
}
d$F_edu <- NA
for (col in c("CF5EDU", "CF4EDU", "CF3EDU", "CF2EDU", "CF1EDU")) {
d$F_edu <- ifelse(is.na(d$F_edu) & d[[col]] >= 1, d[[col]], d$F_edu)
}
d$mid_edu <- rowMeans(cbind(d$F_edu, d$M_edu), na.rm = TRUE)
d$mid_edu[is.na(d$F_edu) & is.na(d$M_edu)] <- NA
d$pov_bl <- ifelse(d$CM1INPOV >= 0, d$CM1INPOV, NA)
d$pov_yr1 <- ifelse(d$CM2POVCO >= 0, d$CM2POVCO, NA)
d$pov_yr3 <- ifelse(d$CM3POVCO >= 0, d$CM3POVCO, NA)
d$pov_yr5 <- ifelse(d$CM4POVCO >= 0, d$CM4POVCO, NA)
d$pov_yr9 <- ifelse(d$CM5POVCO >= 0, d$CM5POVCO, NA)
d$ave_pov <- rowMeans(cbind(d$pov_bl, d$pov_yr1, d$pov_yr3, d$pov_yr5, d$pov_yr9), na.rm = TRUE)
d$ave_pov[is.na(d$pov_bl) & is.na(d$pov_yr1) & is.na(d$pov_yr3) & is.na(d$pov_yr5) & is.na(d$pov_yr9)] <- NA
inc_vars <- c("CM1HHINC", "CM2HHINC", "CM3HHINC", "CM4HHINC", "CM5HHINC")
inc_names <- c("inc_bl", "inc_yr1", "inc_yr3", "inc_yr5", "inc_yr9")
for (i in seq_along(inc_vars)) {
d[[inc_names[i]]] <- ifelse(d[[inc_vars[i]]] >= 0, d[[inc_vars[i]]], NA)
d[[paste0("log_", inc_names[i])]] <- log(pmax(d[[inc_names[i]]], 1))
}
d$ave_log_inc <- rowMeans(d[, paste0("log_", inc_names)], na.rm = TRUE)
d$ave_log_inc[rowSums(!is.na(d[, inc_names])) == 0] <- NA
ses_indicators <- data.frame(
mid_edu = as.vector(scale(d$mid_edu)),
ave_pov = as.vector(scale(d$ave_pov)),
ave_log_inc = as.vector(scale(d$ave_log_inc))
)
n_ses_valid <- rowSums(!is.na(ses_indicators))
set.seed(54321)
ses_imp <- mice(ses_indicators, m = 5, method = "pmm", maxit = 10, printFlag = FALSE)
ses_imputed <- complete(ses_imp, action = 1)
d$SES_raw <- rowMeans(ses_imputed, na.rm = FALSE)
d$SES_raw[n_ses_valid == 0] <- NA
d$wt <- ifelse(d$K5CITYWT > 0, d$K5CITYWT, NA)
# ==============================================================
# 8. Child Cognitive Variables: g (GAM-based residualization)
# ==============================================================
d$PPVT <- ifelse(d$CH5PPVTSS >= 0, d$CH5PPVTSS, NA)
d$WJ9 <- ifelse(d$CH5WJ9SS >= 0, d$CH5WJ9SS, NA)
d$WJ10 <- ifelse(d$CH5WJ10SS >= 0, d$CH5WJ10SS, NA)
d$DS <- ifelse(d$CH5DSSS >= 0, d$CH5DSSS, NA)
# GAM-based residualization: smooth age + parametric sex
resid_fun_gam <- function(y, data) {
fit <- gam(y ~ s(age, k = 5) + sex, data = data,
na.action = na.exclude, method = "REML")
as.vector(residuals(fit))
}
d$PPVT_r <- resid_fun_gam(d$PPVT, d)
d$WJ9_r <- resid_fun_gam(d$WJ9, d)
d$WJ10_r <- resid_fun_gam(d$WJ10, d)
d$DS_r <- resid_fun_gam(d$DS, d)
d$PPVT_z <- as.vector(scale(d$PPVT_r))
d$WJ9_z <- as.vector(scale(d$WJ9_r))
d$WJ10_z <- as.vector(scale(d$WJ10_r))
d$DS_z <- as.vector(scale(d$DS_r))
# 3-test g (DS kept for LSEM/CFA only)
cog_imp_data <- data.frame(
PPVT_z = d$PPVT_z,
WJ9_z = d$WJ9_z,
WJ10_z = d$WJ10_z
)
has_any <- rowSums(!is.na(cog_imp_data)) >= 1
set.seed(12345)
imp <- mice(cog_imp_data, m = 5, method = "norm", maxit = 10, printFlag = FALSE)
cog_imputed <- complete(imp, action = 1)
cog_for_pca <- as.matrix(cog_imputed[has_any, ])
pca_out <- prcomp(cog_for_pca, center = FALSE, scale. = FALSE)
d$g_raw <- NA
d$g_raw[has_any] <- pca_out$x[, 1]
cat("PCA loadings (3-test g):", round(pca_out$rotation[, 1], 3), "\n")
## PCA loadings (3-test g): 0.564 0.588 0.58
cat("Variance explained:", round(summary(pca_out)$importance[2, 1] * 100, 1), "%\n")
## Variance explained: 71.7 %
# ==============================================================
# 9. Parental WAIS: Conservative USB + English Filter
# ==============================================================
d$mother_wais_raw <- ifelse(d$CM3COGSC >= 0 & d$CM3COGSC <= 15, d$CM3COGSC, NA)
d$father_wais_raw <- ifelse(d$CF3COGSC >= 0 & d$CF3COGSC <= 15, d$CF3COGSC, NA)
d$mother_confirmed_fb <- as.numeric(d$M1H2B %in% c(1, 2, 3, 4, 5))
d$father_confirmed_fb <- as.numeric(d$F1H2B %in% c(1, 2, 3, 4, 5))
d$mother_confirmed_span <- as.numeric(d$CM3SPAN == 1)
d$father_confirmed_span <- as.numeric(d$CF3SPAN == 1)
d$mother_wais <- ifelse(d$mother_confirmed_fb == 1 | d$mother_confirmed_span == 1,
NA, d$mother_wais_raw)
d$father_wais <- ifelse(d$father_confirmed_fb == 1 | d$father_confirmed_span == 1,
NA, d$father_wais_raw)
cat("\nMother WAIS: total =", sum(!is.na(d$mother_wais_raw)),
", after filter =", sum(!is.na(d$mother_wais)),
", dropped =", sum(!is.na(d$mother_wais_raw)) - sum(!is.na(d$mother_wais)), "\n")
##
## Mother WAIS: total = 1536 , after filter = 1480 , dropped = 56
cat("Father WAIS: total =", sum(!is.na(d$father_wais_raw)),
", after filter =", sum(!is.na(d$father_wais)),
", dropped =", sum(!is.na(d$father_wais_raw)) - sum(!is.na(d$father_wais)), "\n")
## Father WAIS: total = 1276 , after filter = 1231 , dropped = 45
# Age-residualize
d$mother_age <- ifelse(d$CM3AGE > 0, d$CM3AGE, NA)
d$father_age <- ifelse(d$CF3AGE > 0, d$CF3AGE, NA)
m_mask <- !is.na(d$mother_wais) & !is.na(d$mother_age)
if (sum(m_mask) > 10) {
d$mother_age_c <- d$mother_age - mean(d$mother_age[m_mask])
d$mother_age_c2 <- d$mother_age_c^2
d$M_WAIS_res <- NA_real_
d$M_WAIS_res[m_mask] <- residuals(lm(mother_wais ~ mother_age_c + mother_age_c2,
data = d[m_mask, ]))
}
f_mask <- !is.na(d$father_wais) & !is.na(d$father_age)
if (sum(f_mask) > 10) {
d$father_age_c <- d$father_age - mean(d$father_age[f_mask])
d$father_age_c2 <- d$father_age_c^2
d$F_WAIS_res <- NA_real_
d$F_WAIS_res[f_mask] <- residuals(lm(father_wais ~ father_age_c + father_age_c2,
data = d[f_mask, ]))
}
d$ave_WAIS_raw <- rowMeans(cbind(d$M_WAIS_res, d$F_WAIS_res), na.rm = TRUE)
d$ave_WAIS_raw[is.na(d$M_WAIS_res) & is.na(d$F_WAIS_res)] <- NA
cat("Average WAIS (filtered): N =", sum(!is.na(d$ave_WAIS_raw)), "\n")
## Average WAIS (filtered): N = 1533
# ==============================================================
# 10. Build Subsets (after TVIP exclusion)
# ==============================================================
d <- d[d$took_tvip3 == 0, ]
cat("After TVIP exclusion: N =", nrow(d), "\n")
## After TVIP exclusion: N = 1639
d$PC1A_z <- as.vector(scale(d$K5PC1A))
d_apa <- d[d$any_black_parent & !d$african_origin, ]
d_apa$sire <- droplevels(d_apa$sire)
# Parent self-ID flags for individual WAIS models
d_apa$mother_black_sire <- as.numeric(d_apa$CM1ETHRACE == 2)
d_apa$father_black_sire <- as.numeric(d_apa$CF1ETHRACE == 2)
d_apa_noBW <- d_apa[d_apa$sire != "BW", ]
d_apa_noBW$sire <- droplevels(d_apa_noBW$sire)
d_haa <- d[d$CF1ETHRACE == 2 & d$CM1ETHRACE == 2 &
!d$african_origin & !d$parent_fb & !d$gp_fb, ]
cat("APA:", nrow(d_apa), " | APA excl BW:", nrow(d_apa_noBW), " | HAA:", nrow(d_haa), "\n")
## APA: 1588 | APA excl BW: 1516 | HAA: 1323
# ==============================================================
# 11. LSEM Analysis — Liu et al. (2025) approach
# Progressive invariance testing with bootstrap
# ==============================================================
# DS included for model identification (4 indicators, 2 df)
# Run on post-TVIP APA sample to match regression analyses
d_apa_cog <- d_apa[rowSums(!is.na(d_apa[, c("PPVT_z", "WJ9_z", "WJ10_z", "DS_z")])) >= 1, ]
grid_points <- quantile(d_apa_cog$PC1A_z,
probs = c(0.20, 0.35, 0.50, 0.65, 0.80),
na.rm = TRUE)
lsem_model <- 'g =~ PPVT_z + WJ9_z + WJ10_z + DS_z'
# ── Define invariance constraints ────────────────────────────────────────
# Metric (weak) — loadings invariant
weak_constraints <- c("g=~PPVT_z", "g=~WJ9_z", "g=~WJ10_z", "g=~DS_z")
# Scalar (strong) — loadings + intercepts invariant
strong_constraints <- c(weak_constraints,
"PPVT_z~1", "WJ9_z~1", "WJ10_z~1", "DS_z~1")
# Strict — loadings + intercepts + residual variances
strict_constraints <- c(strong_constraints,
"PPVT_z~~PPVT_z", "WJ9_z~~WJ9_z",
"WJ10_z~~WJ10_z", "DS_z~~DS_z")
# ── Configural ───────────────────────────────────────────────────────────
parallel::detectCores()
## [1] 32
lsem_config <- lsem.estimate(
data = as.data.frame(d_apa_cog),
moderator = "PC1A_z",
moderator.grid = grid_points,
lavmodel = lsem_model,
h = 1.5,
estimator = "ML",
meanstructure = TRUE,
standardized = TRUE,
est_joint = TRUE
)
## ** Residualize Data
## ** Fit lavaan model
## |*****|
## |lavoptions ... done.
## lavdata ... done.
## Number of observations per group:
## Group 1 415
## Group 2 673
## Group 3 782
## Group 4 793
## Group 5 707
## lavsamplestats ... done.
## lavh1 ... start:
## lavh1 ... done.
## lavpartable bounds ... done.
## lavstart ... done.
## lavmodel ... done.
## lavoptim ... start:
## attempt 1 -- default options
## objective function = 0.0725269442612360
## quasi-Newton steps using NLMINB:
## objective function = 0.0725269442612360
## objective function = 0.0588563752561805
## objective function = 0.0270056466184860
## objective function = 0.0210308952064980
## objective function = 0.0162708728292977
## objective function = 0.0144724321967550
## objective function = 0.0133186275293610
## objective function = 0.0119570202883823
## objective function = 0.0110561959205937
## objective function = 0.0103766602343360
## objective function = 0.0099292582693531
## objective function = 0.0095664544279759
## objective function = 0.0092385440551466
## objective function = 0.0088935697491971
## objective function = 0.0087137702334482
## objective function = 0.0085931967573048
## objective function = 0.0085234683929610
## objective function = 0.0084667862856638
## objective function = 0.0084228674348790
## objective function = 0.0083879322704697
## objective function = 0.0083587123599363
## objective function = 0.0083324027104931
## objective function = 0.0083069817621604
## objective function = 0.0082826012564447
## objective function = 0.0082609647475788
## objective function = 0.0082456969064643
## objective function = 0.0082345497386336
## objective function = 0.0082248423976411
## objective function = 0.0082161307427895
## objective function = 0.0082088721104753
## objective function = 0.0082039966403310
## objective function = 0.0082014389211640
## objective function = 0.0081989492124753
## objective function = 0.0081966453118561
## objective function = 0.0081942421699808
## objective function = 0.0081926960503104
## objective function = 0.0081916222565688
## objective function = 0.0081908949101970
## objective function = 0.0081902528938547
## objective function = 0.0081896970215332
## objective function = 0.0081890880719255
## objective function = 0.0081886035419813
## objective function = 0.0081881598792567
## objective function = 0.0081878762016064
## objective function = 0.0081876784286030
## objective function = 0.0081875354060321
## objective function = 0.0081874283577126
## objective function = 0.0081873565508840
## objective function = 0.0081873064149768
## objective function = 0.0081872675392588
## objective function = 0.0081872377310486
## objective function = 0.0081872179803860
## objective function = 0.0081872052015698
## objective function = 0.0081871949613283
## objective function = 0.0081871859911021
## objective function = 0.0081871765635997
## objective function = 0.0081871713970859
## objective function = 0.0081871683223214
## objective function = 0.0081871668703660
## objective function = 0.0081871657129795
## objective function = 0.0081871647409269
## objective function = 0.0081871639258151
## objective function = 0.0081871633041530
## objective function = 0.0081871629658704
## objective function = 0.0081871627502336
## objective function = 0.0081871626297596
## objective function = 0.0081871625383501
## objective function = 0.0081871624631940
## objective function = 0.0081871624040176
## objective function = 0.0081871623645527
## objective function = 0.0081871623447807
## objective function = 0.0081871623315614
## objective function = 0.0081871623232390
## objective function = 0.0081871623175860
## objective function = 0.0081871623140465
## objective function = 0.0081871623114891
## objective function = 0.0081871623097353
## objective function = 0.0081871623088481
## objective function = 0.0081871623088481
## convergence status (0=ok): 0
## nlminb message says: relative convergence (4)
## number of iterations: 77
## number of function evaluations [objective, gradient]: 78 78
## lavoptim ... done.
## lavimplied ... done.
## lavloglik ... done.
## computing TEST for test(s) = standard ... done.
## lavbaseline ... done.
## post check ... done.
##
## Use lsem.bootstrap() for computing standard errors!
##
## -|
## ** Parameter summary
summary(lsem_config)
## -----------------------------------------------------------------
## Local Structural Equation Model
##
## sirt 4.2-133 (2025-09-27 12:57:51)
## lavaan 0.6-19 (2024-09-26 15:50:14 UTC)
##
## R version 4.5.1 (2025-06-13 ucrt) x86_64, mingw32 | nodename=HUANG | login=mh198
##
## Function 'sirt::lsem.estimate', type='LSEM'
##
##
## Call:
## lsem.estimate(data = as.data.frame(d_apa_cog), moderator = "PC1A_z",
## moderator.grid = grid_points, lavmodel = lsem_model, h = 1.5,
## standardized = TRUE, est_joint = TRUE, estimator = "ML",
## meanstructure = TRUE)
##
## Date of Analysis: 2026-06-12 14:08:25.85806
## Time difference of 4.666388 secs
## Computation Time: 4.666388
##
## Number of observations in datasets = 1538
## Used observations in analysis = 1524
## Used sampling weights: FALSE
## Bandwidth factor = 1.5
## Bandwidth = 0.294
## Number of focal points for moderator = 5
##
## Used joint estimation: TRUE
## Used sufficient statistics: TRUE
## Used local linear smoothing: TRUE
## Used pseudo weights: FALSE
## Used lavaan package: TRUE
## Used lavaan.survey package: FALSE
##
## Mean structure modelled: TRUE
##
## lavaan Model
## g =~ PPVT_z + WJ9_z + WJ10_z + DS_z
##
## Global Fit Statistics for Joint Estimation
##
## stat value
## 1 rmsea 0.082
## 2 cfi 0.989
## 3 tli 0.967
## 4 gfi 0.992
## 5 srmr 0.020
##
##
## Parameter Estimate Summary
##
## par parindex M SD MAD Min Max lin_int lin_slo
## 1 g=~PPVT_z 1 1.000 0.000 0.000 1.000 1.000 1.000 0.000
## 2 g=~WJ9_z 2 1.315 0.086 0.078 1.211 1.436 1.413 -0.311
## 3 g=~WJ10_z 3 1.257 0.059 0.045 1.205 1.406 1.326 -0.219
## 4 g=~DS_z 4 0.877 0.030 0.023 0.859 0.959 0.910 -0.105
## 5 PPVT_z~~PPVT_z 5 0.568 0.086 0.074 0.476 0.736 0.671 -0.325
## 6 WJ9_z~~WJ9_z 6 0.322 0.023 0.014 0.301 0.386 0.347 -0.081
## 7 WJ10_z~~WJ10_z 7 0.402 0.021 0.019 0.355 0.421 0.400 0.008
## 8 DS_z~~DS_z 8 0.718 0.024 0.019 0.667 0.741 0.691 0.083
## 9 g~~g 9 0.397 0.036 0.030 0.325 0.434 0.354 0.135
## 10 PPVT_z~1 10 -0.001 0.004 0.003 -0.008 0.004 -0.001 -0.001
## 11 WJ9_z~1 11 0.001 0.003 0.002 -0.002 0.008 0.004 -0.007
## 12 WJ10_z~1 12 -0.002 0.005 0.004 -0.007 0.004 0.001 -0.008
## 13 DS_z~1 13 -0.001 0.002 0.002 -0.004 0.001 -0.001 -0.001
## 14 g~1 14 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## 15 std__g=~PPVT_z 15 0.642 0.045 0.039 0.553 0.690 0.588 0.171
## 16 std__g=~WJ9_z 16 0.824 0.010 0.006 0.796 0.832 0.819 0.017
## 17 std__g=~WJ10_z 17 0.780 0.011 0.009 0.769 0.802 0.780 -0.001
## 18 std__g=~DS_z 18 0.545 0.006 0.005 0.539 0.556 0.545 0.000
## 19 std__PPVT_z~~PPVT_z 19 0.586 0.057 0.050 0.523 0.694 0.654 -0.214
## 20 std__WJ9_z~~WJ9_z 20 0.321 0.016 0.009 0.308 0.366 0.329 -0.028
## 21 std__WJ10_z~~WJ10_z 21 0.392 0.017 0.014 0.356 0.409 0.391 0.002
## 22 std__DS_z~~DS_z 22 0.702 0.006 0.006 0.691 0.710 0.703 0.000
## 23 std__g~~g 23 1.000 0.000 0.000 1.000 1.000 1.000 0.000
## 24 std__PPVT_z~1 24 -0.001 0.004 0.003 -0.008 0.003 -0.001 -0.001
## 25 std__WJ9_z~1 25 0.001 0.003 0.002 -0.002 0.008 0.004 -0.007
## 26 std__WJ10_z~1 26 -0.002 0.005 0.004 -0.007 0.004 0.001 -0.008
## 27 std__DS_z~1 27 -0.001 0.002 0.002 -0.004 0.001 -0.001 0.000
## 28 std__g~1 28 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## SD_nonlin
## 1 0.000
## 2 0.026
## 3 0.014
## 4 0.013
## 5 0.012
## 6 0.010
## 7 0.021
## 8 0.009
## 9 0.005
## 10 0.004
## 11 0.003
## 12 0.004
## 13 0.002
## 14 0.000
## 15 0.006
## 16 0.009
## 17 0.011
## 18 0.006
## 19 0.008
## 20 0.014
## 21 0.017
## 22 0.006
## 23 0.000
## 24 0.003
## 25 0.003
## 26 0.004
## 27 0.002
## 28 0.000
##
## Distribution of Moderator: Density and Effective Sample Size
##
## M=0.062 | SD=0.851
##
## moderator wgt Neff
## 1 -0.246 0.100 419.601
## 2 0.070 0.167 679.165
## 3 0.264 0.233 789.229
## 4 0.453 0.259 800.087
## 5 0.619 0.241 713.065
##
## variable M SD min max
## 1 moderator 0.062 0.851 -4.930 1.448
## 2 wgt 0.200 0.066 0.100 0.259
## 3 Neff 680.229 154.315 419.601 800.087
b_config <- lsem.bootstrap(lsem_config, R = 50, n.core = 14)
## Compute LSEM for bootstrap samples:
## Warning in sqrt(w): NaNs produced
# ── Metric (weak) ────────────────────────────────────────────────────────
lsem_weak <- lsem.estimate(
data = as.data.frame(d_apa_cog),
moderator = "PC1A_z",
moderator.grid = grid_points,
lavmodel = lsem_model,
h = 1.5,
estimator = "ML",
meanstructure = TRUE,
standardized = TRUE,
par_invariant = weak_constraints
)
## Use joint estimation
## ** Residualize Data
## ** Fit lavaan model
## |*****|
## |lavoptions ... done.
## lavdata ... done.
## Number of observations per group:
## Group 1 415
## Group 2 673
## Group 3 782
## Group 4 793
## Group 5 707
## lavsamplestats ... done.
## lavh1 ... start:
## lavh1 ... done.
## lavpartable bounds ... done.
## lavstart ... done.
## lavmodel ... done.
## lavoptim ... start:
## attempt 1 -- default options
## objective function = 0.0725269442612357
## quasi-Newton steps using NLMINB:
## objective function = 0.0725269442612357
## objective function = 0.0589336747621080
## objective function = 0.0266784670386390
## objective function = 0.0210113040431718
## objective function = 0.0165662029281500
## objective function = 0.0144912328134312
## objective function = 0.0134192340063573
## objective function = 0.0120333458203883
## objective function = 0.0111935861154176
## objective function = 0.0105554481992239
## objective function = 0.0101974403139898
## objective function = 0.0099022438237459
## objective function = 0.0096521467408002
## objective function = 0.0094300346431928
## objective function = 0.0093082557354885
## objective function = 0.0092221192245195
## objective function = 0.0091853894974883
## objective function = 0.0091683083600652
## objective function = 0.0091613016088296
## objective function = 0.0091578781532012
## objective function = 0.0091560347242648
## objective function = 0.0091549273880945
## objective function = 0.0091542755388411
## objective function = 0.0091538680637045
## objective function = 0.0091536333198369
## objective function = 0.0091535248273743
## objective function = 0.0091534830232298
## objective function = 0.0091534669100903
## objective function = 0.0091534608246759
## objective function = 0.0091534584698128
## objective function = 0.0091534571024759
## objective function = 0.0091534560292262
## objective function = 0.0091534553902843
## objective function = 0.0091534551407963
## objective function = 0.0091534550543652
## objective function = 0.0091534550114223
## objective function = 0.0091534549884046
## objective function = 0.0091534549787392
## objective function = 0.0091534549747107
## objective function = 0.0091534549725652
## objective function = 0.0091534549713624
## objective function = 0.0091534549713624
## convergence status (0=ok): 0
## nlminb message says: relative convergence (4)
## number of iterations: 40
## number of function evaluations [objective, gradient]: 41 41
## lavoptim ... done.
## lavimplied ... done.
## lavloglik ... done.
## computing TEST for test(s) = standard ... done.
## lavbaseline ... done.
## post check ... done.
##
## Use lsem.bootstrap() for computing standard errors!
##
## -|
## ** Parameter summary
summary(lsem_weak)
## -----------------------------------------------------------------
## Local Structural Equation Model
##
## sirt 4.2-133 (2025-09-27 12:57:51)
## lavaan 0.6-19 (2024-09-26 15:50:14 UTC)
##
## R version 4.5.1 (2025-06-13 ucrt) x86_64, mingw32 | nodename=HUANG | login=mh198
##
## Function 'sirt::lsem.estimate', type='LSEM'
##
##
## Call:
## lsem.estimate(data = as.data.frame(d_apa_cog), moderator = "PC1A_z",
## moderator.grid = grid_points, lavmodel = lsem_model, h = 1.5,
## standardized = TRUE, par_invariant = weak_constraints, estimator = "ML",
## meanstructure = TRUE)
##
## Date of Analysis: 2026-06-12 14:22:05.557419
## Time difference of 6.47238 secs
## Computation Time: 6.47238
##
## Number of observations in datasets = 1538
## Used observations in analysis = 1524
## Used sampling weights: FALSE
## Bandwidth factor = 1.5
## Bandwidth = 0.294
## Number of focal points for moderator = 5
##
## Used joint estimation: TRUE
## Used sufficient statistics: TRUE
## Used local linear smoothing: TRUE
## Used pseudo weights: FALSE
## Used lavaan package: TRUE
## Used lavaan.survey package: FALSE
##
## Mean structure modelled: TRUE
##
## lavaan Model
## g =~ PPVT_z + WJ9_z + WJ10_z + DS_z
##
## Global Fit Statistics for Joint Estimation
##
## stat value
## 1 rmsea 0.052
## 2 cfi 0.990
## 3 tli 0.987
## 4 gfi 0.991
## 5 srmr 0.024
##
##
## Parameter Estimate Summary
##
## par parindex M SD MAD Min Max lin_int lin_slo
## 1 g=~PPVT_z 1 1.000 0.000 0.000 1.000 1.000 1.000 0.000
## 2 g=~WJ9_z 2 1.311 0.000 0.000 1.311 1.311 1.311 0.000
## 3 g=~WJ10_z 3 1.255 0.000 0.000 1.255 1.255 1.255 0.000
## 4 g=~DS_z 4 0.874 0.000 0.000 0.874 0.874 0.874 0.000
## 5 PPVT_z~~PPVT_z 5 0.569 0.079 0.067 0.487 0.723 0.663 -0.296
## 6 WJ9_z~~WJ9_z 6 0.322 0.031 0.027 0.283 0.382 0.359 -0.118
## 7 WJ10_z~~WJ10_z 7 0.402 0.014 0.012 0.370 0.413 0.397 0.017
## 8 DS_z~~DS_z 8 0.718 0.025 0.021 0.669 0.742 0.690 0.088
## 9 g~~g 9 0.396 0.003 0.003 0.388 0.399 0.394 0.006
## 10 PPVT_z~1 10 -0.001 0.004 0.003 -0.008 0.004 -0.001 -0.001
## 11 WJ9_z~1 11 0.001 0.003 0.002 -0.002 0.008 0.004 -0.007
## 12 WJ10_z~1 12 -0.002 0.005 0.004 -0.007 0.004 0.001 -0.008
## 13 DS_z~1 13 -0.001 0.002 0.002 -0.004 0.001 -0.001 -0.001
## 14 g~1 14 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## 15 std__g=~PPVT_z 15 0.642 0.026 0.022 0.591 0.669 0.611 0.098
## 16 std__g=~WJ9_z 16 0.824 0.013 0.011 0.797 0.840 0.808 0.050
## 17 std__g=~WJ10_z 17 0.780 0.004 0.003 0.777 0.789 0.781 -0.004
## 18 std__g=~DS_z 18 0.545 0.006 0.005 0.539 0.554 0.551 -0.021
## 19 std__PPVT_z~~PPVT_z 19 0.587 0.033 0.029 0.552 0.651 0.626 -0.124
## 20 std__WJ9_z~~WJ9_z 20 0.321 0.022 0.019 0.295 0.365 0.346 -0.082
## 21 std__WJ10_z~~WJ10_z 21 0.392 0.006 0.005 0.377 0.397 0.390 0.007
## 22 std__DS_z~~DS_z 22 0.703 0.006 0.006 0.693 0.709 0.696 0.023
## 23 std__g~~g 23 1.000 0.000 0.000 1.000 1.000 1.000 0.000
## 24 std__PPVT_z~1 24 -0.001 0.003 0.003 -0.008 0.003 -0.001 -0.001
## 25 std__WJ9_z~1 25 0.001 0.003 0.002 -0.002 0.008 0.004 -0.007
## 26 std__WJ10_z~1 26 -0.002 0.005 0.004 -0.007 0.004 0.001 -0.008
## 27 std__DS_z~1 27 -0.001 0.002 0.002 -0.004 0.001 -0.001 0.000
## 28 std__g~1 28 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## SD_nonlin
## 1 0.000
## 2 0.000
## 3 0.000
## 4 0.000
## 5 0.012
## 6 0.004
## 7 0.013
## 8 0.009
## 9 0.003
## 10 0.004
## 11 0.003
## 12 0.004
## 13 0.002
## 14 0.000
## 15 0.004
## 16 0.001
## 17 0.004
## 18 0.002
## 19 0.005
## 20 0.002
## 21 0.006
## 22 0.002
## 23 0.000
## 24 0.003
## 25 0.003
## 26 0.004
## 27 0.002
## 28 0.000
##
## Distribution of Moderator: Density and Effective Sample Size
##
## M=0.062 | SD=0.851
##
## moderator wgt Neff
## 1 -0.246 0.100 419.601
## 2 0.070 0.167 679.165
## 3 0.264 0.233 789.229
## 4 0.453 0.259 800.087
## 5 0.619 0.241 713.065
##
## variable M SD min max
## 1 moderator 0.062 0.851 -4.930 1.448
## 2 wgt 0.200 0.066 0.100 0.259
## 3 Neff 680.229 154.315 419.601 800.087
b_weak <- lsem.bootstrap(lsem_weak, R = 50, n.core = 14)
## Compute LSEM for bootstrap samples:
## Warning in sqrt(w): NaNs produced
# ── Scalar (strong) ──────────────────────────────────────────────────────
lsem_strong <- lsem.estimate(
data = as.data.frame(d_apa_cog),
moderator = "PC1A_z",
moderator.grid = grid_points,
lavmodel = lsem_model,
h = 1.5,
estimator = "ML",
meanstructure = TRUE,
standardized = TRUE,
par_invariant = strong_constraints
)
## Use joint estimation
## ** Residualize Data
## ** Fit lavaan model
## |*****|
## |lavoptions ... done.
## lavdata ... done.
## Number of observations per group:
## Group 1 415
## Group 2 673
## Group 3 782
## Group 4 793
## Group 5 707
## lavsamplestats ... done.
## lavh1 ... start:
## lavh1 ... done.
## lavpartable bounds ... done.
## lavstart ... done.
## lavmodel ... done.
## lavoptim ... start:
## attempt 1 -- default options
## objective function = 0.0725269442612355
## quasi-Newton steps using NLMINB:
## objective function = 0.0725269442612355
## objective function = 0.0566749374170481
## objective function = 0.0273724740446447
## objective function = 0.0215787767042751
## objective function = 0.0154079500627411
## objective function = 0.0142616518253383
## objective function = 0.0133152586899534
## objective function = 0.0125822646837887
## objective function = 0.0118615081359439
## objective function = 0.0106721034105831
## objective function = 0.0098718681840125
## objective function = 0.0095931567616584
## objective function = 0.0094202046031605
## objective function = 0.0093063926478368
## objective function = 0.0092455561298926
## objective function = 0.0092140212556277
## objective function = 0.0091950364054556
## objective function = 0.0091833493829232
## objective function = 0.0091776924920001
## objective function = 0.0091752531281609
## objective function = 0.0091741669187282
## objective function = 0.0091736701498349
## objective function = 0.0091734631705778
## objective function = 0.0091733807707360
## objective function = 0.0091733423329667
## objective function = 0.0091733218520558
## objective function = 0.0091733119221287
## objective function = 0.0091733084532284
## objective function = 0.0091733076292169
## objective function = 0.0091733074335419
## objective function = 0.0091733073609105
## objective function = 0.0091733073334915
## objective function = 0.0091733073251845
## objective function = 0.0091733073219766
## objective function = 0.0091733073203694
## objective function = 0.0091733073197574
## objective function = 0.0091733073197574
## convergence status (0=ok): 0
## nlminb message says: relative convergence (4)
## number of iterations: 34
## number of function evaluations [objective, gradient]: 36 35
## lavoptim ... done.
## lavimplied ... done.
## lavloglik ... done.
## computing TEST for test(s) = standard ... done.
## lavbaseline ... done.
## post check ... done.
##
## Use lsem.bootstrap() for computing standard errors!
##
## -|
## ** Parameter summary
summary(lsem_strong)
## -----------------------------------------------------------------
## Local Structural Equation Model
##
## sirt 4.2-133 (2025-09-27 12:57:51)
## lavaan 0.6-19 (2024-09-26 15:50:14 UTC)
##
## R version 4.5.1 (2025-06-13 ucrt) x86_64, mingw32 | nodename=HUANG | login=mh198
##
## Function 'sirt::lsem.estimate', type='LSEM'
##
##
## Call:
## lsem.estimate(data = as.data.frame(d_apa_cog), moderator = "PC1A_z",
## moderator.grid = grid_points, lavmodel = lsem_model, h = 1.5,
## standardized = TRUE, par_invariant = strong_constraints,
## estimator = "ML", meanstructure = TRUE)
##
## Date of Analysis: 2026-06-12 14:36:32.916872
## Time difference of 4.034776 secs
## Computation Time: 4.034776
##
## Number of observations in datasets = 1538
## Used observations in analysis = 1524
## Used sampling weights: FALSE
## Bandwidth factor = 1.5
## Bandwidth = 0.294
## Number of focal points for moderator = 5
##
## Used joint estimation: TRUE
## Used sufficient statistics: TRUE
## Used local linear smoothing: TRUE
## Used pseudo weights: FALSE
## Used lavaan package: TRUE
## Used lavaan.survey package: FALSE
##
## Mean structure modelled: TRUE
##
## lavaan Model
## g =~ PPVT_z + WJ9_z + WJ10_z + DS_z
##
## Global Fit Statistics for Joint Estimation
##
## stat value
## 1 rmsea 0.031
## 2 cfi 0.994
## 3 tli 0.995
## 4 gfi 0.991
## 5 srmr 0.024
##
##
## Parameter Estimate Summary
##
## par parindex M SD MAD Min Max lin_int
## 1 g=~PPVT_z 1 1.000 0.000 0.000 1.000 1.000 1.000
## 2 g=~WJ9_z 2 1.311 0.000 0.000 1.311 1.311 1.311
## 3 g=~WJ10_z 3 1.255 0.000 0.000 1.255 1.255 1.255
## 4 g=~DS_z 4 0.874 0.000 0.000 0.874 0.874 0.874
## 5 PPVT_z~~PPVT_z 5 0.569 0.079 0.067 0.487 0.723 0.663
## 6 WJ9_z~~WJ9_z 6 0.322 0.031 0.027 0.283 0.382 0.359
## 7 WJ10_z~~WJ10_z 7 0.402 0.014 0.012 0.370 0.413 0.397
## 8 DS_z~~DS_z 8 0.718 0.025 0.021 0.669 0.742 0.690
## 9 g~~g 9 0.396 0.003 0.003 0.388 0.399 0.394
## 10 PPVT_z~1 10 -0.001 0.000 0.000 -0.001 -0.001 -0.001
## 11 WJ9_z~1 11 0.002 0.000 0.000 0.002 0.002 0.002
## 12 WJ10_z~1 12 -0.002 0.000 0.000 -0.002 -0.002 -0.002
## 13 DS_z~1 13 -0.001 0.000 0.000 -0.001 -0.001 -0.001
## 14 g~1 14 0.000 0.000 0.000 0.000 0.000 0.000
## 15 std__g=~PPVT_z 15 0.642 0.026 0.022 0.591 0.669 0.611
## 16 std__g=~WJ9_z 16 0.824 0.013 0.011 0.797 0.840 0.808
## 17 std__g=~WJ10_z 17 0.780 0.004 0.003 0.777 0.789 0.781
## 18 std__g=~DS_z 18 0.545 0.006 0.005 0.539 0.554 0.551
## 19 std__PPVT_z~~PPVT_z 19 0.587 0.033 0.029 0.552 0.651 0.626
## 20 std__WJ9_z~~WJ9_z 20 0.321 0.022 0.019 0.295 0.365 0.346
## 21 std__WJ10_z~~WJ10_z 21 0.392 0.006 0.005 0.377 0.397 0.390
## 22 std__DS_z~~DS_z 22 0.703 0.006 0.006 0.693 0.709 0.696
## 23 std__g~~g 23 1.000 0.000 0.000 1.000 1.000 1.000
## 24 std__PPVT_z~1 24 -0.001 0.000 0.000 -0.001 -0.001 -0.001
## 25 std__WJ9_z~1 25 0.002 0.000 0.000 0.002 0.002 0.002
## 26 std__WJ10_z~1 26 -0.002 0.000 0.000 -0.002 -0.002 -0.002
## 27 std__DS_z~1 27 -0.001 0.000 0.000 -0.001 -0.001 -0.001
## 28 std__g~1 28 0.000 0.000 0.000 0.000 0.000 0.000
## lin_slo SD_nonlin
## 1 0.000 0.000
## 2 0.000 0.000
## 3 0.000 0.000
## 4 0.000 0.000
## 5 -0.296 0.012
## 6 -0.118 0.004
## 7 0.017 0.013
## 8 0.088 0.009
## 9 0.006 0.003
## 10 0.000 0.000
## 11 0.000 0.000
## 12 0.000 0.000
## 13 0.000 0.000
## 14 0.000 0.000
## 15 0.098 0.004
## 16 0.050 0.001
## 17 -0.004 0.004
## 18 -0.021 0.002
## 19 -0.124 0.005
## 20 -0.082 0.002
## 21 0.007 0.006
## 22 0.023 0.002
## 23 0.000 0.000
## 24 0.000 0.000
## 25 0.000 0.000
## 26 0.000 0.000
## 27 0.000 0.000
## 28 0.000 0.000
##
## Distribution of Moderator: Density and Effective Sample Size
##
## M=0.062 | SD=0.851
##
## moderator wgt Neff
## 1 -0.246 0.100 419.601
## 2 0.070 0.167 679.165
## 3 0.264 0.233 789.229
## 4 0.453 0.259 800.087
## 5 0.619 0.241 713.065
##
## variable M SD min max
## 1 moderator 0.062 0.851 -4.930 1.448
## 2 wgt 0.200 0.066 0.100 0.259
## 3 Neff 680.229 154.315 419.601 800.087
b_strong <- lsem.bootstrap(lsem_strong, R = 50, n.core = 14)
## Compute LSEM for bootstrap samples:
## Warning in sqrt(w): NaNs produced
# ── Strict ───────────────────────────────────────────────────────────────
lsem_strict <- lsem.estimate(
data = as.data.frame(d_apa_cog),
moderator = "PC1A_z",
moderator.grid = grid_points,
lavmodel = lsem_model,
h = 1.5,
estimator = "ML",
meanstructure = TRUE,
standardized = TRUE,
par_invariant = strict_constraints
)
## Use joint estimation
## ** Residualize Data
## ** Fit lavaan model
## |*****|
## |lavoptions ... done.
## lavdata ... done.
## Number of observations per group:
## Group 1 415
## Group 2 673
## Group 3 782
## Group 4 793
## Group 5 707
## lavsamplestats ... done.
## lavh1 ... start:
## lavh1 ... done.
## lavpartable bounds ... done.
## lavstart ... done.
## lavmodel ... done.
## lavoptim ... start:
## attempt 1 -- default options
## objective function = 0.0725269442612357
## quasi-Newton steps using NLMINB:
## objective function = 0.0725269442612357
## objective function = 0.0598849149521864
## objective function = 0.0456558227751016
## objective function = 0.0280676116603524
## objective function = 0.0233835766382247
## objective function = 0.0203046704361951
## objective function = 0.0172907948008855
## objective function = 0.0166145122110242
## objective function = 0.0152118914595368
## objective function = 0.0149220340810096
## objective function = 0.0146726699976791
## objective function = 0.0145038426628579
## objective function = 0.0143994581838247
## objective function = 0.0143413746876366
## objective function = 0.0143237628606473
## objective function = 0.0143173566351997
## objective function = 0.0143133119428696
## objective function = 0.0143117244210845
## objective function = 0.0143112455442859
## objective function = 0.0143109822916277
## objective function = 0.0143107648550805
## objective function = 0.0143106236481005
## objective function = 0.0143105939605634
## objective function = 0.0143105894862483
## objective function = 0.0143105889594805
## objective function = 0.0143105888699047
## objective function = 0.0143105888529152
## objective function = 0.0143105888463410
## objective function = 0.0143105888416870
## objective function = 0.0143105888387537
## objective function = 0.0143105888370059
## objective function = 0.0143105888370059
## convergence status (0=ok): 0
## nlminb message says: relative convergence (4)
## number of iterations: 30
## number of function evaluations [objective, gradient]: 31 31
## lavoptim ... done.
## lavimplied ... done.
## lavloglik ... done.
## computing TEST for test(s) = standard ... done.
## lavbaseline ... done.
## post check ... done.
##
## Use lsem.bootstrap() for computing standard errors!
##
## -|
## ** Parameter summary
summary(lsem_strict)
## -----------------------------------------------------------------
## Local Structural Equation Model
##
## sirt 4.2-133 (2025-09-27 12:57:51)
## lavaan 0.6-19 (2024-09-26 15:50:14 UTC)
##
## R version 4.5.1 (2025-06-13 ucrt) x86_64, mingw32 | nodename=HUANG | login=mh198
##
## Function 'sirt::lsem.estimate', type='LSEM'
##
##
## Call:
## lsem.estimate(data = as.data.frame(d_apa_cog), moderator = "PC1A_z",
## moderator.grid = grid_points, lavmodel = lsem_model, h = 1.5,
## standardized = TRUE, par_invariant = strict_constraints,
## estimator = "ML", meanstructure = TRUE)
##
## Date of Analysis: 2026-06-12 14:50:23.673098
## Time difference of 8.916642 secs
## Computation Time: 8.916642
##
## Number of observations in datasets = 1538
## Used observations in analysis = 1524
## Used sampling weights: FALSE
## Bandwidth factor = 1.5
## Bandwidth = 0.294
## Number of focal points for moderator = 5
##
## Used joint estimation: TRUE
## Used sufficient statistics: TRUE
## Used local linear smoothing: TRUE
## Used pseudo weights: FALSE
## Used lavaan package: TRUE
## Used lavaan.survey package: FALSE
##
## Mean structure modelled: TRUE
##
## lavaan Model
## g =~ PPVT_z + WJ9_z + WJ10_z + DS_z
##
## Global Fit Statistics for Joint Estimation
##
## stat value
## 1 rmsea 0.034
## 2 cfi 0.990
## 3 tli 0.994
## 4 gfi 0.986
## 5 srmr 0.030
##
##
## Parameter Estimate Summary
##
## par parindex M SD MAD Min Max lin_int
## 1 g=~PPVT_z 1 1.000 0.000 0.000 1.000 1.000 1.000
## 2 g=~WJ9_z 2 1.322 0.000 0.000 1.322 1.322 1.322
## 3 g=~WJ10_z 3 1.260 0.000 0.000 1.260 1.260 1.260
## 4 g=~DS_z 4 0.879 0.000 0.000 0.879 0.879 0.879
## 5 PPVT_z~~PPVT_z 5 0.581 0.000 0.000 0.581 0.581 0.581
## 6 WJ9_z~~WJ9_z 6 0.324 0.000 0.000 0.324 0.324 0.324
## 7 WJ10_z~~WJ10_z 7 0.403 0.000 0.000 0.403 0.403 0.403
## 8 DS_z~~DS_z 8 0.715 0.000 0.000 0.715 0.715 0.715
## 9 g~~g 9 0.391 0.005 0.004 0.383 0.396 0.394
## 10 PPVT_z~1 10 -0.001 0.000 0.000 -0.001 -0.001 -0.001
## 11 WJ9_z~1 11 0.002 0.000 0.000 0.002 0.002 0.002
## 12 WJ10_z~1 12 -0.002 0.000 0.000 -0.002 -0.002 -0.002
## 13 DS_z~1 13 -0.001 0.000 0.000 -0.001 -0.001 -0.001
## 14 g~1 14 0.000 0.000 0.000 0.000 0.000 0.000
## 15 std__g=~PPVT_z 15 0.634 0.002 0.002 0.630 0.637 0.636
## 16 std__g=~WJ9_z 16 0.824 0.002 0.001 0.821 0.825 0.825
## 17 std__g=~WJ10_z 17 0.779 0.002 0.002 0.776 0.781 0.780
## 18 std__g=~DS_z 18 0.545 0.003 0.002 0.541 0.548 0.547
## 19 std__PPVT_z~~PPVT_z 19 0.598 0.003 0.003 0.594 0.603 0.596
## 20 std__WJ9_z~~WJ9_z 20 0.322 0.003 0.002 0.319 0.326 0.320
## 21 std__WJ10_z~~WJ10_z 21 0.394 0.003 0.003 0.390 0.398 0.392
## 22 std__DS_z~~DS_z 22 0.703 0.003 0.002 0.700 0.707 0.701
## 23 std__g~~g 23 1.000 0.000 0.000 1.000 1.000 1.000
## 24 std__PPVT_z~1 24 -0.001 0.000 0.000 -0.001 -0.001 -0.001
## 25 std__WJ9_z~1 25 0.002 0.000 0.000 0.002 0.002 0.002
## 26 std__WJ10_z~1 26 -0.002 0.000 0.000 -0.002 -0.002 -0.002
## 27 std__DS_z~1 27 -0.001 0.000 0.000 -0.001 -0.001 -0.001
## 28 std__g~1 28 0.000 0.000 0.000 0.000 0.000 0.000
## lin_slo SD_nonlin
## 1 0.000 0.000
## 2 0.000 0.000
## 3 0.000 0.000
## 4 0.000 0.000
## 5 0.000 0.000
## 6 0.000 0.000
## 7 0.000 0.000
## 8 0.000 0.000
## 9 -0.010 0.004
## 10 0.000 0.000
## 11 0.000 0.000
## 12 0.000 0.000
## 13 0.000 0.000
## 14 0.000 0.000
## 15 -0.005 0.002
## 16 -0.003 0.002
## 17 -0.004 0.002
## 18 -0.005 0.002
## 19 0.006 0.003
## 20 0.006 0.002
## 21 0.006 0.003
## 22 0.005 0.002
## 23 0.000 0.000
## 24 0.000 0.000
## 25 0.000 0.000
## 26 0.000 0.000
## 27 0.000 0.000
## 28 0.000 0.000
##
## Distribution of Moderator: Density and Effective Sample Size
##
## M=0.062 | SD=0.851
##
## moderator wgt Neff
## 1 -0.246 0.100 419.601
## 2 0.070 0.167 679.165
## 3 0.264 0.233 789.229
## 4 0.453 0.259 800.087
## 5 0.619 0.241 713.065
##
## variable M SD min max
## 1 moderator 0.062 0.851 -4.930 1.448
## 2 wgt 0.200 0.066 0.100 0.259
## 3 Neff 680.229 154.315 419.601 800.087
b_strict <- lsem.bootstrap(lsem_strict, R = 50, n.core = 14)
## Compute LSEM for bootstrap samples:
## Warning in sqrt(w): NaNs produced
# ── Compare fit across invariance levels ─────────────────────────────────
invar_test <- rbind(b_config$fitstats_joint,
b_weak$fitstats_joint,
b_strong$fitstats_joint,
b_strict$fitstats_joint) %>%
mutate(ci_l = value - 1.96 * se,
ci_u = value + 1.96 * se)
invar_test$level <- c(rep("configural", nrow(b_config$fitstats_joint)),
rep("metric", nrow(b_weak$fitstats_joint)),
rep("scalar", nrow(b_strong$fitstats_joint)),
rep("strict", nrow(b_strict$fitstats_joint)))
invar_test$level <- factor(invar_test$level, levels = unique(invar_test$level))
print(invar_test)
## stat value value_bc se ci_l ci_u
## rmsea rmsea 0.08187480 0.07453702 0.022479901 0.037814194 0.12593540
## cfi cfi 0.98886075 0.99148966 0.006223548 0.976662595 1.00105890
## tli tli 0.96658225 0.97446897 0.018670643 0.929987786 1.00317671
## gfi gfi 0.99207405 0.99380322 0.003624108 0.984970803 0.99917731
## srmr srmr 0.01987493 0.01828622 0.004311263 0.011424853 0.02832501
## rmsea1 rmsea 0.05173956 0.04124471 0.015231831 0.021885170 0.08159395
## cfi1 cfi 0.99021359 0.99523860 0.006643360 0.977192600 1.00323457
## tli1 tli 0.98665489 0.99350718 0.009059128 0.968898999 1.00441078
## gfi1 gfi 0.99119013 0.99421749 0.003888456 0.983568756 0.99881151
## srmr1 srmr 0.02425864 0.01626246 0.005372913 0.013727736 0.03478955
## rmsea2 rmsea 0.03050162 0.01947504 0.016554626 -0.001945449 0.06294869
## cfi2 cfi 0.99412531 1.00061827 0.008288920 0.977879029 1.01037159
## tli2 tli 0.99536209 1.00045923 0.006589323 0.982447014 1.00827716
## gfi2 gfi 0.99117039 0.99492553 0.004651413 0.982053621 1.00028716
## srmr2 srmr 0.02432868 0.01727356 0.005738153 0.013081902 0.03557546
## rmsea3 rmsea 0.03415303 0.02385697 0.011172281 0.012255361 0.05605070
## cfi3 cfi 0.98953335 0.99790062 0.009031757 0.971831106 1.00723559
## tli3 tli 0.99418519 0.99883368 0.005017643 0.984350614 1.00401977
## gfi3 gfi 0.98630652 0.99103597 0.004880105 0.976741515 0.99587153
## srmr3 srmr 0.03023610 0.02046019 0.006595546 0.017308827 0.04316337
## level
## rmsea configural
## cfi configural
## tli configural
## gfi configural
## srmr configural
## rmsea1 metric
## cfi1 metric
## tli1 metric
## gfi1 metric
## srmr1 metric
## rmsea2 scalar
## cfi2 scalar
## tli2 scalar
## gfi2 scalar
## srmr2 scalar
## rmsea3 strict
## cfi3 strict
## tli3 strict
## gfi3 strict
## srmr3 strict
# ── Plots ────────────────────────────────────────────────────────────────
p_cfi <- invar_test %>%
filter(stat == "cfi") %>%
ggplot(aes(x = level, y = value)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = ci_l, ymax = ci_u), width = .2) +
scale_y_continuous(limits = c(.9, 1)) +
ylab("CFI") + xlab("Invariance Level") +
theme_minimal()
p_rmsea <- invar_test %>%
filter(stat == "rmsea") %>%
ggplot(aes(x = level, y = value)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = ci_l, ymax = ci_u), width = .2) +
scale_y_continuous(limits = c(0, .15)) +
ylab("RMSEA") + xlab("Invariance Level") +
theme_minimal()
p_srmr <- invar_test %>%
filter(stat == "srmr") %>%
ggplot(aes(x = level, y = value)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = ci_l, ymax = ci_u), width = .2) +
scale_y_continuous(limits = c(0, .15)) +
ylab("SRMR") + xlab("Invariance Level") +
theme_minimal()
print(p_cfi)
print(p_rmsea)
print(p_srmr)
#LevelRMSEACFITLISRMRConfigural.082 (.031–.131).989 (.973–1.005).967 (.919–1.015).020 (.010–.030)
#Metric.051 (.021–.082).990 (.976–1.005).987 (.967–1.006).024 (.013–.035)
#Scalar.030 (.001–.059).994 (.982–1.006).995 (.986–1.005).024 (.014–.034)
#Strict.034 (.010–.058).990 (.973–1.006).994 (.985–1.004).030 (.020–.041)
# ==============================================================
# 11b. Multi-Group CFA: Measurement Invariance (PC1A tertiles)
# ==============================================================
# Direct test of configural -> metric -> scalar -> strict.
# Complements LSEM: LSEM tests continuous parameter moderation
# but doesn't explicitly test intercepts. This confirms scalar
# invariance directly.
cfa_model <- 'g =~ PPVT_z + WJ9_z + WJ10_z + DS_z'
run_mi_cfa <- function(data, pc_var, label) {
dd <- data[complete.cases(data[, c("PPVT_z", "WJ9_z", "WJ10_z", "DS_z", pc_var)]), ]
dd$tert <- cut(dd[[pc_var]],
breaks = quantile(dd[[pc_var]], probs = c(0, 1/3, 2/3, 1)),
include.lowest = TRUE, labels = c("Low", "Mid", "High"))
cat("\n", paste(rep("=", 80), collapse = ""), "\n")
cat(label, "\n")
cat("N =", nrow(dd), " | Tertile Ns:", table(dd$tert), "\n")
cat(paste(rep("=", 80), collapse = ""), "\n")
fit_config <- cfa(cfa_model, data = dd, group = "tert", estimator = "ML",
meanstructure = TRUE)
fit_metric <- cfa(cfa_model, data = dd, group = "tert", estimator = "ML",
meanstructure = TRUE, group.equal = "loadings")
fit_scalar <- cfa(cfa_model, data = dd, group = "tert", estimator = "ML",
meanstructure = TRUE, group.equal = c("loadings", "intercepts"))
fit_strict <- cfa(cfa_model, data = dd, group = "tert", estimator = "ML",
meanstructure = TRUE, group.equal = c("loadings", "intercepts", "residuals"))
extract_fit <- function(fit, lbl) {
fm <- fitmeasures(fit, c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr"))
data.frame(Model = lbl,
chisq = round(fm["chisq"], 3),
df = fm["df"],
p = round(fm["pvalue"], 4),
CFI = round(fm["cfi"], 4),
TLI = round(fm["tli"], 4),
RMSEA = round(fm["rmsea"], 4),
SRMR = round(fm["srmr"], 4),
row.names = NULL)
}
fit_table <- rbind(
extract_fit(fit_config, "1. Configural"),
extract_fit(fit_metric, "2. Metric"),
extract_fit(fit_scalar, "3. Scalar"),
extract_fit(fit_strict, "4. Strict")
)
print(fit_table, row.names = FALSE)
cfis <- sapply(list(fit_config, fit_metric, fit_scalar, fit_strict), fitmeasures, "cfi")
steps <- c("Config->Metric", "Metric->Scalar", "Scalar->Strict")
cat("\n")
for (i in 1:3) {
d_cfi <- cfis[i + 1] - cfis[i]
cat(sprintf(" %-20s ΔCFI = %+.4f %s\n", steps[i], d_cfi,
ifelse(abs(d_cfi) < 0.01, "[PASS]", "[FAIL]")))
}
cat("\n")
}
run_mi_cfa(d_apa, "K5PC1A", "APA x African Ancestry PC (PC1A) tertiles")
##
## ================================================================================
## APA x African Ancestry PC (PC1A) tertiles
## N = 1507 | Tertile Ns: 503 502 502
## ================================================================================
## Model chisq df p CFI TLI RMSEA SRMR
## 1. Configural 27.049 6 0.0001 0.9882 0.9647 0.0836 0.0202
## 2. Metric 30.072 12 0.0027 0.9899 0.9849 0.0548 0.0245
## 3. Scalar 44.348 18 0.0005 0.9853 0.9853 0.0540 0.0316
## 4. Strict 64.096 26 0.0000 0.9787 0.9853 0.0540 0.0424
##
## Config->Metric ΔCFI = +0.0017 [PASS]
## Metric->Scalar ΔCFI = -0.0046 [PASS]
## Scalar->Strict ΔCFI = -0.0066 [PASS]
run_mi_cfa(d_haa, "K5PC1A", "HAA x African Ancestry PC (PC1A) tertiles")
##
## ================================================================================
## HAA x African Ancestry PC (PC1A) tertiles
## N = 1261 | Tertile Ns: 421 420 420
## ================================================================================
## Model chisq df p CFI TLI RMSEA SRMR
## 1. Configural 24.688 6 0.0004 0.9877 0.9630 0.0861 0.0209
## 2. Metric 27.769 12 0.0060 0.9896 0.9844 0.0559 0.0258
## 3. Scalar 35.893 18 0.0073 0.9882 0.9882 0.0486 0.0302
## 4. Strict 66.951 26 0.0000 0.9730 0.9813 0.0612 0.0439
##
## Config->Metric ΔCFI = +0.0019 [PASS]
## Metric->Scalar ΔCFI = -0.0014 [PASS]
## Scalar->Strict ΔCFI = -0.0152 [FAIL]
# ==============================================================
# 11c. LSEM Bandwidth Sensitivity
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("ROBUSTNESS: LSEM bandwidth sensitivity (configural only)\n")
## ROBUSTNESS: LSEM bandwidth sensitivity (configural only)
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
for (h_val in c(1.0, 2.0)) {
lsem_h <- lsem.estimate(
data = as.data.frame(d_apa_cog),
moderator = "PC1A_z",
moderator.grid = grid_points,
lavmodel = lsem_model,
h = h_val,
estimator = "ML",
meanstructure = TRUE,
standardized = TRUE,
est_joint = TRUE
)
cat("\nh =", h_val, "\n")
print(summary(lsem_h)$fitstat_joint)
}
## ** Residualize Data
## ** Fit lavaan model
## |*****|
## |lavoptions ... done.
## lavdata ... done.
## Number of observations per group:
## Group 1 260
## Group 2 469
## Group 3 568
## Group 4 592
## Group 5 526
## lavsamplestats ... done.
## lavh1 ... start:
## lavh1 ... done.
## lavpartable bounds ... done.
## lavstart ... done.
## lavmodel ... done.
## lavoptim ... start:
## attempt 1 -- default options
## objective function = 0.4754973947033598
## quasi-Newton steps using NLMINB:
## objective function = 0.4754973947033598
## objective function = 0.4310407909366635
## objective function = 0.3143510820102698
## objective function = 0.1542020094541372
## objective function = 1.1001605993795394
## objective function = 0.1304131229638222
## objective function = 0.0935235360325728
## objective function = 0.0744449729434417
## objective function = 0.1136588580712203
## objective function = 0.0577928740295930
## objective function = 0.0501947218570290
## objective function = 0.0362155729052460
## objective function = 0.0274337453316155
## objective function = 0.0228531030155932
## objective function = 0.0196114566771164
## objective function = 0.0175016767418163
## objective function = 0.0157431992278658
## objective function = 0.0140802066247961
## objective function = 0.0130883347630556
## objective function = 0.0122756731364266
## objective function = 0.0116828929902572
## objective function = 0.0111388413994315
## objective function = 0.0106484804178050
## objective function = 0.0103234912375580
## objective function = 0.0101092501167388
## objective function = 0.0099305264135944
## objective function = 0.0096936322734374
## objective function = 0.0096126045050253
## objective function = 0.0095237734414318
## objective function = 0.0094500176843905
## objective function = 0.0093797700575670
## objective function = 0.0093163318601537
## objective function = 0.0092630988053402
## objective function = 0.0092249482170084
## objective function = 0.0092017357033331
## objective function = 0.0091843593734885
## objective function = 0.0091725302905539
## objective function = 0.0091635319056707
## objective function = 0.0091576760816100
## objective function = 0.0091538825444771
## objective function = 0.0091511195656637
## objective function = 0.0091487887009347
## objective function = 0.0091462160392724
## objective function = 0.0091440639201876
## objective function = 0.0091419113645792
## objective function = 0.0091400849735483
## objective function = 0.0091386297166956
## objective function = 0.0091374061436663
## objective function = 0.0091365038108760
## objective function = 0.0091358786871588
## objective function = 0.0091354631727697
## objective function = 0.0091351619084838
## objective function = 0.0091348781332385
## objective function = 0.0091346136740585
## objective function = 0.0091344016601473
## objective function = 0.0091342626298315
## objective function = 0.0091341741928641
## objective function = 0.0091341078379438
## objective function = 0.0091340492025364
## objective function = 0.0091340040948203
## objective function = 0.0091339679742639
## objective function = 0.0091339419812930
## objective function = 0.0091339244067327
## objective function = 0.0091339121810573
## objective function = 0.0091339050142846
## objective function = 0.0091339000311443
## objective function = 0.0091338964401904
## objective function = 0.0091338939518031
## objective function = 0.0091338924274136
## objective function = 0.0091338915867166
## objective function = 0.0091338912090904
## objective function = 0.0091338910702548
## objective function = 0.0091338910107239
## objective function = 0.0091338909707185
## objective function = 0.0091338909304676
## objective function = 0.0091338909029609
## objective function = 0.0091338908791812
## objective function = 0.0091338908622966
## objective function = 0.0091338908502639
## objective function = 0.0091338908443868
## objective function = 0.0091338908419006
## objective function = 0.0091338908404552
## objective function = 0.0091338908393481
## objective function = 0.0091338908381443
## objective function = 0.0091338908376613
## objective function = 0.0091338908376613
## convergence status (0=ok): 0
## nlminb message says: relative convergence (4)
## number of iterations: 80
## number of function evaluations [objective, gradient]: 85 81
## lavoptim ... done.
## lavimplied ... done.
## lavloglik ... done.
## computing TEST for test(s) = standard ... done.
## lavbaseline ... done.
## post check ... done.
##
## Use lsem.bootstrap() for computing standard errors!
##
## -|
## ** Parameter summary
##
##
## h = 1
## -----------------------------------------------------------------
## Local Structural Equation Model
##
## sirt 4.2-133 (2025-09-27 12:57:51)
## lavaan 0.6-19 (2024-09-26 15:50:14 UTC)
##
## R version 4.5.1 (2025-06-13 ucrt) x86_64, mingw32 | nodename=HUANG | login=mh198
##
## Function 'sirt::lsem.estimate', type='LSEM'
##
##
## Call:
## lsem.estimate(data = as.data.frame(d_apa_cog), moderator = "PC1A_z",
## moderator.grid = grid_points, lavmodel = lsem_model, h = h_val,
## standardized = TRUE, est_joint = TRUE, estimator = "ML",
## meanstructure = TRUE)
##
## Date of Analysis: 2026-06-12 15:04:32.017989
## Time difference of 4.17801 secs
## Computation Time: 4.17801
##
## Number of observations in datasets = 1538
## Used observations in analysis = 1524
## Used sampling weights: FALSE
## Bandwidth factor = 1
## Bandwidth = 0.196
## Number of focal points for moderator = 5
##
## Used joint estimation: TRUE
## Used sufficient statistics: TRUE
## Used local linear smoothing: TRUE
## Used pseudo weights: FALSE
## Used lavaan package: TRUE
## Used lavaan.survey package: FALSE
##
## Mean structure modelled: TRUE
##
## lavaan Model
## g =~ PPVT_z + WJ9_z + WJ10_z + DS_z
##
## Global Fit Statistics for Joint Estimation
##
## stat value
## 1 rmsea 0.084
## 2 cfi 0.988
## 3 tli 0.965
## 4 gfi 0.991
## 5 srmr 0.021
##
##
## Parameter Estimate Summary
##
## par parindex M SD MAD Min Max lin_int lin_slo
## 1 g=~PPVT_z 1 1.000 0.000 0.000 1.000 1.000 1.000 0.000
## 2 g=~WJ9_z 2 1.317 0.099 0.090 1.214 1.485 1.419 -0.325
## 3 g=~WJ10_z 3 1.256 0.060 0.048 1.209 1.407 1.324 -0.214
## 4 g=~DS_z 4 0.879 0.037 0.027 0.848 0.975 0.911 -0.102
## 5 PPVT_z~~PPVT_z 5 0.563 0.096 0.082 0.469 0.733 0.674 -0.352
## 6 WJ9_z~~WJ9_z 6 0.316 0.014 0.011 0.294 0.341 0.330 -0.047
## 7 WJ10_z~~WJ10_z 7 0.409 0.018 0.016 0.373 0.422 0.407 0.005
## 8 DS_z~~DS_z 8 0.723 0.034 0.031 0.662 0.760 0.685 0.120
## 9 g~~g 9 0.397 0.038 0.032 0.327 0.431 0.353 0.140
## 10 PPVT_z~1 10 0.001 0.006 0.004 -0.008 0.013 0.002 -0.001
## 11 WJ9_z~1 11 0.002 0.005 0.004 -0.002 0.013 0.003 -0.004
## 12 WJ10_z~1 12 0.000 0.004 0.004 -0.005 0.005 0.001 -0.005
## 13 DS_z~1 13 -0.001 0.002 0.002 -0.003 0.003 0.000 -0.002
## 14 g~1 14 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## 15 std__g=~PPVT_z 15 0.644 0.050 0.043 0.555 0.692 0.586 0.183
## 16 std__g=~WJ9_z 16 0.827 0.008 0.007 0.809 0.838 0.825 0.004
## 17 std__g=~WJ10_z 17 0.777 0.011 0.009 0.764 0.797 0.776 0.004
## 18 std__g=~DS_z 18 0.545 0.009 0.008 0.534 0.565 0.547 -0.005
## 19 std__PPVT_z~~PPVT_z 19 0.583 0.063 0.054 0.521 0.692 0.655 -0.230
## 20 std__WJ9_z~~WJ9_z 20 0.317 0.014 0.011 0.298 0.346 0.318 -0.006
## 21 std__WJ10_z~~WJ10_z 21 0.396 0.016 0.014 0.366 0.416 0.398 -0.006
## 22 std__DS_z~~DS_z 22 0.703 0.010 0.009 0.680 0.715 0.701 0.005
## 23 std__g~~g 23 1.000 0.000 0.000 1.000 1.000 1.000 0.000
## 24 std__PPVT_z~1 24 0.001 0.006 0.004 -0.008 0.012 0.002 -0.001
## 25 std__WJ9_z~1 25 0.002 0.005 0.004 -0.002 0.012 0.003 -0.003
## 26 std__WJ10_z~1 26 0.000 0.004 0.003 -0.005 0.005 0.001 -0.005
## 27 std__DS_z~1 27 -0.001 0.002 0.002 -0.003 0.003 0.000 -0.002
## 28 std__g~1 28 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## SD_nonlin
## 1 0.000
## 2 0.051
## 3 0.022
## 4 0.025
## 5 0.024
## 6 0.006
## 7 0.018
## 8 0.014
## 9 0.011
## 10 0.006
## 11 0.005
## 12 0.004
## 13 0.002
## 14 0.000
## 15 0.013
## 16 0.008
## 17 0.010
## 18 0.009
## 19 0.016
## 20 0.014
## 21 0.016
## 22 0.010
## 23 0.000
## 24 0.005
## 25 0.005
## 26 0.004
## 27 0.002
## 28 0.000
##
## Distribution of Moderator: Density and Effective Sample Size
##
## M=0.062 | SD=0.851
##
## moderator wgt Neff
## 1 -0.246 0.100 262.694
## 2 0.070 0.167 473.933
## 3 0.264 0.233 572.817
## 4 0.453 0.259 597.446
## 5 0.619 0.241 531.044
##
## variable M SD min max
## 1 moderator 0.062 0.851 -4.930 1.448
## 2 wgt 0.200 0.066 0.100 0.259
## 3 Neff 487.587 134.150 262.694 597.446
## NULL
## ** Residualize Data
## ** Fit lavaan model
## |*****|
## |lavoptions ... done.
## lavdata ... done.
## Number of observations per group:
## Group 1 575
## Group 2 840
## Group 3 936
## Group 4 934
## Group 5 846
## lavsamplestats ... done.
## lavh1 ... start:
## lavh1 ... done.
## lavpartable bounds ... done.
## lavstart ... done.
## lavmodel ... done.
## lavoptim ... start:
## attempt 1 -- default options
## objective function = 0.0363330532126399
## quasi-Newton steps using NLMINB:
## objective function = 0.0363330532126399
## objective function = 0.0279370167711291
## objective function = 0.0153857529769107
## objective function = 0.0133272104962442
## objective function = 0.0110416922998723
## objective function = 0.0102301237819656
## objective function = 0.0095860777667646
## objective function = 0.0090282965794805
## objective function = 0.0086508329545345
## objective function = 0.0084524089896739
## objective function = 0.0082043021034368
## objective function = 0.0079977449504251
## objective function = 0.0078213564944291
## objective function = 0.0076833600586591
## objective function = 0.0076053558631472
## objective function = 0.0075428939581485
## objective function = 0.0074975157835587
## objective function = 0.0074639541003761
## objective function = 0.0074402058494830
## objective function = 0.0074234696281078
## objective function = 0.0074105597352218
## objective function = 0.0073997951802440
## objective function = 0.0073912982889315
## objective function = 0.0073844255708944
## objective function = 0.0073782324906850
## objective function = 0.0073724830020676
## objective function = 0.0073661733332528
## objective function = 0.0073617432153425
## objective function = 0.0073585362794566
## objective function = 0.0073566714062932
## objective function = 0.0073551657051661
## objective function = 0.0073538140686926
## objective function = 0.0073526451094608
## objective function = 0.0073516220349611
## objective function = 0.0073511537137943
## objective function = 0.0073508302576434
## objective function = 0.0073506373226652
## objective function = 0.0073504693542614
## objective function = 0.0073502867918172
## objective function = 0.0073501464427732
## objective function = 0.0073500134091112
## objective function = 0.0073498980684367
## objective function = 0.0073498192366753
## objective function = 0.0073497578053155
## objective function = 0.0073497099858648
## objective function = 0.0073496775789855
## objective function = 0.0073496594640100
## objective function = 0.0073496505572901
## objective function = 0.0073496452634557
## objective function = 0.0073496411753814
## objective function = 0.0073496379519627
## objective function = 0.0073496355783971
## objective function = 0.0073496338176442
## objective function = 0.0073496324457824
## objective function = 0.0073496314035622
## objective function = 0.0073496306685951
## objective function = 0.0073496301655286
## objective function = 0.0073496298094182
## objective function = 0.0073496295472912
## objective function = 0.0073496293551493
## objective function = 0.0073496292246976
## objective function = 0.0073496291466901
## objective function = 0.0073496291015490
## objective function = 0.0073496290698576
## objective function = 0.0073496290427915
## objective function = 0.0073496290200234
## objective function = 0.0073496290026585
## objective function = 0.0073496289948418
## objective function = 0.0073496289916240
## objective function = 0.0073496289905358
## objective function = 0.0073496289900420
## objective function = 0.0073496289900420
## convergence status (0=ok): 0
## nlminb message says: relative convergence (4)
## number of iterations: 70
## number of function evaluations [objective, gradient]: 71 71
## lavoptim ... done.
## lavimplied ... done.
## lavloglik ... done.
## computing TEST for test(s) = standard ... done.
## lavbaseline ... done.
## post check ... done.
##
## Use lsem.bootstrap() for computing standard errors!
##
## -|
## ** Parameter summary
##
##
## h = 2
## -----------------------------------------------------------------
## Local Structural Equation Model
##
## sirt 4.2-133 (2025-09-27 12:57:51)
## lavaan 0.6-19 (2024-09-26 15:50:14 UTC)
##
## R version 4.5.1 (2025-06-13 ucrt) x86_64, mingw32 | nodename=HUANG | login=mh198
##
## Function 'sirt::lsem.estimate', type='LSEM'
##
##
## Call:
## lsem.estimate(data = as.data.frame(d_apa_cog), moderator = "PC1A_z",
## moderator.grid = grid_points, lavmodel = lsem_model, h = h_val,
## standardized = TRUE, est_joint = TRUE, estimator = "ML",
## meanstructure = TRUE)
##
## Date of Analysis: 2026-06-12 15:04:36.554083
## Time difference of 4.390869 secs
## Computation Time: 4.390869
##
## Number of observations in datasets = 1538
## Used observations in analysis = 1524
## Used sampling weights: FALSE
## Bandwidth factor = 2
## Bandwidth = 0.392
## Number of focal points for moderator = 5
##
## Used joint estimation: TRUE
## Used sufficient statistics: TRUE
## Used local linear smoothing: TRUE
## Used pseudo weights: FALSE
## Used lavaan package: TRUE
## Used lavaan.survey package: FALSE
##
## Mean structure modelled: TRUE
##
## lavaan Model
## g =~ PPVT_z + WJ9_z + WJ10_z + DS_z
##
## Global Fit Statistics for Joint Estimation
##
## stat value
## 1 rmsea 0.078
## 2 cfi 0.990
## 3 tli 0.969
## 4 gfi 0.993
## 5 srmr 0.019
##
##
## Parameter Estimate Summary
##
## par parindex M SD MAD Min Max lin_int lin_slo
## 1 g=~PPVT_z 1 1.000 0.000 0.000 1.000 1.000 1.000 0.000
## 2 g=~WJ9_z 2 1.308 0.079 0.071 1.208 1.433 1.400 -0.291
## 3 g=~WJ10_z 3 1.255 0.054 0.042 1.200 1.380 1.318 -0.202
## 4 g=~DS_z 4 0.874 0.027 0.021 0.847 0.940 0.906 -0.102
## 5 PPVT_z~~PPVT_z 5 0.569 0.077 0.067 0.482 0.720 0.661 -0.292
## 6 WJ9_z~~WJ9_z 6 0.326 0.028 0.019 0.304 0.403 0.357 -0.097
## 7 WJ10_z~~WJ10_z 7 0.398 0.019 0.016 0.356 0.416 0.395 0.009
## 8 DS_z~~DS_z 8 0.715 0.017 0.014 0.674 0.730 0.695 0.061
## 9 g~~g 9 0.398 0.034 0.029 0.330 0.437 0.357 0.129
## 10 PPVT_z~1 10 -0.002 0.002 0.002 -0.004 0.000 -0.001 -0.003
## 11 WJ9_z~1 11 0.000 0.003 0.003 -0.003 0.005 0.004 -0.011
## 12 WJ10_z~1 12 -0.002 0.004 0.004 -0.007 0.004 0.002 -0.012
## 13 DS_z~1 13 -0.002 0.001 0.001 -0.005 0.000 -0.002 0.001
## 14 g~1 14 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## 15 std__g=~PPVT_z 15 0.642 0.042 0.036 0.560 0.690 0.592 0.157
## 16 std__g=~WJ9_z 16 0.822 0.010 0.006 0.792 0.827 0.814 0.024
## 17 std__g=~WJ10_z 17 0.781 0.008 0.007 0.773 0.799 0.782 -0.001
## 18 std__g=~DS_z 18 0.545 0.003 0.003 0.541 0.550 0.544 0.004
## 19 std__PPVT_z~~PPVT_z 19 0.586 0.052 0.046 0.525 0.686 0.649 -0.198
## 20 std__WJ9_z~~WJ9_z 20 0.325 0.016 0.010 0.315 0.373 0.337 -0.039
## 21 std__WJ10_z~~WJ10_z 21 0.389 0.013 0.011 0.362 0.403 0.389 0.002
## 22 std__DS_z~~DS_z 22 0.703 0.004 0.003 0.698 0.707 0.704 -0.004
## 23 std__g~~g 23 1.000 0.000 0.000 1.000 1.000 1.000 0.000
## 24 std__PPVT_z~1 24 -0.002 0.002 0.002 -0.005 0.000 -0.001 -0.004
## 25 std__WJ9_z~1 25 0.000 0.003 0.003 -0.003 0.005 0.004 -0.011
## 26 std__WJ10_z~1 26 -0.002 0.004 0.004 -0.007 0.003 0.002 -0.012
## 27 std__DS_z~1 27 -0.002 0.002 0.001 -0.005 0.000 -0.002 0.001
## 28 std__g~1 28 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## SD_nonlin
## 1 0.000
## 2 0.018
## 3 0.007
## 4 0.005
## 5 0.008
## 6 0.011
## 7 0.018
## 8 0.007
## 9 0.002
## 10 0.002
## 11 0.001
## 12 0.002
## 13 0.001
## 14 0.000
## 15 0.004
## 16 0.008
## 17 0.008
## 18 0.003
## 19 0.006
## 20 0.013
## 21 0.013
## 22 0.003
## 23 0.000
## 24 0.002
## 25 0.001
## 26 0.002
## 27 0.001
## 28 0.000
##
## Distribution of Moderator: Density and Effective Sample Size
##
## M=0.062 | SD=0.851
##
## moderator wgt Neff
## 1 -0.246 0.100 580.209
## 2 0.070 0.167 848.175
## 3 0.264 0.233 945.139
## 4 0.453 0.259 942.705
## 5 0.619 0.241 853.379
##
## variable M SD min max
## 1 moderator 0.062 0.851 -4.930 1.448
## 2 wgt 0.200 0.066 0.100 0.259
## 3 Neff 833.921 149.294 580.209 945.139
## NULL
# run invariance models for comparison
cat("\nScalar invariance at h = 1.0:\n")
##
## Scalar invariance at h = 1.0:
lsem_strong_h1 <- lsem.estimate(
data = as.data.frame(d_apa_cog),
moderator = "PC1A_z",
moderator.grid = grid_points,
lavmodel = lsem_model,
h = 1.0,
estimator = "ML",
meanstructure = TRUE,
standardized = TRUE,
par_invariant = strong_constraints
)
## Use joint estimation
## ** Residualize Data
## ** Fit lavaan model
## |*****|
## |lavoptions ... done.
## lavdata ... done.
## Number of observations per group:
## Group 1 260
## Group 2 469
## Group 3 568
## Group 4 592
## Group 5 526
## lavsamplestats ... done.
## lavh1 ... start:
## lavh1 ... done.
## lavpartable bounds ... done.
## lavstart ... done.
## lavmodel ... done.
## lavoptim ... start:
## attempt 1 -- default options
## objective function = 0.4754973947033595
## quasi-Newton steps using NLMINB:
## objective function = 0.4754973947033595
## objective function = 0.4323832378287140
## objective function = 0.3874701232351613
## objective function = 0.3482134483756066
## objective function = 0.2804427797398291
## objective function = Inf
## objective function = 0.2379622115766932
## objective function = 0.1730091282698487
## objective function = 0.1680302255292255
## objective function = 0.1567538520268486
## objective function = 0.1247766203257206
## objective function = 0.2271158877186178
## objective function = 0.1086158335287673
## objective function = 0.0957725969251728
## objective function = 0.1052869419504699
## objective function = 0.0745525163747240
## objective function = 0.0731171375826762
## objective function = 0.0554442546532254
## objective function = 0.0435283489168887
## objective function = 0.0347739538808998
## objective function = 0.0320784728835493
## objective function = 0.0230142477661597
## objective function = 0.0194795132922565
## objective function = 0.0175309096612321
## objective function = 0.0150897371076818
## objective function = 0.0136454140263130
## objective function = 0.0124600870458997
## objective function = 0.0115087286161267
## objective function = 0.0109267888074380
## objective function = 0.0107942172087586
## objective function = 0.0106899229950802
## objective function = 0.0106549826098450
## objective function = 0.0106364690668438
## objective function = 0.0106296204274507
## objective function = 0.0106267595722763
## objective function = 0.0106251438164016
## objective function = 0.0106239152061211
## objective function = 0.0106235657222502
## objective function = 0.0106234409695969
## objective function = 0.0106233851458748
## objective function = 0.0106233491320922
## objective function = 0.0106233278272785
## objective function = 0.0106233164413848
## objective function = 0.0106233122049660
## objective function = 0.0106233107152364
## objective function = 0.0106233100159490
## objective function = 0.0106233096438381
## objective function = 0.0106233094940255
## objective function = 0.0106233094441505
## objective function = 0.0106233094262164
## objective function = 0.0106233094179794
## objective function = 0.0106233094140130
## objective function = 0.0106233094124133
## objective function = 0.0106233094119274
## objective function = 0.0106233094119274
## convergence status (0=ok): 0
## nlminb message says: relative convergence (4)
## number of iterations: 47
## number of function evaluations [objective, gradient]: 54 48
## lavoptim ... done.
## lavimplied ... done.
## lavloglik ... done.
## computing TEST for test(s) = standard ... done.
## lavbaseline ... done.
## post check ... done.
##
## Use lsem.bootstrap() for computing standard errors!
##
## -|
## ** Parameter summary
print(summary(lsem_strong_h1)$fitstat_joint)
## -----------------------------------------------------------------
## Local Structural Equation Model
##
## sirt 4.2-133 (2025-09-27 12:57:51)
## lavaan 0.6-19 (2024-09-26 15:50:14 UTC)
##
## R version 4.5.1 (2025-06-13 ucrt) x86_64, mingw32 | nodename=HUANG | login=mh198
##
## Function 'sirt::lsem.estimate', type='LSEM'
##
##
## Call:
## lsem.estimate(data = as.data.frame(d_apa_cog), moderator = "PC1A_z",
## moderator.grid = grid_points, lavmodel = lsem_model, h = 1,
## standardized = TRUE, par_invariant = strong_constraints,
## estimator = "ML", meanstructure = TRUE)
##
## Date of Analysis: 2026-06-12 15:04:41.856415
## Time difference of 5.148766 secs
## Computation Time: 5.148766
##
## Number of observations in datasets = 1538
## Used observations in analysis = 1524
## Used sampling weights: FALSE
## Bandwidth factor = 1
## Bandwidth = 0.196
## Number of focal points for moderator = 5
##
## Used joint estimation: TRUE
## Used sufficient statistics: TRUE
## Used local linear smoothing: TRUE
## Used pseudo weights: FALSE
## Used lavaan package: TRUE
## Used lavaan.survey package: FALSE
##
## Mean structure modelled: TRUE
##
## lavaan Model
## g =~ PPVT_z + WJ9_z + WJ10_z + DS_z
##
## Global Fit Statistics for Joint Estimation
##
## stat value
## 1 rmsea 0.027
## 2 cfi 0.995
## 3 tli 0.996
## 4 gfi 0.990
## 5 srmr 0.026
##
##
## Parameter Estimate Summary
##
## par parindex M SD MAD Min Max lin_int
## 1 g=~PPVT_z 1 1.000 0.000 0.000 1.000 1.000 1.000
## 2 g=~WJ9_z 2 1.307 0.000 0.000 1.307 1.307 1.307
## 3 g=~WJ10_z 3 1.250 0.000 0.000 1.250 1.250 1.250
## 4 g=~DS_z 4 0.873 0.000 0.000 0.873 0.873 0.873
## 5 PPVT_z~~PPVT_z 5 0.564 0.088 0.075 0.478 0.720 0.666
## 6 WJ9_z~~WJ9_z 6 0.317 0.031 0.028 0.275 0.364 0.346
## 7 WJ10_z~~WJ10_z 7 0.409 0.013 0.011 0.394 0.431 0.403
## 8 DS_z~~DS_z 8 0.723 0.036 0.034 0.669 0.762 0.684
## 9 g~~g 9 0.398 0.005 0.004 0.387 0.404 0.395
## 10 PPVT_z~1 10 0.002 0.000 0.000 0.002 0.002 0.002
## 11 WJ9_z~1 11 0.002 0.000 0.000 0.002 0.002 0.002
## 12 WJ10_z~1 12 0.000 0.000 0.000 0.000 0.000 0.000
## 13 DS_z~1 13 -0.001 0.000 0.000 -0.001 -0.001 -0.001
## 14 g~1 14 0.000 0.000 0.000 0.000 0.000 0.000
## 15 std__g=~PPVT_z 15 0.645 0.030 0.026 0.591 0.673 0.611
## 16 std__g=~WJ9_z 16 0.826 0.013 0.012 0.807 0.843 0.813
## 17 std__g=~WJ10_z 17 0.777 0.004 0.003 0.771 0.780 0.778
## 18 std__g=~DS_z 18 0.544 0.009 0.008 0.535 0.557 0.553
## 19 std__PPVT_z~~PPVT_z 19 0.583 0.038 0.033 0.547 0.651 0.626
## 20 std__WJ9_z~~WJ9_z 20 0.317 0.021 0.020 0.289 0.349 0.339
## 21 std__WJ10_z~~WJ10_z 21 0.397 0.006 0.005 0.391 0.406 0.395
## 22 std__DS_z~~DS_z 22 0.704 0.009 0.009 0.690 0.713 0.694
## 23 std__g~~g 23 1.000 0.000 0.000 1.000 1.000 1.000
## 24 std__PPVT_z~1 24 0.002 0.000 0.000 0.001 0.002 0.001
## 25 std__WJ9_z~1 25 0.002 0.000 0.000 0.002 0.002 0.002
## 26 std__WJ10_z~1 26 0.000 0.000 0.000 0.000 0.000 0.000
## 27 std__DS_z~1 27 -0.001 0.000 0.000 -0.001 -0.001 -0.001
## 28 std__g~1 28 0.000 0.000 0.000 0.000 0.000 0.000
## lin_slo SD_nonlin
## 1 0.000 0.000
## 2 0.000 0.000
## 3 0.000 0.000
## 4 0.000 0.000
## 5 -0.323 0.022
## 6 -0.090 0.020
## 7 0.019 0.012
## 8 0.124 0.015
## 9 0.010 0.004
## 10 0.000 0.000
## 11 0.000 0.000
## 12 0.000 0.000
## 13 0.000 0.000
## 14 0.000 0.000
## 15 0.109 0.007
## 16 0.041 0.007
## 17 -0.003 0.003
## 18 -0.029 0.004
## 19 -0.139 0.009
## 20 -0.068 0.012
## 21 0.005 0.005
## 22 0.031 0.004
## 23 0.000 0.000
## 24 0.000 0.000
## 25 0.000 0.000
## 26 0.000 0.000
## 27 0.000 0.000
## 28 0.000 0.000
##
## Distribution of Moderator: Density and Effective Sample Size
##
## M=0.062 | SD=0.851
##
## moderator wgt Neff
## 1 -0.246 0.100 262.694
## 2 0.070 0.167 473.933
## 3 0.264 0.233 572.817
## 4 0.453 0.259 597.446
## 5 0.619 0.241 531.044
##
## variable M SD min max
## 1 moderator 0.062 0.851 -4.930 1.448
## 2 wgt 0.200 0.066 0.100 0.259
## 3 Neff 487.587 134.150 262.694 597.446
## NULL
cat("\nScalar invariance at h = 2.0:\n")
##
## Scalar invariance at h = 2.0:
lsem_strong_h2 <- lsem.estimate(
data = as.data.frame(d_apa_cog),
moderator = "PC1A_z",
moderator.grid = grid_points,
lavmodel = lsem_model,
h = 2.0,
estimator = "ML",
meanstructure = TRUE,
standardized = TRUE,
par_invariant = strong_constraints
)
## Use joint estimation
## ** Residualize Data
## ** Fit lavaan model
## |*****|
## |lavoptions ... done.
## lavdata ... done.
## Number of observations per group:
## Group 1 575
## Group 2 840
## Group 3 936
## Group 4 934
## Group 5 846
## lavsamplestats ... done.
## lavh1 ... start:
## lavh1 ... done.
## lavpartable bounds ... done.
## lavstart ... done.
## lavmodel ... done.
## lavoptim ... start:
## attempt 1 -- default options
## objective function = 0.0363330532126398
## quasi-Newton steps using NLMINB:
## objective function = 0.0363330532126398
## objective function = 0.0277690480209134
## objective function = 0.0163251982990948
## objective function = 0.0142496399345276
## objective function = 0.0128700138865049
## objective function = 0.0120467045640627
## objective function = 0.0106292798373284
## objective function = 0.0099049490033018
## objective function = 0.0092688999741999
## objective function = 0.0087938395377918
## objective function = 0.0085739861569305
## objective function = 0.0084008216551196
## objective function = 0.0082741335982951
## objective function = 0.0081913254427921
## objective function = 0.0081527112279510
## objective function = 0.0081361921129442
## objective function = 0.0081284871198465
## objective function = 0.0081248618935237
## objective function = 0.0081232184407308
## objective function = 0.0081224615251444
## objective function = 0.0081220217220510
## objective function = 0.0081217527192213
## objective function = 0.0081216095271687
## objective function = 0.0081215389270073
## objective function = 0.0081215085968742
## objective function = 0.0081214992573260
## objective function = 0.0081214971770618
## objective function = 0.0081214966605514
## objective function = 0.0081214964828304
## objective function = 0.0081214964300407
## objective function = 0.0081214964161902
## objective function = 0.0081214964114861
## objective function = 0.0081214964095177
## objective function = 0.0081214964087900
## objective function = 0.0081214964087900
## convergence status (0=ok): 0
## nlminb message says: relative convergence (4)
## number of iterations: 33
## number of function evaluations [objective, gradient]: 34 34
## lavoptim ... done.
## lavimplied ... done.
## lavloglik ... done.
## computing TEST for test(s) = standard ... done.
## lavbaseline ... done.
## post check ... done.
##
## Use lsem.bootstrap() for computing standard errors!
##
## -|
## ** Parameter summary
print(summary(lsem_strong_h2)$fitstat_joint)
## -----------------------------------------------------------------
## Local Structural Equation Model
##
## sirt 4.2-133 (2025-09-27 12:57:51)
## lavaan 0.6-19 (2024-09-26 15:50:14 UTC)
##
## R version 4.5.1 (2025-06-13 ucrt) x86_64, mingw32 | nodename=HUANG | login=mh198
##
## Function 'sirt::lsem.estimate', type='LSEM'
##
##
## Call:
## lsem.estimate(data = as.data.frame(d_apa_cog), moderator = "PC1A_z",
## moderator.grid = grid_points, lavmodel = lsem_model, h = 2,
## standardized = TRUE, par_invariant = strong_constraints,
## estimator = "ML", meanstructure = TRUE)
##
## Date of Analysis: 2026-06-12 15:04:46.372483
## Time difference of 4.361023 secs
## Computation Time: 4.361023
##
## Number of observations in datasets = 1538
## Used observations in analysis = 1524
## Used sampling weights: FALSE
## Bandwidth factor = 2
## Bandwidth = 0.392
## Number of focal points for moderator = 5
##
## Used joint estimation: TRUE
## Used sufficient statistics: TRUE
## Used local linear smoothing: TRUE
## Used pseudo weights: FALSE
## Used lavaan package: TRUE
## Used lavaan.survey package: FALSE
##
## Mean structure modelled: TRUE
##
## lavaan Model
## g =~ PPVT_z + WJ9_z + WJ10_z + DS_z
##
## Global Fit Statistics for Joint Estimation
##
## stat value
## 1 rmsea 0.030
## 2 cfi 0.994
## 3 tli 0.995
## 4 gfi 0.992
## 5 srmr 0.023
##
##
## Parameter Estimate Summary
##
## par parindex M SD MAD Min Max lin_int
## 1 g=~PPVT_z 1 1.000 0.000 0.000 1.000 1.000 1.000
## 2 g=~WJ9_z 2 1.309 0.000 0.000 1.309 1.309 1.309
## 3 g=~WJ10_z 3 1.255 0.000 0.000 1.255 1.255 1.255
## 4 g=~DS_z 4 0.873 0.000 0.000 0.873 0.873 0.873
## 5 PPVT_z~~PPVT_z 5 0.570 0.070 0.060 0.493 0.708 0.654
## 6 WJ9_z~~WJ9_z 6 0.326 0.035 0.028 0.288 0.406 0.368
## 7 WJ10_z~~WJ10_z 7 0.398 0.013 0.011 0.364 0.409 0.392
## 8 DS_z~~DS_z 8 0.715 0.018 0.015 0.674 0.731 0.694
## 9 g~~g 9 0.396 0.004 0.002 0.386 0.398 0.393
## 10 PPVT_z~1 10 -0.002 0.000 0.000 -0.002 -0.002 -0.002
## 11 WJ9_z~1 11 0.001 0.000 0.000 0.001 0.001 0.001
## 12 WJ10_z~1 12 -0.002 0.000 0.000 -0.002 -0.002 -0.002
## 13 DS_z~1 13 -0.002 0.000 0.000 -0.002 -0.002 -0.002
## 14 g~1 14 0.000 0.000 0.000 0.000 0.000 0.000
## 15 std__g=~PPVT_z 15 0.641 0.024 0.020 0.594 0.667 0.613
## 16 std__g=~WJ9_z 16 0.822 0.015 0.012 0.787 0.838 0.804
## 17 std__g=~WJ10_z 17 0.781 0.004 0.003 0.778 0.791 0.783
## 18 std__g=~DS_z 18 0.545 0.004 0.003 0.541 0.551 0.549
## 19 std__PPVT_z~~PPVT_z 19 0.588 0.030 0.026 0.555 0.648 0.624
## 20 std__WJ9_z~~WJ9_z 20 0.324 0.024 0.019 0.298 0.381 0.353
## 21 std__WJ10_z~~WJ10_z 21 0.389 0.006 0.005 0.375 0.395 0.387
## 22 std__DS_z~~DS_z 22 0.703 0.004 0.004 0.697 0.707 0.699
## 23 std__g~~g 23 1.000 0.000 0.000 1.000 1.000 1.000
## 24 std__PPVT_z~1 24 -0.002 0.000 0.000 -0.002 -0.002 -0.002
## 25 std__WJ9_z~1 25 0.001 0.000 0.000 0.001 0.001 0.001
## 26 std__WJ10_z~1 26 -0.002 0.000 0.000 -0.002 -0.002 -0.002
## 27 std__DS_z~1 27 -0.002 0.000 0.000 -0.002 -0.002 -0.002
## 28 std__g~1 28 0.000 0.000 0.000 0.000 0.000 0.000
## lin_slo SD_nonlin
## 1 0.000 0.000
## 2 0.000 0.000
## 3 0.000 0.000
## 4 0.000 0.000
## 5 -0.264 0.008
## 6 -0.132 0.003
## 7 0.018 0.012
## 8 0.065 0.006
## 9 0.008 0.003
## 10 0.000 0.000
## 11 0.000 0.000
## 12 0.000 0.000
## 13 0.000 0.000
## 14 0.000 0.000
## 15 0.089 0.002
## 16 0.056 0.002
## 17 -0.004 0.004
## 18 -0.014 0.001
## 19 -0.113 0.003
## 20 -0.091 0.003
## 21 0.006 0.006
## 22 0.015 0.001
## 23 0.000 0.000
## 24 0.000 0.000
## 25 0.000 0.000
## 26 0.000 0.000
## 27 0.000 0.000
## 28 0.000 0.000
##
## Distribution of Moderator: Density and Effective Sample Size
##
## M=0.062 | SD=0.851
##
## moderator wgt Neff
## 1 -0.246 0.100 580.209
## 2 0.070 0.167 848.175
## 3 0.264 0.233 945.139
## 4 0.453 0.259 942.705
## 5 0.619 0.241 853.379
##
## variable M SD min max
## 1 moderator 0.062 0.851 -4.930 1.448
## 2 wgt 0.200 0.066 0.100 0.259
## 3 Neff 833.921 149.294 580.209 945.139
## NULL
# ==============================================================
# 12. Helper: WLS with Within-Model Standardization
# ==============================================================
run_model <- function(formula, data, label) {
all_vars <- all.vars(formula)
var_map <- c(
"g" = "g_raw",
"PC1A_z" = "K5PC1A",
"age_c" = "age",
"ave_WAIS_res_z" = "ave_WAIS_raw",
"M_WAIS_z" = "M_WAIS_res",
"F_WAIS_z" = "F_WAIS_res",
"SES" = "SES_raw"
)
raw_vars <- all_vars
for (i in seq_along(raw_vars)) {
if (raw_vars[i] %in% names(var_map)) {
raw_vars[i] <- var_map[raw_vars[i]]
}
}
d_complete <- data[complete.cases(data[, raw_vars]) & !is.na(data$wt), ]
if ("PC1A_z" %in% all_vars) d_complete$PC1A_z <- as.vector(scale(d_complete$K5PC1A))
if ("g" %in% all_vars) d_complete$g <- as.vector(scale(d_complete$g_raw))
if ("age_c" %in% all_vars) d_complete$age_c <- as.vector(scale(d_complete$age))
if ("ave_WAIS_res_z" %in% all_vars) d_complete$ave_WAIS_res_z <- as.vector(scale(d_complete$ave_WAIS_raw))
if ("M_WAIS_z" %in% all_vars) d_complete$M_WAIS_z <- as.vector(scale(d_complete$M_WAIS_res))
if ("F_WAIS_z" %in% all_vars) d_complete$F_WAIS_z <- as.vector(scale(d_complete$F_WAIS_res))
if ("SES" %in% all_vars) d_complete$SES <- as.vector(scale(d_complete$SES_raw))
for (v in c("PPVT", "WJ9", "WJ10")) {
if (v %in% all_vars && v %in% names(d_complete)) {
d_complete[[v]] <- as.vector(scale(d_complete[[v]]))
}
}
model <- lm(formula, data = d_complete, weights = wt)
robust <- coeftest(model, vcov = vcovHC(model, type = "HC1"))
cat("\n", paste(rep("=", 80), collapse = ""), "\n")
cat(label, "\n")
cat("R2 =", round(summary(model)$r.squared, 4),
" Adj-R2 =", round(summary(model)$adj.r.squared, 4), "\n")
cat(paste(rep("=", 80), collapse = ""), "\n")
cat(sprintf("%20s %9s %8s %8s %10s\n", "Variable", "B", "SE", "t", "p"))
cat(paste(rep("-", 60), collapse = ""), "\n")
for (i in 1:nrow(robust)) {
var_name <- rownames(robust)[i]
b <- robust[i, 1]; se <- robust[i, 2]; t <- robust[i, 3]; p <- robust[i, 4]
p_str <- ifelse(p >= 0.0001, sprintf("%.4f", p), sprintf("%.2e", p))
cat(sprintf("%20s %+9.4f %8.4f %+8.3f %10s\n", var_name, b, se, t, p_str))
}
cat("N =", nrow(d_complete), "\n")
invisible(model)
}
# ==============================================================
# 13. APA: Child g Models
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("AFRO-PAN-AMERICAN: Child g as DV\n")
## AFRO-PAN-AMERICAN: Child g as DV
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
run_model(g ~ PC1A_z + sire + gen_2nd + gen_3rd + sex + age_c, d_apa,
"APA M1: PC1A + SIRE + Generation + Demographics")
##
## ================================================================================
## APA M1: PC1A + SIRE + Generation + Demographics
## R2 = 0.0194 Adj-R2 = 0.0142
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) +0.0061 0.0553 +0.110 0.9121
## PC1A_z -0.1298 0.0474 -2.737 0.0063
## sireBW -0.3055 0.2117 -1.443 0.1491
## sireBH +0.1142 0.1866 +0.612 0.5407
## sireBO +0.0059 0.1435 +0.041 0.9673
## gen_2ndTRUE +0.4220 0.2193 +1.925 0.0545
## gen_3rdTRUE +0.2412 0.2361 +1.022 0.3070
## sex -0.0111 0.0679 -0.164 0.8699
## age_c -0.0178 0.0305 -0.584 0.5591
## N = 1538
run_model(g ~ PC1A_z + sire + gen_2nd + gen_3rd + sex + age_c + ave_WAIS_res_z, d_apa,
"APA M2: + Parental WAIS")
##
## ================================================================================
## APA M2: + Parental WAIS
## R2 = 0.0396 Adj-R2 = 0.0336
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) +0.0115 0.0576 +0.199 0.8424
## PC1A_z -0.1041 0.0480 -2.166 0.0304
## sireBW -0.3094 0.2209 -1.401 0.1616
## sireBH +0.1641 0.2119 +0.774 0.4390
## sireBO +0.0945 0.1384 +0.683 0.4946
## gen_2ndTRUE +0.4598 0.3051 +1.507 0.1320
## gen_3rdTRUE +0.3124 0.2421 +1.290 0.1971
## sex -0.0172 0.0712 -0.242 0.8091
## age_c -0.0206 0.0313 -0.659 0.5100
## ave_WAIS_res_z +0.1514 0.0345 +4.391 1.21e-05
## N = 1450
run_model(g ~ PC1A_z + sire + gen_2nd + gen_3rd + sex + age_c + ave_WAIS_res_z + single + SES, d_apa,
"APA M3: + SES + Single Parent")
##
## ================================================================================
## APA M3: + SES + Single Parent
## R2 = 0.1366 Adj-R2 = 0.13
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0032 0.0698 -0.046 0.9635
## PC1A_z -0.0812 0.0470 -1.729 0.0841
## sireBW -0.1897 0.2131 -0.890 0.3735
## sireBH +0.3436 0.1946 +1.765 0.0777
## sireBO +0.1111 0.1284 +0.865 0.3873
## gen_2ndTRUE +0.2863 0.2525 +1.134 0.2571
## gen_3rdTRUE +0.2614 0.2659 +0.983 0.3259
## sex +0.0054 0.0671 +0.081 0.9354
## age_c -0.0052 0.0296 -0.174 0.8618
## ave_WAIS_res_z +0.0695 0.0345 +2.011 0.0445
## single -0.0461 0.0642 -0.719 0.4724
## SES +0.3185 0.0311 +10.228 9.48e-24
## N = 1450
# ==============================================================
# 13b. Unweighted Model Sensitivity
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("ROBUSTNESS: Unweighted models\n")
## ROBUSTNESS: Unweighted models
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
d_uw <- d_apa[complete.cases(d_apa[, c("g_raw", "K5PC1A", "sire",
"gen_2nd", "gen_3rd", "sex", "age")]), ]
d_uw$PC1A_z <- as.vector(scale(d_uw$K5PC1A))
d_uw$g <- as.vector(scale(d_uw$g_raw))
d_uw$age_c <- as.vector(scale(d_uw$age))
m_uw <- lm(g ~ PC1A_z + sire + gen_2nd + gen_3rd + sex + age_c, data = d_uw)
robust_uw <- coeftest(m_uw, vcov = vcovHC(m_uw, type = "HC1"))
cat("APA unweighted M1: PC1A coefficient:\n")
## APA unweighted M1: PC1A coefficient:
print(robust_uw["PC1A_z", ])
## Estimate Std. Error t value Pr(>|t|)
## -0.119666159 0.037005695 -3.233722753 0.001247956
cat("N =", nrow(d_uw), "\n")
## N = 1538
# ==============================================================
# 14. APA: Average Parental WAIS as DV (excluding BW)
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("APA (excl BW): Parental WAIS as DV\n")
## APA (excl BW): Parental WAIS as DV
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
run_model(ave_WAIS_res_z ~ PC1A_z + sire + gen_2nd + gen_3rd, d_apa_noBW,
"APA WAIS M1: PC1A + SIRE + Generation")
##
## ================================================================================
## APA WAIS M1: PC1A + SIRE + Generation
## R2 = 0.0151 Adj-R2 = 0.0116
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) +0.0558 0.0392 +1.423 0.1549
## PC1A_z -0.1165 0.0358 -3.259 0.0011
## sireBH -0.2848 0.1385 -2.057 0.0399
## sireBO -0.2225 0.1519 -1.465 0.1433
## gen_2ndTRUE +0.0986 0.1309 +0.754 0.4512
## gen_3rdTRUE -0.4637 0.2154 -2.153 0.0315
## N = 1393
# ==============================================================
# 15. APA: Mother & Father WAIS separately (full APA, parent self-ID)
# ==============================================================
# Use parent's own racial self-identification, not couple SIRE.
# Full APA sample (including BW) — child ancestry is the predictor,
# so both parents' WAIS scores are informative in their own models.
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("APA: Mother & Father WAIS separately (full sample, parent self-ID)\n")
## APA: Mother & Father WAIS separately (full sample, parent self-ID)
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
run_model(M_WAIS_z ~ PC1A_z + mother_black_sire + gen_2nd + gen_3rd, d_apa,
"APA Mother WAIS M1 (full sample, parent self-ID)")
##
## ================================================================================
## APA Mother WAIS M1 (full sample, parent self-ID)
## R2 = 0.0102 Adj-R2 = 0.0074
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) +0.0564 0.1570 +0.359 0.7194
## PC1A_z -0.1002 0.0494 -2.028 0.0428
## mother_black_sire -0.0273 0.1663 -0.164 0.8697
## gen_2ndTRUE +0.0956 0.2313 +0.413 0.6794
## gen_3rdTRUE -0.0389 0.3212 -0.121 0.9035
## N = 1410
run_model(F_WAIS_z ~ PC1A_z + father_black_sire + gen_2nd + gen_3rd, d_apa,
"APA Father WAIS M1 (full sample, parent self-ID)")
##
## ================================================================================
## APA Father WAIS M1 (full sample, parent self-ID)
## R2 = 0.0116 Adj-R2 = 0.0079
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) +0.0502 0.1120 +0.448 0.6541
## PC1A_z -0.0816 0.0315 -2.590 0.0097
## father_black_sire -0.0259 0.1179 -0.220 0.8260
## gen_2ndTRUE -0.0114 0.1315 -0.086 0.9311
## gen_3rdTRUE -0.7400 0.3857 -1.919 0.0553
## N = 1076
# ==============================================================
# 16. APA: SES as DV (excluding BW)
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("APA (excl BW): SES as DV\n")
## APA (excl BW): SES as DV
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
run_model(SES ~ PC1A_z + sire + gen_2nd + gen_3rd, d_apa_noBW,
"APA SES M1")
##
## ================================================================================
## APA SES M1
## R2 = 0.0228 Adj-R2 = 0.0195
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) +0.0701 0.0359 +1.953 0.0510
## PC1A_z -0.0980 0.0354 -2.768 0.0057
## sireBH -0.5652 0.1762 -3.208 0.0014
## sireBO -0.1239 0.1089 -1.138 0.2554
## gen_2ndTRUE +0.5087 0.1679 +3.030 0.0025
## gen_3rdTRUE +0.0188 0.1853 +0.101 0.9194
## N = 1477
# ==============================================================
# 17. APA: Individual Tests (Model 1 spec)
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("APA: Individual Tests (Model 1 spec)\n")
## APA: Individual Tests (Model 1 spec)
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
run_model(PPVT ~ PC1A_z + sire + gen_2nd + gen_3rd + sex + age_c, d_apa, "APA: PPVT")
##
## ================================================================================
## APA: PPVT
## R2 = 0.0243 Adj-R2 = 0.0192
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0165 0.0526 -0.314 0.7537
## PC1A_z -0.1618 0.0508 -3.188 0.0015
## sireBW -0.3483 0.2267 -1.537 0.1246
## sireBH +0.1403 0.2376 +0.590 0.5551
## sireBO -0.0264 0.1595 -0.166 0.8684
## gen_2ndTRUE +0.3827 0.2128 +1.799 0.0723
## gen_3rdTRUE +0.0645 0.3968 +0.163 0.8709
## sex +0.0210 0.0661 +0.319 0.7501
## age_c -0.0286 0.0316 -0.905 0.3658
## N = 1531
run_model(WJ9 ~ PC1A_z + sire + gen_2nd + gen_3rd + sex + age_c, d_apa, "APA: WJ Applied Problems")
##
## ================================================================================
## APA: WJ Applied Problems
## R2 = 0.0315 Adj-R2 = 0.0264
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.1509 0.0538 -2.803 0.0051
## PC1A_z -0.0791 0.0452 -1.750 0.0803
## sireBW -0.2980 0.2103 -1.417 0.1567
## sireBH -0.0106 0.1556 -0.068 0.9456
## sireBO +0.0998 0.1559 +0.640 0.5221
## gen_2ndTRUE +0.2954 0.1849 +1.597 0.1104
## gen_3rdTRUE +0.1111 0.1836 +0.605 0.5453
## sex +0.3031 0.0641 +4.731 2.44e-06
## age_c -0.0496 0.0306 -1.618 0.1059
## N = 1523
run_model(WJ10 ~ PC1A_z + sire + gen_2nd + gen_3rd + sex + age_c, d_apa, "APA: WJ Letter-Word")
##
## ================================================================================
## APA: WJ Letter-Word
## R2 = 0.0222 Adj-R2 = 0.0171
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0372 0.0569 -0.655 0.5128
## PC1A_z -0.1014 0.0438 -2.315 0.0208
## sireBW -0.1976 0.1991 -0.992 0.3211
## sireBH +0.0681 0.1506 +0.452 0.6511
## sireBO -0.0490 0.1343 -0.365 0.7151
## gen_2ndTRUE +0.3742 0.2112 +1.772 0.0766
## gen_3rdTRUE +0.4718 0.1862 +2.533 0.0114
## sex +0.1105 0.0686 +1.612 0.1072
## age_c -0.0881 0.0293 -3.004 0.0027
## N = 1524
# ==============================================================
# 18. HAA: Child g Models
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("HISTORICAL AFRICAN AMERICAN: Child g as DV\n")
## HISTORICAL AFRICAN AMERICAN: Child g as DV
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
run_model(g ~ PC1A_z + sex + age_c, d_haa, "HAA M1: PC1A + Demographics")
##
## ================================================================================
## HAA M1: PC1A + Demographics
## R2 = 0.009 Adj-R2 = 0.0066
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) +0.0034 0.0589 +0.058 0.9540
## PC1A_z -0.1012 0.0317 -3.193 0.0014
## sex +0.0263 0.0727 +0.362 0.7177
## age_c +0.0062 0.0332 +0.185 0.8529
## N = 1287
run_model(g ~ PC1A_z + sex + age_c + ave_WAIS_res_z, d_haa, "HAA M2: + Parental WAIS")
##
## ================================================================================
## HAA M2: + Parental WAIS
## R2 = 0.036 Adj-R2 = 0.0329
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0136 0.0613 -0.221 0.8252
## PC1A_z -0.0934 0.0327 -2.858 0.0043
## sex +0.0559 0.0757 +0.738 0.4607
## age_c +0.0032 0.0334 +0.097 0.9226
## ave_WAIS_res_z +0.1675 0.0368 +4.548 5.95e-06
## N = 1234
run_model(g ~ PC1A_z + sex + age_c + ave_WAIS_res_z + single + SES, d_haa, "HAA M3: + SES + Single")
##
## ================================================================================
## HAA M3: + SES + Single
## R2 = 0.1217 Adj-R2 = 0.1174
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0062 0.0775 -0.080 0.9359
## PC1A_z -0.0786 0.0323 -2.434 0.0151
## sex +0.0605 0.0731 +0.827 0.4081
## age_c +0.0187 0.0325 +0.574 0.5661
## ave_WAIS_res_z +0.0877 0.0365 +2.403 0.0164
## single -0.0642 0.0710 -0.904 0.3659
## SES +0.2925 0.0300 +9.764 9.69e-22
## N = 1234
# ==============================================================
# 19. HAA: Parental WAIS as DV
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("HAA: Parental WAIS as DV\n")
## HAA: Parental WAIS as DV
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
run_model(ave_WAIS_res_z ~ PC1A_z, d_haa, "HAA WAIS M1")
##
## ================================================================================
## HAA WAIS M1
## R2 = 0.0091 Adj-R2 = 0.0083
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) +0.0456 0.0399 +1.144 0.2530
## PC1A_z -0.0979 0.0389 -2.514 0.0121
## N = 1238
# ==============================================================
# 20. HAA: SES as DV
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("HAA: SES as DV\n")
## HAA: SES as DV
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
run_model(SES ~ PC1A_z, d_haa, "HAA SES M1")
##
## ================================================================================
## HAA SES M1
## R2 = 0.0062 Adj-R2 = 0.0055
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) +0.0787 0.0371 +2.122 0.0340
## PC1A_z -0.0836 0.0400 -2.090 0.0368
## N = 1292
# ==============================================================
# 21. HAA: Individual Tests (Model 1 spec)
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("HAA: Individual Tests (Model 1 spec)\n")
## HAA: Individual Tests (Model 1 spec)
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
run_model(PPVT ~ PC1A_z + sex + age_c, d_haa, "HAA: PPVT")
##
## ================================================================================
## HAA: PPVT
## R2 = 0.0158 Adj-R2 = 0.0135
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0442 0.0558 -0.791 0.4290
## PC1A_z -0.1229 0.0317 -3.873 0.0001
## sex +0.1079 0.0719 +1.501 0.1337
## age_c -0.0139 0.0347 -0.400 0.6893
## N = 1281
run_model(WJ9 ~ PC1A_z + sex + age_c, d_haa, "HAA: WJ Applied Problems")
##
## ================================================================================
## HAA: WJ Applied Problems
## R2 = 0.0233 Adj-R2 = 0.021
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.1436 0.0565 -2.544 0.0111
## PC1A_z -0.0578 0.0321 -1.803 0.0717
## sex +0.3011 0.0689 +4.372 1.33e-05
## age_c -0.0193 0.0332 -0.579 0.5625
## N = 1273
run_model(WJ10 ~ PC1A_z + sex + age_c, d_haa, "HAA: WJ Letter-Word")
##
## ================================================================================
## HAA: WJ Letter-Word
## R2 = 0.013 Adj-R2 = 0.0107
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0305 0.0605 -0.504 0.6145
## PC1A_z -0.0788 0.0311 -2.534 0.0114
## sex +0.1241 0.0740 +1.677 0.0937
## age_c -0.0770 0.0317 -2.426 0.0154
## N = 1277
# ==============================================================
# 21b. Mediation Analyses
# ==============================================================
# Mediation 1: PC1A -> SES -> child g (SES mediates ancestry-cognition)
# Mediation 2: PC1A -> parental WAIS -> parental SES
# (Mediation 3: PC1A -> parental WAIS -> child g is in Section 21c)
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("MEDIATION ANALYSES\n")
## MEDIATION ANALYSES
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
# Mediation 1: PC1A -> SES -> child g (APA excluding BW)
d_med1_apa <- d_apa_noBW[complete.cases(d_apa_noBW[, c(
"g_raw", "K5PC1A", "SES_raw", "sire", "gen_2nd", "gen_3rd", "sex", "age"
)]) & !is.na(d_apa_noBW$wt), ]
d_med1_apa$PC1A_z <- as.vector(scale(d_med1_apa$K5PC1A))
d_med1_apa$g <- as.vector(scale(d_med1_apa$g_raw))
d_med1_apa$SES <- as.vector(scale(d_med1_apa$SES_raw))
d_med1_apa$age_c <- as.vector(scale(d_med1_apa$age))
# M: PC1A -> SES
med_M1_apa <- lm(SES ~ PC1A_z + sire + gen_2nd + gen_3rd + sex + age_c,
data = d_med1_apa, weights = wt)
# Y: SES -> g, controlling for PC1A
out_Y1_apa <- lm(g ~ SES + PC1A_z + sire + gen_2nd + gen_3rd + sex + age_c,
data = d_med1_apa, weights = wt)
set.seed(11111)
med1_apa <- mediate(med_M1_apa, out_Y1_apa, treat = "PC1A_z",
mediator = "SES", boot = TRUE, sims = 1000)
## Running nonparametric bootstrap
summary(med1_apa)
##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.0312264 -0.0571902 -0.0067020 0.012 *
## ADE -0.0619238 -0.1256069 0.0019901 0.062 .
## Total Effect -0.0931502 -0.1610801 -0.0261930 0.008 **
## Prop. Mediated 0.3352262 0.0785313 1.0107897 0.016 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 1472
##
##
## Simulations: 1000
# Mediation 1: PC1A -> SES -> child g (HAA)
d_med1_haa <- d_haa[complete.cases(d_haa[, c(
"g_raw", "K5PC1A", "SES_raw", "sex", "age"
)]) & !is.na(d_haa$wt), ]
d_med1_haa$PC1A_z <- as.vector(scale(d_med1_haa$K5PC1A))
d_med1_haa$g <- as.vector(scale(d_med1_haa$g_raw))
d_med1_haa$SES <- as.vector(scale(d_med1_haa$SES_raw))
d_med1_haa$age_c <- as.vector(scale(d_med1_haa$age))
med_M1_haa <- lm(SES ~ PC1A_z + sex + age_c,
data = d_med1_haa, weights = wt)
out_Y1_haa <- lm(g ~ SES + PC1A_z + sex + age_c,
data = d_med1_haa, weights = wt)
set.seed(11111)
med1_haa <- mediate(med_M1_haa, out_Y1_haa, treat = "PC1A_z",
mediator = "SES", boot = TRUE, sims = 1000)
## Running nonparametric bootstrap
summary(med1_haa)
##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.0251325 -0.0525039 0.0030279 0.074 .
## ADE -0.0760313 -0.1349419 -0.0130390 0.016 *
## Total Effect -0.1011638 -0.1640583 -0.0368243 <2e-16 ***
## Prop. Mediated 0.2484338 -0.0354906 0.7008143 0.074 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 1287
##
##
## Simulations: 1000
# Mediation 2: PC1A -> parental WAIS -> parental SES (APA excluding BW)
d_med2_apa <- d_apa_noBW[complete.cases(d_apa_noBW[, c(
"ave_WAIS_raw", "K5PC1A", "SES_raw", "sire", "gen_2nd", "gen_3rd", "sex", "age"
)]) & !is.na(d_apa_noBW$wt), ]
d_med2_apa$PC1A_z <- as.vector(scale(d_med2_apa$K5PC1A))
d_med2_apa$SES <- as.vector(scale(d_med2_apa$SES_raw))
d_med2_apa$W <- as.vector(scale(d_med2_apa$ave_WAIS_raw))
d_med2_apa$age_c <- as.vector(scale(d_med2_apa$age))
# M: PC1A -> parental WAIS
med_M2_apa <- lm(W ~ PC1A_z + sire + gen_2nd + gen_3rd + sex + age_c,
data = d_med2_apa, weights = wt)
# Y: parental WAIS -> SES, controlling for PC1A
out_Y2_apa <- lm(SES ~ W + PC1A_z + sire + gen_2nd + gen_3rd + sex + age_c,
data = d_med2_apa, weights = wt)
set.seed(22222)
med2_apa <- mediate(med_M2_apa, out_Y2_apa, treat = "PC1A_z",
mediator = "W", boot = TRUE, sims = 1000)
## Running nonparametric bootstrap
summary(med2_apa)
##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.0276413 -0.0463945 -0.0101489 0.006 **
## ADE -0.0487482 -0.1218704 0.0218145 0.168
## Total Effect -0.0763895 -0.1534454 -0.0017425 0.046 *
## Prop. Mediated 0.3618467 0.0102994 1.4462108 0.048 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 1392
##
##
## Simulations: 1000
# Mediation 2: PC1A -> parental WAIS -> parental SES (HAA)
d_med2_haa <- d_haa[complete.cases(d_haa[, c(
"ave_WAIS_raw", "K5PC1A", "SES_raw", "sex", "age"
)]) & !is.na(d_haa$wt), ]
d_med2_haa$PC1A_z <- as.vector(scale(d_med2_haa$K5PC1A))
d_med2_haa$SES <- as.vector(scale(d_med2_haa$SES_raw))
d_med2_haa$W <- as.vector(scale(d_med2_haa$ave_WAIS_raw))
d_med2_haa$age_c <- as.vector(scale(d_med2_haa$age))
med_M2_haa <- lm(W ~ PC1A_z + sex + age_c,
data = d_med2_haa, weights = wt)
out_Y2_haa <- lm(SES ~ W + PC1A_z + sex + age_c,
data = d_med2_haa, weights = wt)
set.seed(22222)
med2_haa <- mediate(med_M2_haa, out_Y2_haa, treat = "PC1A_z",
mediator = "W", boot = TRUE, sims = 1000)
## Running nonparametric bootstrap
summary(med2_haa)
##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.022874 -0.044549 -0.002494 0.022 *
## ADE -0.051226 -0.126610 0.035647 0.216
## Total Effect -0.074100 -0.153312 0.013584 0.086 .
## Prop. Mediated 0.308693 -0.766854 1.955594 0.104
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 1237
##
##
## Simulations: 1000
# ==============================================================
# 21c. Mediation: PC1A -> Parental WAIS -> Child g
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("MEDIATION: PC1A -> Parental WAIS -> Child g\n")
## MEDIATION: PC1A -> Parental WAIS -> Child g
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
# APA (exclude BW)
d_med_pw_apa <- d_apa_noBW[complete.cases(d_apa_noBW[, c(
"g_raw", "K5PC1A", "ave_WAIS_raw", "sire", "gen_2nd", "gen_3rd", "sex", "age"
)]) & !is.na(d_apa_noBW$wt), ]
d_med_pw_apa$PC1A_z <- as.vector(scale(d_med_pw_apa$K5PC1A))
d_med_pw_apa$g <- as.vector(scale(d_med_pw_apa$g_raw))
d_med_pw_apa$W <- as.vector(scale(d_med_pw_apa$ave_WAIS_raw))
d_med_pw_apa$age_c <- as.vector(scale(d_med_pw_apa$age))
# M: PC1A -> parental WAIS
med_M_pw_apa <- lm(W ~ PC1A_z + sire + gen_2nd + gen_3rd + sex + age_c,
data = d_med_pw_apa, weights = wt)
# Y: parental WAIS -> g, controlling for PC1A
out_Y_pw_apa <- lm(g ~ W + PC1A_z + sire + gen_2nd + gen_3rd + sex + age_c,
data = d_med_pw_apa, weights = wt)
set.seed(88888)
med_pw_apa <- mediate(med_M_pw_apa, out_Y_pw_apa, treat = "PC1A_z",
mediator = "W", boot = TRUE, sims = 1000)
## Running nonparametric bootstrap
cat("\nAPA (excl BW): PC1A -> Parental WAIS -> Child g\n")
##
## APA (excl BW): PC1A -> Parental WAIS -> Child g
summary(med_pw_apa)
##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.0168691 -0.0319366 -0.0050132 0.002 **
## ADE -0.0727504 -0.1380107 -0.0035731 0.032 *
## Total Effect -0.0896195 -0.1572542 -0.0202900 0.016 *
## Prop. Mediated 0.1882302 0.0525785 0.7121987 0.018 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 1389
##
##
## Simulations: 1000
# HAA
d_med_pw_haa <- d_haa[complete.cases(d_haa[, c(
"g_raw", "K5PC1A", "ave_WAIS_raw", "sex", "age"
)]) & !is.na(d_haa$wt), ]
d_med_pw_haa$PC1A_z <- as.vector(scale(d_med_pw_haa$K5PC1A))
d_med_pw_haa$g <- as.vector(scale(d_med_pw_haa$g_raw))
d_med_pw_haa$W <- as.vector(scale(d_med_pw_haa$ave_WAIS_raw))
d_med_pw_haa$age_c <- as.vector(scale(d_med_pw_haa$age))
med_M_pw_haa <- lm(W ~ PC1A_z + sex + age_c, data = d_med_pw_haa, weights = wt)
out_Y_pw_haa <- lm(g ~ W + PC1A_z + sex + age_c, data = d_med_pw_haa, weights = wt)
set.seed(88888)
med_pw_haa <- mediate(med_M_pw_haa, out_Y_pw_haa, treat = "PC1A_z",
mediator = "W", boot = TRUE, sims = 1000)
## Running nonparametric bootstrap
cat("\nHAA: PC1A -> Parental WAIS -> Child g\n")
##
## HAA: PC1A -> Parental WAIS -> Child g
summary(med_pw_haa)
##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.0143197 -0.0284459 -0.0019528 0.028 *
## ADE -0.0933909 -0.1583070 -0.0267268 0.008 **
## Total Effect -0.1077106 -0.1757673 -0.0427747 <2e-16 ***
## Prop. Mediated 0.1329459 0.0166515 0.4002888 0.028 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 1234
##
##
## Simulations: 1000
# ==============================================================
# 21d. Non-Linearity Test: Quadratic Ancestry Term
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("ROBUSTNESS: Quadratic ancestry term\n")
## ROBUSTNESS: Quadratic ancestry term
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
run_model(g ~ PC1A_z + I(PC1A_z^2) + sire + gen_2nd + gen_3rd + sex + age_c,
d_apa, "APA: g ~ PC1A + PC1A² + controls")
##
## ================================================================================
## APA: g ~ PC1A + PC1A² + controls
## R2 = 0.0203 Adj-R2 = 0.0145
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0076 0.0579 -0.131 0.8955
## PC1A_z -0.1040 0.0566 -1.838 0.0662
## I(PC1A_z^2) +0.0236 0.0267 +0.885 0.3764
## sireBW -0.4793 0.2744 -1.747 0.0809
## sireBH +0.1061 0.1914 +0.554 0.5795
## sireBO -0.0081 0.1396 -0.058 0.9537
## gen_2ndTRUE +0.4113 0.2191 +1.877 0.0607
## gen_3rdTRUE +0.2496 0.2353 +1.061 0.2891
## sex -0.0134 0.0680 -0.197 0.8437
## age_c -0.0195 0.0307 -0.635 0.5253
## N = 1538
run_model(g ~ PC1A_z + I(PC1A_z^2) + sex + age_c,
d_haa, "HAA: g ~ PC1A + PC1A² + controls")
##
## ================================================================================
## HAA: g ~ PC1A + PC1A² + controls
## R2 = 0.009 Adj-R2 = 0.0059
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) +0.0014 0.0632 +0.023 0.9817
## PC1A_z -0.0984 0.0433 -2.272 0.0233
## I(PC1A_z^2) +0.0022 0.0158 +0.137 0.8907
## sex +0.0258 0.0722 +0.358 0.7206
## age_c +0.0057 0.0332 +0.172 0.8635
## N = 1287
# Also for PPVT (strongest individual test effect)
run_model(PPVT ~ PC1A_z + I(PC1A_z^2) + sire + gen_2nd + gen_3rd + sex + age_c,
d_apa, "APA: PPVT ~ PC1A + PC1A² + controls")
##
## ================================================================================
## APA: PPVT ~ PC1A + PC1A² + controls
## R2 = 0.0262 Adj-R2 = 0.0205
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0360 0.0559 -0.644 0.5194
## PC1A_z -0.1250 0.0551 -2.270 0.0233
## I(PC1A_z^2) +0.0335 0.0332 +1.011 0.3122
## sireBW -0.5960 0.3582 -1.664 0.0963
## sireBH +0.1288 0.2472 +0.521 0.6023
## sireBO -0.0463 0.1560 -0.297 0.7666
## gen_2ndTRUE +0.3674 0.2137 +1.720 0.0857
## gen_3rdTRUE +0.0764 0.3956 +0.193 0.8470
## sex +0.0178 0.0660 +0.269 0.7877
## age_c -0.0310 0.0316 -0.980 0.3271
## N = 1531
# ==============================================================
# 21f. Complete-case sensitivity check
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("ROBUSTNESS: Complete-case analysis (no imputation)\n")
## ROBUSTNESS: Complete-case analysis (no imputation)
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
# Identify cases with complete raw (non-imputed) data
d_apa_cc <- d_apa[!is.na(d_apa$PPVT) & !is.na(d_apa$WJ9) & !is.na(d_apa$WJ10) &
!is.na(d_apa$K5PC1A) & !is.na(d_apa$sex) & !is.na(d_apa$age) &
!is.na(d_apa$wt) & !is.na(d_apa$M_edu) & !is.na(d_apa$ave_log_inc), ]
# Build g score without imputation (PCA on complete cases)
cog_raw <- as.matrix(d_apa_cc[, c("PPVT_r", "WJ9_r", "WJ10_r")])
cog_raw_z <- apply(cog_raw, 2, scale)
pca_cc <- prcomp(cog_raw_z, center = FALSE, scale. = FALSE)
d_apa_cc$g_cc <- as.vector(pca_cc$x[, 1])
d_apa_cc$PC1A_z <- as.vector(scale(d_apa_cc$K5PC1A))
d_apa_cc$g <- as.vector(scale(d_apa_cc$g_cc))
d_apa_cc$age_c <- as.vector(scale(d_apa_cc$age))
m_cc <- lm(g ~ PC1A_z + sire + gen_2nd + gen_3rd + sex + age_c,
data = d_apa_cc, weights = wt)
robust_cc <- coeftest(m_cc, vcov = vcovHC(m_cc, type = "HC1"))
cat("APA complete-case: PC1A coefficient:\n")
## APA complete-case: PC1A coefficient:
print(robust_cc["PC1A_z", ])
## Estimate Std. Error t value Pr(>|t|)
## -0.138021014 0.048548511 -2.842950509 0.004530434
cat("N =", nrow(d_apa_cc), "\n")
## N = 1507
# ==============================================================
# 22. Descriptive Table
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("TABLE: Sample Composition\n")
## TABLE: Sample Composition
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
wt_mean <- function(x, w) {
keep <- !is.na(x) & !is.na(w) & w > 0
if (sum(keep) == 0) return(list(m = NA, n = 0))
list(m = weighted.mean(x[keep], w[keep]), n = sum(keep))
}
print_row <- function(label, subset, apa_status, haa_status) {
n <- nrow(subset)
pct <- 100 * n / nrow(d)
pc1a_m <- mean(subset$K5PC1A, na.rm = TRUE)
pc1a_sd <- sd(subset$K5PC1A, na.rm = TRUE)
g_res <- wt_mean(subset$g_raw, subset$wt)
w_res <- wt_mean(subset$ave_WAIS_raw, subset$wt)
g_str <- ifelse(is.na(g_res$m), " --", sprintf("%+.3f", g_res$m))
w_str <- ifelse(is.na(w_res$m), " --", sprintf("%.2f", w_res$m))
sd_str <- ifelse(is.na(pc1a_sd) | n < 2, " --", sprintf("%.4f", pc1a_sd))
cat(sprintf(" %-25s %5d %5.1f%% %+.4f %7s %4s %4s %5d %7s %5d %7s\n",
label, n, pct, pc1a_m, sd_str, apa_status, haa_status,
g_res$n, g_str, w_res$n, w_str))
}
cat(sprintf("\n%-27s %5s %6s %7s %7s %4s %4s %5s %7s %5s %7s\n",
"Category", "N", "%", "PC1A", "SD", "APA", "HAA", "N(g)", "Wt.g", "N(W)", "Wt.W"))
##
## Category N % PC1A SD APA HAA N(g) Wt.g N(W) Wt.W
cat(paste(rep("-", 105), collapse = ""), "\n")
## ---------------------------------------------------------------------------------------------------------
cat("Panel A: Analytic Sample\n")
## Panel A: Analytic Sample
bb_apa <- d[d$any_black_parent & !d$african_origin & d$sire == "BB", ]
print_row("Black x Black", bb_apa, "Yes", "Yes*")
## Black x Black 1418 86.5% +0.0058 0.0128 Yes Yes* 1376 -0.019 1303 0.04
bw_apa <- d[d$any_black_parent & !d$african_origin & d$sire == "BW", ]
print_row("Black x White", bw_apa, "Yes", "Excl")
## Black x White 72 4.4% -0.0671 0.0178 Yes Excl 66 +0.191 62 0.62
bh_apa <- d[d$any_black_parent & !d$african_origin & d$sire == "BH", ]
print_row("Black x Hispanic", bh_apa, "Yes", "Excl")
## Black x Hispanic 32 2.0% -0.0137 0.0304 Yes Excl 32 +0.495 29 -0.26
bo_apa <- d[d$any_black_parent & !d$african_origin & d$sire == "BO", ]
print_row("Black x Other/Missing", bo_apa, "Yes", "Excl")
## Black x Other/Missing 66 4.0% -0.0064 0.0266 Yes Excl 64 +0.177 61 -0.14
cat("\nPanel B: Continental African Origin\n")
##
## Panel B: Continental African Origin
afr <- d[d$any_black_parent & d$african_origin, ]
if (nrow(afr) > 0) print_row("African origin (excl)", afr, "Excl", "Excl")
## African origin (excl) 21 1.3% +0.0031 0.0414 Excl Excl 18 +0.106 11 0.71
cat("\nPanel C: No Black Parent\n")
##
## Panel C: No Black Parent
no_black <- d[!d$any_black_parent, ]
print_row("No Black parent", no_black, "Excl", "Excl")
## No Black parent 30 1.8% -0.0818 0.0418 Excl Excl 28 +0.547 22 0.14
cat(sprintf("\n*HAA = BB 3rd+ gen only: N=%d\n", nrow(d_haa)))
##
## *HAA = BB 3rd+ gen only: N=1323
# ==============================================================
# 23. Regression Plots
# ==============================================================
partial_resid <- function(model, data, xvar) {
r <- residuals(model)
b <- coef(model)[xvar]
r + b * data[[xvar]]
}
plot_ancestry <- function(model, data, xvar, color, ylab, main_title) {
pr <- partial_resid(model, data, xvar)
x <- data[[xvar]]
plot(x, pr, pch = 16, cex = 0.4,
col = adjustcolor("gray30", alpha.f = 0.2),
xlab = "African Ancestry PC (std.)", ylab = ylab, main = main_title)
lo <- loess(pr ~ x, span = 0.75)
x_seq <- seq(min(x), max(x), length.out = 200)
lo_pred <- predict(lo, newdata = data.frame(x = x_seq), se = TRUE)
polygon(c(x_seq, rev(x_seq)),
c(lo_pred$fit + 1.96 * lo_pred$se.fit, rev(lo_pred$fit - 1.96 * lo_pred$se.fit)),
col = adjustcolor(color, alpha.f = 0.15), border = NA)
lines(x_seq, lo_pred$fit, col = color, lwd = 2.5)
abline(lm(pr ~ x), col = color, lwd = 1.5, lty = 2)
b <- coef(model)[xvar]
p <- summary(model)$coefficients[xvar, "Pr(>|t|)"]
p_str <- ifelse(p < 0.001, "p < .001", sprintf("p = %.3f", p))
legend("topright", legend = bquote(beta == .(sprintf("%.3f", b)) ~ "," ~ .(p_str)),
bty = "n", cex = 0.9)
}
d_apa_g <- d_apa[complete.cases(d_apa[, c("g_raw", "K5PC1A", "sire", "gen_2nd", "gen_3rd", "sex", "age")]) & !is.na(d_apa$wt), ]
d_apa_g$g <- as.vector(scale(d_apa_g$g_raw)); d_apa_g$PC1A_z <- as.vector(scale(d_apa_g$K5PC1A)); d_apa_g$age_c <- as.vector(scale(d_apa_g$age))
m_apa_g <- lm(g ~ PC1A_z + sire + gen_2nd + gen_3rd + sex + age_c, data = d_apa_g, weights = d_apa_g$wt)
d_haa_g <- d_haa[complete.cases(d_haa[, c("g_raw", "K5PC1A", "sex", "age")]) & !is.na(d_haa$wt), ]
d_haa_g$g <- as.vector(scale(d_haa_g$g_raw)); d_haa_g$PC1A_z <- as.vector(scale(d_haa_g$K5PC1A)); d_haa_g$age_c <- as.vector(scale(d_haa_g$age))
m_haa_g <- lm(g ~ PC1A_z + sex + age_c, data = d_haa_g, weights = d_haa_g$wt)
d_apa_w <- d_apa_noBW[complete.cases(d_apa_noBW[, c("ave_WAIS_raw", "K5PC1A", "sire", "gen_2nd", "gen_3rd")]) & !is.na(d_apa_noBW$wt), ]
d_apa_w$ave_WAIS_res_z <- as.vector(scale(d_apa_w$ave_WAIS_raw)); d_apa_w$PC1A_z <- as.vector(scale(d_apa_w$K5PC1A))
m_apa_w <- lm(ave_WAIS_res_z ~ PC1A_z + sire + gen_2nd + gen_3rd, data = d_apa_w, weights = d_apa_w$wt)
d_haa_w <- d_haa[complete.cases(d_haa[, c("ave_WAIS_raw", "K5PC1A")]) & !is.na(d_haa$wt), ]
d_haa_w$ave_WAIS_res_z <- as.vector(scale(d_haa_w$ave_WAIS_raw)); d_haa_w$PC1A_z <- as.vector(scale(d_haa_w$K5PC1A))
m_haa_w <- lm(ave_WAIS_res_z ~ PC1A_z, data = d_haa_w, weights = d_haa_w$wt)
pdf("admixture_regression_plots.pdf", width = 10, height = 8)
par(mfrow = c(2, 2), mar = c(4.5, 4.5, 3, 1.5))
plot_ancestry(m_apa_g, d_apa_g, "PC1A_z", "blue", "Child g (z)", "A. Pan-African American: Ancestry -> Child g")
plot_ancestry(m_haa_g, d_haa_g, "PC1A_z", "blue", "Child g (z)", "B. Historical AA: Ancestry -> Child g")
plot_ancestry(m_apa_w, d_apa_w, "PC1A_z", "darkred", "Parental WAIS (z)", "C. Pan-African American: Ancestry -> Parental WAIS")
plot_ancestry(m_haa_w, d_haa_w, "PC1A_z", "darkred", "Parental WAIS (z)", "D. Historical AA: Ancestry -> Parental WAIS")
dev.off()
## png
## 2
cat("\nPlots saved to: admixture_regression_plots.pdf\n")
##
## Plots saved to: admixture_regression_plots.pdf
#_______________________________________________________________________________________
#Exploratory analysis
# ==============================================================
# Epigenetic Clock Analysis
# ==============================================================
epi <- read_sav("C:/Users/mh198/Documents/Data/FFCWS/DS0008/31622-0008-Data.sav")
# Merge with main dataset
d <- merge(d, epi[, c("IDNUM", "K5MK_POAM45", "K5MK_POAM38",
"K5MK_PHENOAGE", "K5MK_GRIM", "K5MK_HORVATH",
"K5MK_SKINBLOOD", "K5MK_PEDBE", "K5MK_HANNUM",
"K5MK_AGE", "K5MK_BATCH")],
by = "IDNUM", all.x = TRUE)
# Clean: set negative codes to NA
epi_vars <- c("K5MK_POAM45", "K5MK_POAM38", "K5MK_PHENOAGE", "K5MK_GRIM",
"K5MK_HORVATH", "K5MK_SKINBLOOD", "K5MK_PEDBE", "K5MK_HANNUM",
"K5MK_AGE", "K5MK_BATCH")
for (v in epi_vars) {
d[[v]] <- ifelse(d[[v]] < -5, NA, d[[v]])
}
cat("Epigenetic data merged. DunedinPACE available: N =", sum(!is.na(d$K5MK_POAM45)), "\n")
## Epigenetic data merged. DunedinPACE available: N = 485
# Rebuild subsets with epi data
d_apa <- d[d$any_black_parent & !d$african_origin, ]
d_apa$sire <- droplevels(d_apa$sire)
d_apa$mother_black_sire <- as.numeric(d_apa$CM1ETHRACE == 2)
d_apa$father_black_sire <- as.numeric(d_apa$CF1ETHRACE == 2)
d_apa_noBW <- d_apa[d_apa$sire != "BW", ]
d_apa_noBW$sire <- droplevels(d_apa_noBW$sire)
d_haa <- d[d$CF1ETHRACE == 2 & d$CM1ETHRACE == 2 &
!d$african_origin & !d$parent_fb & !d$gp_fb, ]
cat("\nDunedinPACE by sample:\n")
##
## DunedinPACE by sample:
cat(sprintf(" APA: N = %d\n", sum(!is.na(d_apa$K5MK_POAM45))))
## APA: N = 469
cat(sprintf(" APA excl BW: N = %d\n", sum(!is.na(d_apa_noBW$K5MK_POAM45))))
## APA excl BW: N = 457
cat(sprintf(" HAA: N = %d\n", sum(!is.na(d_haa$K5MK_POAM45))))
## HAA: N = 400
# ==============================================================
# Full Diagnostics: All Clocks
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("Epigenetic Clocks: Full Diagnostics\n")
## Epigenetic Clocks: Full Diagnostics
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
clock_vars <- c("K5MK_POAM45", "K5MK_POAM38", "K5MK_PHENOAGE", "K5MK_GRIM",
"K5MK_HORVATH", "K5MK_SKINBLOOD", "K5MK_PEDBE", "K5MK_HANNUM")
clock_labels <- c("DunedinPACE", "DunedinPoAm38", "PhenoAge", "GrimAge",
"Horvath", "SkinBlood", "PedBE", "Hannum")
# APA
cat("\nCorrelations: Clocks with Ancestry and Outcomes (APA)\n")
##
## Correlations: Clocks with Ancestry and Outcomes (APA)
cat(sprintf("%-15s %5s %9s %9s %9s %9s %9s %9s\n",
"Clock", "N", "r(PC1A)", "r(g)", "r(PPVT)", "r(WAIS)", "r(SES)", "r(WJ9)"))
## Clock N r(PC1A) r(g) r(PPVT) r(WAIS) r(SES) r(WJ9)
cat(paste(rep("-", 80), collapse = ""), "\n")
## --------------------------------------------------------------------------------
for (i in seq_along(clock_vars)) {
v <- clock_vars[i]
mask_pc <- !is.na(d_apa[[v]]) & !is.na(d_apa$K5PC1A)
mask_g <- !is.na(d_apa[[v]]) & !is.na(d_apa$g_raw)
mask_ppvt <- !is.na(d_apa[[v]]) & !is.na(d_apa$PPVT)
mask_wais <- !is.na(d_apa[[v]]) & !is.na(d_apa$ave_WAIS_raw)
mask_ses <- !is.na(d_apa[[v]]) & !is.na(d_apa$SES_raw)
mask_wj9 <- !is.na(d_apa[[v]]) & !is.na(d_apa$WJ9)
n <- sum(mask_pc)
if (n < 10) next
r_pc <- cor(d_apa[[v]][mask_pc], d_apa$K5PC1A[mask_pc])
r_g <- cor(d_apa[[v]][mask_g], d_apa$g_raw[mask_g])
r_ppvt <- cor(d_apa[[v]][mask_ppvt], d_apa$PPVT[mask_ppvt])
r_wais <- cor(d_apa[[v]][mask_wais], d_apa$ave_WAIS_raw[mask_wais])
r_ses <- cor(d_apa[[v]][mask_ses], d_apa$SES_raw[mask_ses])
r_wj9 <- cor(d_apa[[v]][mask_wj9], d_apa$WJ9[mask_wj9])
cat(sprintf("%-15s %5d %+9.4f %+9.4f %+9.4f %+9.4f %+9.4f %+9.4f\n",
clock_labels[i], n, r_pc, r_g, r_ppvt, r_wais, r_ses, r_wj9))
}
## DunedinPACE 469 +0.1137 +0.0239 +0.0240 -0.0152 -0.0032 -0.0216
## DunedinPoAm38 469 +0.1263 -0.0306 -0.0835 -0.0516 -0.0690 -0.0087
## PhenoAge 146 +0.1252 -0.0189 +0.0486 -0.0096 -0.0241 -0.0142
## GrimAge 469 +0.0509 +0.0325 +0.0539 -0.0018 +0.0066 -0.0270
## Horvath 469 +0.0023 -0.0426 -0.0148 -0.0163 +0.0216 -0.0532
## SkinBlood 469 +0.0584 +0.0090 +0.0161 +0.0119 +0.0255 +0.0338
## PedBE 469 +0.0210 +0.0136 +0.0481 -0.0053 +0.0576 +0.0214
## Hannum 469 +0.0282 +0.0551 +0.0539 -0.0226 +0.0146 -0.0015
# HAA
cat("\nCorrelations: Clocks with Ancestry and Outcomes (HAA)\n")
##
## Correlations: Clocks with Ancestry and Outcomes (HAA)
cat(sprintf("%-15s %5s %9s %9s %9s %9s %9s %9s\n",
"Clock", "N", "r(PC1A)", "r(g)", "r(PPVT)", "r(WAIS)", "r(SES)", "r(WJ9)"))
## Clock N r(PC1A) r(g) r(PPVT) r(WAIS) r(SES) r(WJ9)
cat(paste(rep("-", 80), collapse = ""), "\n")
## --------------------------------------------------------------------------------
for (i in seq_along(clock_vars)) {
v <- clock_vars[i]
mask_pc <- !is.na(d_haa[[v]]) & !is.na(d_haa$K5PC1A)
mask_g <- !is.na(d_haa[[v]]) & !is.na(d_haa$g_raw)
mask_ppvt <- !is.na(d_haa[[v]]) & !is.na(d_haa$PPVT)
mask_wais <- !is.na(d_haa[[v]]) & !is.na(d_haa$ave_WAIS_raw)
mask_ses <- !is.na(d_haa[[v]]) & !is.na(d_haa$SES_raw)
mask_wj9 <- !is.na(d_haa[[v]]) & !is.na(d_haa$WJ9)
n <- sum(mask_pc)
if (n < 10) next
r_pc <- cor(d_haa[[v]][mask_pc], d_haa$K5PC1A[mask_pc])
r_g <- cor(d_haa[[v]][mask_g], d_haa$g_raw[mask_g])
r_ppvt <- cor(d_haa[[v]][mask_ppvt], d_haa$PPVT[mask_ppvt])
r_wais <- cor(d_haa[[v]][mask_wais], d_haa$ave_WAIS_raw[mask_wais])
r_ses <- cor(d_haa[[v]][mask_ses], d_haa$SES_raw[mask_ses])
r_wj9 <- cor(d_haa[[v]][mask_wj9], d_haa$WJ9[mask_wj9])
cat(sprintf("%-15s %5d %+9.4f %+9.4f %+9.4f %+9.4f %+9.4f %+9.4f\n",
clock_labels[i], n, r_pc, r_g, r_ppvt, r_wais, r_ses, r_wj9))
}
## DunedinPACE 400 +0.0907 +0.0322 +0.0410 -0.0125 +0.0004 -0.0218
## DunedinPoAm38 400 +0.0581 -0.0004 -0.0272 -0.0518 -0.0659 +0.0002
## PhenoAge 122 +0.0485 +0.0236 +0.0695 -0.0396 -0.0024 +0.0431
## GrimAge 400 +0.0511 +0.0266 +0.0424 +0.0085 +0.0278 -0.0216
## Horvath 400 -0.0658 -0.0367 -0.0261 -0.0012 +0.0246 -0.0343
## SkinBlood 400 +0.0498 +0.0179 +0.0219 +0.0184 +0.0457 +0.0518
## PedBE 400 +0.0022 -0.0018 +0.0482 -0.0017 +0.0373 +0.0048
## Hannum 400 +0.0337 +0.0570 +0.0568 -0.0110 +0.0207 +0.0067
# ==============================================================
# Mediation: DunedinPACE
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("Mediation: Does DunedinPACE attenuate ancestry-cognition?\n")
## Mediation: Does DunedinPACE attenuate ancestry-cognition?
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
d_apa$PACE_z <- as.vector(scale(d_apa$K5MK_POAM45))
d_haa$PACE_z <- as.vector(scale(d_haa$K5MK_POAM45))
# APA: without PACE (restricted to epi sample)
run_model(g ~ PC1A_z + sire + gen_2nd + gen_3rd + sex + age_c,
d_apa[!is.na(d_apa$K5MK_POAM45), ],
"APA M1 (no PACE, epi subsample)")
##
## ================================================================================
## APA M1 (no PACE, epi subsample)
## R2 = 0.0377 Adj-R2 = 0.0209
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0272 0.1163 -0.234 0.8152
## PC1A_z -0.1603 0.0831 -1.929 0.0543
## sireBW -0.0971 0.4033 -0.241 0.8100
## sireBH +0.3796 0.2193 +1.731 0.0841
## sireBO +0.1759 0.2403 +0.732 0.4644
## gen_2ndTRUE +0.6602 0.3467 +1.904 0.0575
## gen_3rdTRUE +0.0065 0.2571 +0.025 0.9797
## sex +0.0980 0.1255 +0.781 0.4355
## age_c +0.0753 0.0527 +1.429 0.1538
## N = 467
# APA: with PACE
run_model(g ~ PC1A_z + PACE_z + sire + gen_2nd + gen_3rd + sex + age_c, d_apa,
"APA M1 + DunedinPACE")
##
## ================================================================================
## APA M1 + DunedinPACE
## R2 = 0.0382 Adj-R2 = 0.0193
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0309 0.1109 -0.279 0.7806
## PC1A_z -0.1554 0.0839 -1.852 0.0647
## PACE_z -0.0256 0.0845 -0.303 0.7618
## sireBW -0.0851 0.4008 -0.212 0.8320
## sireBH +0.3812 0.2173 +1.754 0.0801
## sireBO +0.1878 0.2329 +0.806 0.4206
## gen_2ndTRUE +0.6770 0.3501 +1.934 0.0537
## gen_3rdTRUE +0.0027 0.2662 +0.010 0.9918
## sex +0.1010 0.1220 +0.828 0.4081
## age_c +0.0757 0.0530 +1.427 0.1543
## N = 467
# HAA: without PACE (restricted to epi sample)
run_model(g ~ PC1A_z + sex + age_c,
d_haa[!is.na(d_haa$K5MK_POAM45), ],
"HAA M1 (no PACE, epi subsample)")
##
## ================================================================================
## HAA M1 (no PACE, epi subsample)
## R2 = 0.0172 Adj-R2 = 0.0097
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0250 0.1190 -0.210 0.8335
## PC1A_z -0.1089 0.0574 -1.896 0.0586
## sex +0.1471 0.1348 +1.091 0.2758
## age_c +0.0907 0.0566 +1.601 0.1102
## N = 398
# HAA: with PACE
run_model(g ~ PC1A_z + PACE_z + sex + age_c, d_haa,
"HAA M1 + DunedinPACE")
##
## ================================================================================
## HAA M1 + DunedinPACE
## R2 = 0.0189 Adj-R2 = 0.0089
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0319 0.1119 -0.285 0.7757
## PC1A_z -0.1048 0.0580 -1.808 0.0714
## PACE_z -0.0461 0.0944 -0.488 0.6255
## sex +0.1563 0.1274 +1.227 0.2207
## age_c +0.0910 0.0568 +1.603 0.1098
## N = 398
# ==============================================================
# Mediation: GrimAge
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("Mediation: Does GrimAge attenuate ancestry-cognition?\n")
## Mediation: Does GrimAge attenuate ancestry-cognition?
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
d_apa$GRIM_z <- as.vector(scale(d_apa$K5MK_GRIM))
d_haa$GRIM_z <- as.vector(scale(d_haa$K5MK_GRIM))
run_model(g ~ PC1A_z + GRIM_z + sire + gen_2nd + gen_3rd + sex + age_c, d_apa,
"APA M1 + GrimAge")
##
## ================================================================================
## APA M1 + GrimAge
## R2 = 0.0377 Adj-R2 = 0.0188
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0273 0.1178 -0.231 0.8171
## PC1A_z -0.1605 0.0834 -1.925 0.0548
## GRIM_z +0.0012 0.0674 +0.018 0.9854
## sireBW -0.0977 0.4027 -0.243 0.8084
## sireBH +0.3796 0.2197 +1.728 0.0847
## sireBO +0.1754 0.2381 +0.737 0.4617
## gen_2ndTRUE +0.6601 0.3467 +1.904 0.0575
## gen_3rdTRUE +0.0063 0.2577 +0.024 0.9806
## sex +0.0982 0.1310 +0.750 0.4538
## age_c +0.0752 0.0531 +1.416 0.1576
## N = 467
run_model(g ~ PC1A_z + GRIM_z + sex + age_c, d_haa,
"HAA M1 + GrimAge")
##
## ================================================================================
## HAA M1 + GrimAge
## R2 = 0.0172 Adj-R2 = 0.0072
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0251 0.1214 -0.207 0.8364
## PC1A_z -0.1090 0.0576 -1.892 0.0592
## GRIM_z +0.0008 0.0739 +0.010 0.9919
## sex +0.1473 0.1397 +1.054 0.2925
## age_c +0.0906 0.0571 +1.587 0.1134
## N = 398
# ==============================================================
# Mediation: All clocks (loop)
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("Mediation Summary: All Clocks (HAA)\n")
## Mediation Summary: All Clocks (HAA)
cat("PC1A_z coefficient with and without each clock\n")
## PC1A_z coefficient with and without each clock
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
cat(sprintf("\n%-15s %5s %9s %9s %9s\n",
"Clock", "N", "PC1A(no)", "PC1A(+clk)", "Change"))
##
## Clock N PC1A(no) PC1A(+clk) Change
cat(paste(rep("-", 50), collapse = ""), "\n")
## --------------------------------------------------
for (i in seq_along(clock_vars)) {
v <- clock_vars[i]
dd <- d_haa[!is.na(d_haa[[v]]) & !is.na(d_haa$g_raw) &
!is.na(d_haa$K5PC1A) & !is.na(d_haa$age) &
!is.na(d_haa$sex) & !is.na(d_haa$wt), ]
if (nrow(dd) < 30) next
dd$PC1A_z <- as.vector(scale(dd$K5PC1A))
dd$g <- as.vector(scale(dd$g_raw))
dd$age_c <- as.vector(scale(dd$age))
dd$clock_z <- as.vector(scale(dd[[v]]))
m0 <- lm(g ~ PC1A_z + sex + age_c, data = dd, weights = wt)
m1 <- lm(g ~ PC1A_z + clock_z + sex + age_c, data = dd, weights = wt)
b0 <- coef(m0)["PC1A_z"]
b1 <- coef(m1)["PC1A_z"]
change_pct <- 100 * (b1 - b0) / b0
cat(sprintf("%-15s %5d %+9.4f %+9.4f %+8.1f%%\n",
clock_labels[i], nrow(dd), b0, b1, change_pct))
}
## DunedinPACE 398 -0.1089 -0.1048 -3.8%
## DunedinPoAm38 398 -0.1089 -0.1077 -1.1%
## PhenoAge 121 -0.1080 -0.1085 +0.5%
## GrimAge 398 -0.1089 -0.1090 +0.0%
## Horvath 398 -0.1089 -0.1150 +5.6%
## SkinBlood 398 -0.1089 -0.1087 -0.2%
## PedBE 398 -0.1089 -0.1094 +0.4%
## Hannum 398 -0.1089 -0.1090 +0.0%
# ==============================================================
# Epigenetic Formal Mediation: PC1A -> DunedinPACE -> g
# ==============================================================
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("FORMAL MEDIATION: PC1A -> DunedinPACE -> g (HAA)\n")
## FORMAL MEDIATION: PC1A -> DunedinPACE -> g (HAA)
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
d_epi_haa <- d_haa[!is.na(d_haa$K5MK_POAM45) & !is.na(d_haa$g_raw) &
!is.na(d_haa$K5PC1A) & !is.na(d_haa$age) &
!is.na(d_haa$sex) & !is.na(d_haa$wt), ]
d_epi_haa$PC1A_z <- as.vector(scale(d_epi_haa$K5PC1A))
d_epi_haa$g <- as.vector(scale(d_epi_haa$g_raw))
d_epi_haa$PACE_z <- as.vector(scale(d_epi_haa$K5MK_POAM45))
d_epi_haa$age_c <- as.vector(scale(d_epi_haa$age))
med_M_epi <- lm(PACE_z ~ PC1A_z + sex + age_c, data = d_epi_haa, weights = wt)
out_Y_epi <- lm(g ~ PACE_z + PC1A_z + sex + age_c, data = d_epi_haa, weights = wt)
set.seed(77777)
med_epi <- mediate(med_M_epi, out_Y_epi, treat = "PC1A_z", mediator = "PACE_z",
boot = TRUE, sims = 1000)
## Running nonparametric bootstrap
summary(med_epi)
##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.0041069 -0.0300480 0.0137365 0.680
## ADE -0.1048220 -0.2179324 0.0066505 0.068 .
## Total Effect -0.1089289 -0.2197694 0.0050707 0.060 .
## Prop. Mediated 0.0377025 -0.2252410 0.5003497 0.684
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 398
##
##
## Simulations: 1000
################
#"Epigenetic pace-of-aging measures (DunedinPACE, DunedinPoAm38, GrimAge, Horvath, and others) did not attenuate the ancestry-cognition association (maximum change: 3.6%), and showed near-zero correlations with cognitive outcomes (|r| < .05 for all clocks).
#These results do not support epigenetic aging as a mediator of ancestry-related cognitive differences in this sample."
# ==============================================================
# SM-A. Supplementary Material: Moderation Analyses
# ==============================================================
# Note: These use individual subtests (PPVT, WJ9, WJ10) as outcomes to avoid measurement invariance assumptions.
cat("\n", paste(rep("#", 80), collapse = ""), "\n")
##
## ################################################################################
cat("SUPPLEMENTARY MATERIAL: Moderation Analyses\n")
## SUPPLEMENTARY MATERIAL: Moderation Analyses
cat("(Using individual subtests to avoid MI assumptions)\n")
## (Using individual subtests to avoid MI assumptions)
cat(paste(rep("#", 80), collapse = ""), "\n")
## ################################################################################
# SM-A1: Sex × Ancestry on individual tests
cat("\n--- SM-A1: Sex × Ancestry (APA) ---\n")
##
## --- SM-A1: Sex × Ancestry (APA) ---
for (test in c("PPVT", "WJ9", "WJ10")) {
run_model(
as.formula(paste(test, "~ PC1A_z * sex + sire + gen_2nd + gen_3rd + age_c")),
d_apa,
paste("APA:", test, "~ PC1A × Sex")
)
}
##
## ================================================================================
## APA: PPVT ~ PC1A × Sex
## R2 = 0.0309 Adj-R2 = 0.0252
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0198 0.0524 -0.379 0.7051
## PC1A_z -0.2424 0.0594 -4.083 4.68e-05
## sex +0.0161 0.0659 +0.244 0.8075
## sireBW -0.3294 0.2261 -1.457 0.1453
## sireBH +0.1358 0.2210 +0.615 0.5389
## sireBO -0.0315 0.1551 -0.203 0.8389
## gen_2ndTRUE +0.4156 0.2080 +1.998 0.0459
## gen_3rdTRUE +0.0574 0.3918 +0.147 0.8835
## age_c -0.0292 0.0312 -0.935 0.3501
## PC1A_z:sex +0.1706 0.0648 +2.631 0.0086
## N = 1531
##
## ================================================================================
## APA: WJ9 ~ PC1A × Sex
## R2 = 0.0316 Adj-R2 = 0.0258
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.1513 0.0539 -2.808 0.0050
## PC1A_z -0.0884 0.0567 -1.559 0.1192
## sex +0.3026 0.0640 +4.730 2.45e-06
## sireBW -0.2953 0.2105 -1.402 0.1610
## sireBH -0.0111 0.1546 -0.072 0.9429
## sireBO +0.0993 0.1559 +0.637 0.5242
## gen_2ndTRUE +0.2991 0.1850 +1.617 0.1061
## gen_3rdTRUE +0.1102 0.1838 +0.600 0.5486
## age_c -0.0496 0.0306 -1.620 0.1054
## PC1A_z:sex +0.0198 0.0555 +0.357 0.7208
## N = 1523
##
## ================================================================================
## APA: WJ10 ~ PC1A × Sex
## R2 = 0.024 Adj-R2 = 0.0182
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0388 0.0569 -0.683 0.4948
## PC1A_z -0.1437 0.0549 -2.619 0.0089
## sex +0.1079 0.0682 +1.580 0.1142
## sireBW -0.1853 0.1964 -0.943 0.3456
## sireBH +0.0660 0.1500 +0.440 0.6600
## sireBO -0.0518 0.1332 -0.389 0.6975
## gen_2ndTRUE +0.3916 0.2110 +1.856 0.0637
## gen_3rdTRUE +0.4688 0.1826 +2.567 0.0103
## age_c -0.0884 0.0292 -3.024 0.0025
## PC1A_z:sex +0.0901 0.0561 +1.606 0.1084
## N = 1524
cat("\n--- SM-A1: Sex × Ancestry (HAA) ---\n")
##
## --- SM-A1: Sex × Ancestry (HAA) ---
for (test in c("PPVT", "WJ9", "WJ10")) {
run_model(
as.formula(paste(test, "~ PC1A_z * sex + age_c")),
d_haa,
paste("HAA:", test, "~ PC1A × Sex")
)
}
##
## ================================================================================
## HAA: PPVT ~ PC1A × Sex
## R2 = 0.0191 Adj-R2 = 0.016
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0483 0.0559 -0.865 0.3872
## PC1A_z -0.1780 0.0443 -4.018 6.22e-05
## sex +0.1032 0.0719 +1.435 0.1515
## age_c -0.0161 0.0344 -0.466 0.6410
## PC1A_z:sex +0.1206 0.0623 +1.938 0.0529
## N = 1281
##
## ================================================================================
## HAA: WJ9 ~ PC1A × Sex
## R2 = 0.0234 Adj-R2 = 0.0203
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.1429 0.0562 -2.541 0.0112
## PC1A_z -0.0484 0.0494 -0.981 0.3269
## sex +0.3018 0.0689 +4.382 1.27e-05
## age_c -0.0189 0.0333 -0.569 0.5693
## PC1A_z:sex -0.0207 0.0623 -0.333 0.7394
## N = 1273
##
## ================================================================================
## HAA: WJ10 ~ PC1A × Sex
## R2 = 0.0161 Adj-R2 = 0.013
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0342 0.0603 -0.568 0.5699
## PC1A_z -0.1331 0.0457 -2.916 0.0036
## sex +0.1195 0.0736 +1.623 0.1047
## age_c -0.0790 0.0315 -2.506 0.0123
## PC1A_z:sex +0.1198 0.0601 +1.992 0.0466
## N = 1277
# SM-A2: Age × Ancestry on individual tests
# Note: Age is already controlled for linearly in the main models.
# Here we test whether the ancestry effect varies with age.
# And we use the raw (non-centered) age to make the interaction interpretable.
cat("\n--- SM-A2: Age × Ancestry (APA) ---\n")
##
## --- SM-A2: Age × Ancestry (APA) ---
for (test in c("PPVT", "WJ9", "WJ10")) {
run_model(
as.formula(paste(test, "~ PC1A_z * age_c + sire + gen_2nd + gen_3rd + sex")),
d_apa,
paste("APA:", test, "~ PC1A × Age")
)
}
##
## ================================================================================
## APA: PPVT ~ PC1A × Age
## R2 = 0.0246 Adj-R2 = 0.0189
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0159 0.0526 -0.301 0.7631
## PC1A_z -0.1636 0.0510 -3.206 0.0014
## age_c -0.0286 0.0318 -0.900 0.3683
## sireBW -0.3650 0.2272 -1.606 0.1084
## sireBH +0.1366 0.2376 +0.575 0.5653
## sireBO -0.0259 0.1606 -0.161 0.8719
## gen_2ndTRUE +0.3846 0.2133 +1.803 0.0715
## gen_3rdTRUE +0.0633 0.3953 +0.160 0.8729
## sex +0.0211 0.0661 +0.319 0.7495
## PC1A_z:age_c +0.0175 0.0309 +0.566 0.5715
## N = 1531
##
## ================================================================================
## APA: WJ9 ~ PC1A × Age
## R2 = 0.032 Adj-R2 = 0.0263
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.1500 0.0539 -2.781 0.0055
## PC1A_z -0.0815 0.0453 -1.799 0.0723
## age_c -0.0494 0.0306 -1.614 0.1067
## sireBW -0.3210 0.2091 -1.535 0.1250
## sireBH -0.0155 0.1554 -0.100 0.9207
## sireBO +0.1005 0.1556 +0.646 0.5184
## gen_2ndTRUE +0.2977 0.1854 +1.605 0.1086
## gen_3rdTRUE +0.1093 0.1847 +0.592 0.5540
## sex +0.3032 0.0641 +4.730 2.45e-06
## PC1A_z:age_c +0.0237 0.0240 +0.987 0.3236
## N = 1523
##
## ================================================================================
## APA: WJ10 ~ PC1A × Age
## R2 = 0.0224 Adj-R2 = 0.0166
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0379 0.0569 -0.666 0.5058
## PC1A_z -0.0999 0.0444 -2.249 0.0247
## age_c -0.0881 0.0292 -3.014 0.0026
## sireBW -0.1829 0.2003 -0.913 0.3615
## sireBH +0.0713 0.1505 +0.474 0.6358
## sireBO -0.0495 0.1348 -0.367 0.7134
## gen_2ndTRUE +0.3727 0.2102 +1.773 0.0765
## gen_3rdTRUE +0.4732 0.1872 +2.528 0.0116
## sex +0.1105 0.0686 +1.612 0.1072
## PC1A_z:age_c -0.0153 0.0228 -0.673 0.5013
## N = 1524
cat("\n--- SM-A2: Age × Ancestry (HAA) ---\n")
##
## --- SM-A2: Age × Ancestry (HAA) ---
for (test in c("PPVT", "WJ9", "WJ10")) {
run_model(
as.formula(paste(test, "~ PC1A_z * age_c + sex")),
d_haa,
paste("HAA:", test, "~ PC1A × Age")
)
}
##
## ================================================================================
## HAA: PPVT ~ PC1A × Age
## R2 = 0.0158 Adj-R2 = 0.0127
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0442 0.0559 -0.791 0.4290
## PC1A_z -0.1227 0.0316 -3.887 0.0001
## age_c -0.0138 0.0347 -0.397 0.6914
## sex +0.1079 0.0719 +1.502 0.1334
## PC1A_z:age_c -0.0024 0.0266 -0.090 0.9281
## N = 1281
##
## ================================================================================
## HAA: WJ9 ~ PC1A × Age
## R2 = 0.0236 Adj-R2 = 0.0205
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.1437 0.0565 -2.545 0.0111
## PC1A_z -0.0566 0.0323 -1.754 0.0796
## age_c -0.0189 0.0334 -0.565 0.5720
## sex +0.3013 0.0689 +4.373 1.33e-05
## PC1A_z:age_c -0.0162 0.0264 -0.611 0.5410
## N = 1273
##
## ================================================================================
## HAA: WJ10 ~ PC1A × Age
## R2 = 0.0139 Adj-R2 = 0.0108
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0306 0.0606 -0.506 0.6130
## PC1A_z -0.0766 0.0312 -2.450 0.0144
## age_c -0.0761 0.0315 -2.416 0.0158
## sex +0.1249 0.0741 +1.686 0.0921
## PC1A_z:age_c -0.0298 0.0280 -1.065 0.2873
## N = 1277
# SM-A3: Parental SES × Ancestry on individual tests
# SES is continuous (standardized). Test whether ancestry effects differ across the SES distribution.
cat("\n--- SM-A3: SES × Ancestry (APA, excl BW) ---\n")
##
## --- SM-A3: SES × Ancestry (APA, excl BW) ---
for (test in c("PPVT", "WJ9", "WJ10")) {
run_model(
as.formula(paste(test, "~ PC1A_z * SES + sire + gen_2nd + gen_3rd + sex + age_c")),
d_apa_noBW,
paste("APA (excl BW):", test, "~ PC1A × SES")
)
}
##
## ================================================================================
## APA (excl BW): PPVT ~ PC1A × SES
## R2 = 0.1289 Adj-R2 = 0.1235
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0496 0.0502 -0.989 0.3230
## PC1A_z -0.0971 0.0334 -2.904 0.0037
## SES +0.3189 0.0309 +10.326 3.59e-24
## sireBH +0.3249 0.2372 +1.369 0.1711
## sireBO +0.0253 0.1484 +0.170 0.8648
## gen_2ndTRUE +0.2068 0.1939 +1.067 0.2862
## gen_3rdTRUE +0.0495 0.4177 +0.118 0.9058
## sex +0.0392 0.0650 +0.602 0.5472
## age_c -0.0074 0.0309 -0.238 0.8119
## PC1A_z:SES +0.0516 0.0344 +1.499 0.1340
## N = 1466
##
## ================================================================================
## APA (excl BW): WJ9 ~ PC1A × SES
## R2 = 0.0921 Adj-R2 = 0.0864
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.1829 0.0514 -3.558 0.0004
## PC1A_z -0.0317 0.0341 -0.929 0.3530
## SES +0.2502 0.0300 +8.336 1.77e-16
## sireBH +0.1219 0.1585 +0.769 0.4422
## sireBO +0.1163 0.1499 +0.775 0.4382
## gen_2ndTRUE +0.1671 0.1663 +1.005 0.3150
## gen_3rdTRUE +0.0974 0.1980 +0.492 0.6230
## sex +0.3143 0.0634 +4.956 8.05e-07
## age_c -0.0295 0.0303 -0.974 0.3303
## PC1A_z:SES -0.0153 0.0260 -0.587 0.5570
## N = 1457
##
## ================================================================================
## APA (excl BW): WJ10 ~ PC1A × SES
## R2 = 0.0991 Adj-R2 = 0.0935
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0695 0.0528 -1.316 0.1883
## PC1A_z -0.0498 0.0317 -1.569 0.1169
## SES +0.2823 0.0315 +8.958 9.95e-19
## sireBH +0.2222 0.1427 +1.558 0.1195
## sireBO -0.0205 0.1278 -0.160 0.8727
## gen_2ndTRUE +0.2258 0.1887 +1.197 0.2316
## gen_3rdTRUE +0.4599 0.2009 +2.289 0.0222
## sex +0.1280 0.0658 +1.946 0.0518
## age_c -0.0723 0.0285 -2.541 0.0112
## PC1A_z:SES +0.0068 0.0248 +0.274 0.7841
## N = 1458
cat("\n--- SM-A3: SES × Ancestry (HAA) ---\n")
##
## --- SM-A3: SES × Ancestry (HAA) ---
for (test in c("PPVT", "WJ9", "WJ10")) {
run_model(
as.formula(paste(test, "~ PC1A_z * SES + sex + age_c")),
d_haa,
paste("HAA:", test, "~ PC1A × SES")
)
}
##
## ================================================================================
## HAA: PPVT ~ PC1A × SES
## R2 = 0.1188 Adj-R2 = 0.1153
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0740 0.0530 -1.396 0.1630
## PC1A_z -0.0983 0.0298 -3.298 0.0010
## SES +0.3161 0.0306 +10.319 4.93e-24
## sex +0.1225 0.0691 +1.771 0.0767
## age_c +0.0053 0.0333 +0.161 0.8725
## PC1A_z:SES +0.0015 0.0289 +0.053 0.9575
## N = 1281
##
## ================================================================================
## HAA: WJ9 ~ PC1A × SES
## R2 = 0.0795 Adj-R2 = 0.0758
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.1659 0.0547 -3.033 0.0025
## PC1A_z -0.0366 0.0329 -1.110 0.2673
## SES +0.2373 0.0307 +7.731 2.17e-14
## sex +0.3107 0.0672 +4.620 4.23e-06
## age_c -0.0042 0.0330 -0.127 0.8989
## PC1A_z:SES -0.0131 0.0258 -0.507 0.6124
## N = 1273
##
## ================================================================================
## HAA: WJ10 ~ PC1A × SES
## R2 = 0.0893 Adj-R2 = 0.0857
## ================================================================================
## Variable B SE t p
## ------------------------------------------------------------
## (Intercept) -0.0549 0.0558 -0.983 0.3259
## PC1A_z -0.0573 0.0302 -1.902 0.0574
## SES +0.2781 0.0325 +8.551 3.46e-17
## sex +0.1350 0.0701 +1.927 0.0542
## age_c -0.0591 0.0307 -1.923 0.0547
## PC1A_z:SES +0.0020 0.0242 +0.083 0.9337
## N = 1277