Indian IQ in the GSS (Wordsum)
01 September, 2025
Init
options(digits = 2)
#remotes::install_github("tidyverse/haven")
library(kirkegaard)## Loading required package: tidyverse## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: magrittr
##
##
## Attaching package: 'magrittr'
##
##
## The following object is masked from 'package:purrr':
##
## set_names
##
##
## The following object is masked from 'package:tidyr':
##
## extract
##
##
## Loading required package: weights
##
## Loading required package: Hmisc
##
##
## Attaching package: 'Hmisc'
##
##
## The following objects are masked from 'package:dplyr':
##
## src, summarize
##
##
## The following objects are masked from 'package:base':
##
## format.pval, units
##
##
## Loading required package: assertthat
##
##
## Attaching package: 'assertthat'
##
##
## The following object is masked from 'package:tibble':
##
## has_name
##
##
## Loading required package: psych
##
##
## Attaching package: 'psych'
##
##
## The following object is masked from 'package:Hmisc':
##
## describe
##
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
##
##
##
## Attaching package: 'kirkegaard'
##
##
## The following object is masked from 'package:psych':
##
## rescale
##
##
## The following object is masked from 'package:assertthat':
##
## are_equal
##
##
## The following object is masked from 'package:purrr':
##
## is_logical
##
##
## The following object is masked from 'package:base':
##
## +load_packages(
haven,
rms,
mirt,
ggeffects
)## Loading required package: stats4
## Loading required package: latticetheme_set(theme_bw())Functions
Data
d = read_stata("data/GSS_stata/gss7224_r1.dta")
d_vars = df_var_table(d)Recode
#race
d$race %>% as_factor() %>% table2(include_NA = F)d$race_3way = d$race %>% as_factor()
#ethnic origin
d$ethnic2 = d$ethnic %>% as_factor()
d$ethnic2 %>% table2() %>% print(n=20)## # A tibble: 115 × 3
## Group Count Percent
## <chr> <dbl> <dbl>
## 1 uncodeable 13257 17.5
## 2 germany 9936 13.1
## 3 england & wales 7883 10.4
## 4 ireland 7147 9.44
## 5 africa (general) 4779 6.31
## 6 italy 3407 4.50
## 7 mexico 3011 3.98
## 8 american indian 2764 3.65
## 9 scotland 2021 2.67
## 10 poland 1763 2.33
## 11 american only 1460 1.93
## 12 france 1234 1.63
## 13 norway 1048 1.38
## 14 sweden 951 1.26
## 15 netherlands 887 1.17
## 16 russia 840 1.11
## 17 other spanish 729 0.963
## 18 czechoslovakia 723 0.955
## 19 don't know 722 0.954
## 20 puerto rico 716 0.946
## # ℹ 95 more rows#sex
d %<>% mutate(
sex = sex %>% as_factor(),
sex_2way = case_when(
sex %in% c("female", "male") ~ sex,
.default = NA
) %>% fct_drop()
)
#year and age
d$year %<>% as.numeric()
d$age %<>% as.numeric()Analysis
Score WORDSUM
#wordsum items
wordsum_items = d %>% select(worda:wordj)
#IRT fit
wordsum_irt = mirt(
wordsum_items %>% map_df(as.numeric)
)## Warning: data contains response patterns with only NAs## Iteration: 1, Log-Lik: -151624.585, Max-Change: 0.91761Iteration: 2, Log-Lik: -147684.585, Max-Change: 0.47336Iteration: 3, Log-Lik: -146652.257, Max-Change: 0.20919Iteration: 4, Log-Lik: -146509.074, Max-Change: 0.11647Iteration: 5, Log-Lik: -146472.989, Max-Change: 0.06375Iteration: 6, Log-Lik: -146462.250, Max-Change: 0.04367Iteration: 7, Log-Lik: -146458.415, Max-Change: 0.03077Iteration: 8, Log-Lik: -146456.888, Max-Change: 0.01850Iteration: 9, Log-Lik: -146456.321, Max-Change: 0.00929Iteration: 10, Log-Lik: -146456.156, Max-Change: 0.00533Iteration: 11, Log-Lik: -146456.090, Max-Change: 0.00244Iteration: 12, Log-Lik: -146456.064, Max-Change: 0.00209Iteration: 13, Log-Lik: -146456.042, Max-Change: 0.00068Iteration: 14, Log-Lik: -146456.036, Max-Change: 0.00068Iteration: 15, Log-Lik: -146456.034, Max-Change: 0.00051Iteration: 16, Log-Lik: -146456.033, Max-Change: 0.00037Iteration: 17, Log-Lik: -146456.032, Max-Change: 0.00029Iteration: 18, Log-Lik: -146456.031, Max-Change: 0.00028Iteration: 19, Log-Lik: -146456.031, Max-Change: 0.00026Iteration: 20, Log-Lik: -146456.030, Max-Change: 0.00024Iteration: 21, Log-Lik: -146456.030, Max-Change: 0.00025Iteration: 22, Log-Lik: -146456.029, Max-Change: 0.00023Iteration: 23, Log-Lik: -146456.029, Max-Change: 0.00023Iteration: 24, Log-Lik: -146456.029, Max-Change: 0.00021Iteration: 25, Log-Lik: -146456.029, Max-Change: 0.00021Iteration: 26, Log-Lik: -146456.028, Max-Change: 0.00019Iteration: 27, Log-Lik: -146456.028, Max-Change: 0.00019Iteration: 28, Log-Lik: -146456.028, Max-Change: 0.00018Iteration: 29, Log-Lik: -146456.028, Max-Change: 0.00018Iteration: 30, Log-Lik: -146456.027, Max-Change: 0.00016Iteration: 31, Log-Lik: -146456.027, Max-Change: 0.00016Iteration: 32, Log-Lik: -146456.027, Max-Change: 0.00015Iteration: 33, Log-Lik: -146456.027, Max-Change: 0.00015Iteration: 34, Log-Lik: -146456.027, Max-Change: 0.00014Iteration: 35, Log-Lik: -146456.027, Max-Change: 0.00014Iteration: 36, Log-Lik: -146456.027, Max-Change: 0.00012Iteration: 37, Log-Lik: -146456.027, Max-Change: 0.00012Iteration: 38, Log-Lik: -146456.026, Max-Change: 0.00011Iteration: 39, Log-Lik: -146456.026, Max-Change: 0.00011Iteration: 40, Log-Lik: -146456.026, Max-Change: 0.00011Iteration: 41, Log-Lik: -146456.026, Max-Change: 0.00010Iteration: 42, Log-Lik: -146456.026, Max-Change: 0.00010#scores
wordsum_irt_scores = fscores(wordsum_irt)
#save the white standard scores
d$wordsum_g = wordsum_irt_scores[, 1] %>% standardize()Indian IQ
#catch indians in Ethnic variable
#but we can't use this for whites, too many countries
d$Indian = (d$ethnic2 == "india")
d$Indian %>% sum()## [1] 406#examine the missing data structure
d %>% select(
Indian, race_3way, age, born, wordsum_g, sex_2way, year
) %>%
miss_patterns()#then subset to variables with no missing
d_indian = d %>% select(
Indian, race_3way, age, born, wordsum_g, year, sex_2way
) %>%
miss_filter() %>%
filter(
Indian | race_3way == "white"
)
#we need special encoding for race here
d_indian$Indian2 = tibble(
Indian = d_indian$Indian,
White = d_indian$race_3way=="white",
US_born = d_indian$born==1
) %>%
encode_combinations() %>%
fct_relevel("White, US_born")
#counts
d_indian$Indian2 %>% table2()#age std
wordsum_age_adj = make_norms(
d_indian$wordsum_g,
age = d_indian$age,
norm_group = d_indian$Indian2=="White, US_born"
)## Detected linear effect of age on the score (p = <0.001***). Model used.## Detected variance effect of age on the score (p = <0.001***). Model used.d_indian$wordsum_IQ = wordsum_age_adj$data$IQ
#desc by group
describe2(
d_indian$wordsum_IQ,
d_indian$Indian2
)## New names:
## • `` -> `...1`#ols
ols_list = list(
ols(wordsum_IQ ~ Indian2, data = d_indian),
ols(wordsum_IQ ~ Indian2 + sex_2way, data = d_indian),
ols(wordsum_IQ ~ Indian2 + year + sex_2way, data = d_indian)
)
ols_list %>% summarize_models()ols_list[[3]] %>%
ggaverage(terms = c("Indian2")) %>%
plot()
Meta
#versions
write_sessioninfo()## R version 4.5.1 (2025-06-13)
## Platform: x86_64-pc-linux-gnu
## Running under: Linux Mint 21.1
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0 LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_DK.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_DK.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_DK.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Brussels
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] ggeffects_2.2.1 mirt_1.44.0 lattice_0.22-5
## [4] rms_8.0-0 haven_2.5.4 kirkegaard_2025-09-01
## [7] psych_2.5.3 assertthat_0.2.1 weights_1.0.4
## [10] Hmisc_5.2-3 magrittr_2.0.3 lubridate_1.9.4
## [13] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
## [16] purrr_1.0.4 readr_2.1.5 tidyr_1.3.1
## [19] tibble_3.2.1 ggplot2_3.5.2 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] rstudioapi_0.17.1 audio_0.1-11 jsonlite_2.0.0
## [4] shape_1.4.6.1 datawizard_1.2.0 TH.data_1.1-3
## [7] jomo_2.7-6 farver_2.1.2 nloptr_2.2.1
## [10] rmarkdown_2.29 vctrs_0.6.5 minqa_1.2.8
## [13] base64enc_0.1-3 htmltools_0.5.8.1 polspline_1.1.25
## [16] broom_1.0.8 Formula_1.2-5 mitml_0.4-5
## [19] dcurver_0.9.2 sass_0.4.10 parallelly_1.43.0
## [22] bslib_0.9.0 htmlwidgets_1.6.4 plyr_1.8.9
## [25] sandwich_3.1-1 testthat_3.2.3 zoo_1.8-14
## [28] cachem_1.1.0 lifecycle_1.0.4 iterators_1.0.14
## [31] pkgconfig_2.0.3 Matrix_1.7-3 R6_2.6.1
## [34] fastmap_1.2.0 rbibutils_2.3 future_1.40.0
## [37] digest_0.6.37 colorspace_2.1-1 vegan_2.6-10
## [40] labeling_0.4.3 progressr_0.15.1 timechange_0.3.0
## [43] gdata_3.0.1 mgcv_1.9-1 compiler_4.5.1
## [46] withr_3.0.2 htmlTable_2.4.3 backports_1.5.0
## [49] R.utils_2.13.0 pan_1.9 MASS_7.3-65
## [52] quantreg_6.1 sessioninfo_1.2.3 GPArotation_2025.3-1
## [55] gtools_3.9.5 permute_0.9-7 tools_4.5.1
## [58] foreign_0.8-90 future.apply_1.11.3 nnet_7.3-20
## [61] R.oo_1.27.0 glue_1.8.0 nlme_3.1-168
## [64] grid_4.5.1 checkmate_2.3.2 cluster_2.1.8.1
## [67] generics_0.1.3 gtable_0.3.6 tzdb_0.5.0
## [70] R.methodsS3_1.8.2 data.table_1.17.0 hms_1.1.3
## [73] Deriv_4.1.6 utf8_1.2.4 foreach_1.5.2
## [76] pillar_1.10.2 splines_4.5.1 survival_3.8-3
## [79] SparseM_1.84-2 tidyselect_1.2.1 pbapply_1.7-2
## [82] knitr_1.50 reformulas_0.4.0 gridExtra_2.3
## [85] xfun_0.52 brio_1.1.5 stringi_1.8.7
## [88] yaml_2.3.10 boot_1.3-31 evaluate_1.0.3
## [91] codetools_0.2-19 beepr_2.0 cli_3.6.4
## [94] rpart_4.1.24 Rdpack_2.6.4 munsell_0.5.1
## [97] jquerylib_0.1.4 Rcpp_1.0.14 globals_0.17.0
## [100] parallel_4.5.1 MatrixModels_0.5-4 marginaleffects_0.28.0
## [103] lme4_1.1-37 listenv_0.9.1 glmnet_4.1-8
## [106] mvtnorm_1.3-3 SimDesign_2.19.2 scales_1.3.0
## [109] insight_1.3.1 rlang_1.1.6 multcomp_1.4-28
## [112] mnormt_2.1.1 mice_3.17.0#write data to file for reuse
# d %>% write_rds("data/data_for_reuse.rds")
#OSF
if (F) {
library(osfr)
#login
osf_auth(readr::read_lines("~/.config/osf_token"))
#the project we will use
osf_proj = osf_retrieve_node("https://osf.io/XXX/")
#upload all files in project
#overwrite existing (versioning)
osf_upload(
osf_proj,
path = c("data", "figures", "papers", "notebook.Rmd", "notebook.html", "sessions_info.txt"),
conflicts = "overwrite"
)
}