Theoretical motivation

Research Questions

Are item difficulties consistent in English and Vietnamese? Is this mediated by item-specific characteristics, e.g., noun/verb, concreteness, socialness, phonetics etc. How consistent are item difficulty ranking and age of acquisition norms? AoA norms are difficult to estimate for bilingual children, because they interact with age of acquisition onset (Schulz & Grimm 2019). So it’s possible that item difficulty ranking is a better measure for trajectory. Are there demographic factors that determine item difficulty (see Kachergis paper) Possible big question: How much of bilingual word learning trajectory is determined by difficulty of learning underlying concepts vs. language-specific inputs?

Overview

Three questions about item difficulty, in increasing order of how much cross-language machinery they need:

  1. Part 1 – Item difficulty trajectories. Within each language, which kinds of words (by CDI category, item_kind) tend to be produced earliest/most easily? Does the same category-level pattern show up in both English and Vietnamese? This only needs item_kind, which is already shared across the two fields files – no item-level translation matching required.
  2. Part 2 – Demographic bias (DIF). Following Kachergis, Francis, & Frank (2022, Estimating demographic bias on tests of children’s early vocabulary), do individual items show differential item functioning by sex, parent education, or number of languages spoken – within English and within Vietnamese separately?
  3. Part 3 – Cross-language item-level ranking. Do the specific words that are hard in English tend to be the specific words that are hard in Vietnamese (translation-equivalent pairs)? This needs the item-level EN-VI crosswalk and was the original content of this file; kept below largely as-is.
noun_kinds <- c("animals", "body_parts", "clothing", "food_drink", "toys", "vehicles",
                "household", "furniture_rooms", "outdoor")
verb_kinds <- c("action_words")

lang_config <- list(
  en = list(
    data_file   = "data/EnglishAmericanWS_Bui_data_redact.csv",
    fields_file = "data/EnglishAmericanWS_Bui_fields.csv",
    prefix      = "cat_eng_"
  ),
  vi = list(
    data_file   = "data/VietnameseWS_Bui_data.csv",
    fields_file = "data/VietnameseWS_Bui_fields.csv",
    prefix      = "cat_vie_"
  )
)

# Loads the item-response matrix (d_mat), demographics (d_demo), item
# metadata (d_items), and the original wide response data (d_wide) for one
# language. d_mat/d_demo/d_wide all share row order (all three are derived
# from the same read_csv() call on the same file), which the rest of this
# Rmd relies on for position-based subsetting -- there is no join key back to
# d_mat's rows other than row position.
#
# d_demo additionally derives n_languages_household (count of non-missing
# language_1..language_4 slots, i.e. how many languages the child is
# reported to be exposed to) and a <=2-vs->2-languages grouping from it,
# used in Part 2's demographic-bias analysis. A monolingual/multilingual (1
# vs. 2+) cut was considered first but rejected: only 11/142 children are
# truly monolingual, too few to fit a per-group model reliably. The
# distribution is 11/89/34/8 for 1/2/3/4 languages, so <=2 vs. >2 gives a
# much more workable 100/42 split.
load_lang_data <- function(lang) {
  cfg <- lang_config[[lang]]

  d_wide <- read_csv(cfg$data_file, show_col_types = FALSE) %>%
    select(response_id, starts_with(cfg$prefix)) %>%
    mutate(across(everything(), ~replace_na(.x, 0)))

  d_demo <- read_csv(cfg$data_file, show_col_types = FALSE) %>%
    select(row_id, response_id, age, sex, ppl_1_edu,
           language_1, language_2, language_3, language_4) %>%
    mutate(
      n_languages_household = rowSums(!is.na(across(c(language_1, language_2, language_3, language_4)))),
      n_languages_group = factor(if_else(n_languages_household > 2, "more_than_2", "2_or_fewer"),
                                  levels = c("2_or_fewer", "more_than_2"))
    )

  d_items <- read_csv(cfg$fields_file, show_col_types = FALSE) %>%
    filter(group == "item", type == "word") %>%
    mutate(
      definition = make.names(column),
      item_kind = str_remove(column, paste0("^", cfg$prefix)) %>% str_remove("___.*$"),
      item_definition = str_remove(column, "^.*?___"),
      domain = case_when(
        item_kind %in% verb_kinds ~ "verb",
        item_kind %in% noun_kinds ~ "noun",
        TRUE ~ "other"
      )
    ) %>%
    rename(item_id = field) %>%
    select(-c(group, type))

  d_mat <- d_wide %>% data.frame() %>% select(-response_id) %>% data.matrix()
  stopifnot(all(colnames(d_mat) == d_items$definition))

  list(d_mat = d_mat, d_demo = d_demo, d_items = d_items, d_wide = d_wide)
}

en_data <- load_lang_data("en")
vi_data <- load_lang_data("vi")

Part 1: Item difficulty trajectories by category (within language)

Two complementary views of “what gets learned first, and does it match across languages”:

  1. IRT difficulty ranking – reusing each language’s already-fit unidimensional 2PL model (models/en/mod_2pl.Rds, models/vi/mod_2pl.Rds, from 02_IRT_analysis_en.Rmd/02_IRT_analysis_vi.Rmd), extract each item’s difficulty (b) and look at how it’s distributed within and across item_kind categories. Lower b = easier / produced earlier.
  2. Empirical production-rate trajectories by age – for each item, the raw (model-free) proportion of children who produce it, plotted against age and aggregated to the category level.

IRT-based difficulty by category

get_item_difficulty <- function(mod, d_items) {
  as_tibble(coef(mod, simplify = TRUE, IRTpars = TRUE)$items) %>%
    mutate(definition = rownames(coef(mod, simplify = TRUE)$items)) %>%
    left_join(d_items, by = "definition") %>%
    select(definition, item_kind, domain, item_definition, a, b) %>%
    arrange(b) %>%
    mutate(difficulty_rank = row_number())  # 1 = easiest
}

load("models/en/mod_2pl.Rds"); mod_2pl_en <- mod_2pl; rm(mod_2pl)
load("models/vi/mod_2pl.Rds"); mod_2pl_vi <- mod_2pl; rm(mod_2pl)

diff_en <- get_item_difficulty(mod_2pl_en, en_data$d_items) %>% mutate(language = "English")
diff_vi <- get_item_difficulty(mod_2pl_vi, vi_data$d_items) %>% mutate(language = "Vietnamese")
# Category-level (item_kind) average difficulty, ranked -- "which kinds of
# words tend to be produced earliest" within each language.
category_difficulty <- bind_rows(diff_en, diff_vi) %>%
  group_by(language, item_kind, domain) %>%
  summarise(n_items = n(), mean_b = mean(b), median_b = median(b), .groups = "drop") %>%
  arrange(language, mean_b)

category_difficulty %>%
  filter(language == "English") %>%
  select(-language) %>%
  kable(digits = 2, caption = "English: item categories ranked from easiest (most negative mean b) to hardest.") %>%
  html_table_width(c(140, 80, 60, 70, 70))
English: item categories ranked from easiest (most negative mean b) to hardest.
item_kind domain n_items mean_b median_b
sounds other 12 0.91 0.88
games_routines other 25 1.38 1.48
vehicles noun 14 1.44 1.41
toys noun 18 1.46 1.57
body_parts noun 27 1.47 1.35
animals noun 43 1.49 1.45
food_drink noun 68 1.76 1.78
question_words other 7 1.81 1.73
action_words verb 103 1.90 1.90
descriptive_words other 63 1.93 1.95
outdoor noun 31 1.96 1.92
time_words other 12 2.02 2.04
clothing noun 28 2.03 2.14
people other 29 2.06 2.10
furniture_rooms noun 34 2.08 1.99
household noun 50 2.10 1.99
places other 22 2.14 2.12
pronouns other 25 2.23 2.29
locations other 26 2.37 2.30
quantifiers other 17 2.46 2.44
helping_verbs other 21 2.54 2.45
connecting_words other 6 2.70 2.71 
category_difficulty %>%
  filter(language == "Vietnamese") %>%
  select(-language) %>%
  kable(digits = 2, caption = "Vietnamese: item categories ranked from easiest (most negative mean b) to hardest.") %>%
  html_table_width(c(140, 80, 60, 70, 70))
Vietnamese: item categories ranked from easiest (most negative mean b) to hardest.
item_kind domain n_items mean_b median_b
body_parts noun 29 0.52 0.38
people other 29 0.67 1.02
animals noun 46 0.96 1.09
games_routines other 25 0.97 0.94
question_words other 7 1.04 1.24
vehicles noun 15 1.09 1.03
toys noun 19 1.13 1.18
action_words verb 103 1.14 1.13
household noun 58 1.22 1.18
outdoor noun 35 1.32 1.29
food_drink noun 69 1.39 1.58
furniture_rooms noun 33 1.42 1.39
descriptive_words other 62 1.42 1.37
locations other 24 1.48 1.57
sounds other 12 1.55 2.30
time_words other 12 1.67 1.65
clothing noun 28 1.69 1.81
places other 23 1.71 1.73
quantifiers other 16 1.76 1.66
helping_verbs other 12 1.97 1.99
pronouns other 24 2.18 2.11
connecting_words other 6 2.28 2.21

Note: b is on each language’s own, independently-calibrated scale (2PL models fit separately per language, each anchored to that language’s own sample mean-0/variance-1 identification) – a b of 1.0 in English and 1.0 in Vietnamese are not directly comparable magnitudes, only within-language relative positions are meaningful. Both plots below therefore give each language its own free b-axis (scales = "free", not "free_x") so the panels aren’t visually implying a shared scale that doesn’t exist; see the cross-language section further down for how the two languages’ difficulties are actually compared (within-language z-scoring, not raw b).

# Items ordered by difficulty within each language, colored by category -- a
# visual "trajectory" of which categories cluster early (easy/low b) vs. late
# (hard/high b) in the difficulty-implied developmental order.
bind_rows(diff_en, diff_vi) %>%
  ggplot(aes(x = difficulty_rank, y = b, color = item_kind)) +
  geom_point(alpha = 0.6, size = 1) +
  facet_wrap(~language, scales = "free") +
  theme_classic() +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(nrow = 3)) +
  xlab("Difficulty rank (1 = easiest)") +
  ylab("IRT difficulty (b)")

# Same information as a category-ordered boxplot (categories ordered by
# median difficulty within each language) -- easier to read the within-category
# spread than the scatter above.
bind_rows(diff_en, diff_vi) %>%
  mutate(item_kind = fct_reorder(item_kind, b, .fun = median)) %>%
  ggplot(aes(x = item_kind, y = b, fill = domain)) +
  geom_boxplot(outlier.size = 0.5) +
  facet_wrap(~language, ncol = 1, scales = "free") +
  coord_flip() +
  theme_classic() +
  xlab(NULL) + ylab("IRT difficulty (b)")

Age of acquisition (AoA) by item

A second, independent difficulty measure: each item’s empirical age of acquisition (AoA), estimated the standard Wordbank/Braginsky-et-al. way – fit produces ~ age as a logistic regression per item, and take AoA as the age at which the fitted curve crosses 50% production (-intercept/slope). This only uses each item’s own production-by-age pattern, with no IRT model involved, so it’s a useful independent check on whether the IRT difficulty (b) computed above is actually capturing “when children learn this word” or something else.

Coverage note, checked empirically: with only ~140 children per language and a ~14-36 month age range, most items’ 50%-crossing point falls outside the observed age range (i.e. the item is already produced by >50% of even the youngest sampled children, or still produced by <50% of the oldest) – extrapolating the logistic curve to estimate AoA anyway would be unreliable, so those items get AoA = NA rather than a number. In practice this leaves only ~20% of items with an estimable AoA (146/681 for English in a quick check), so the tables and correlation below necessarily cover a minority of the full item set, not “all words.”

# Items with no variance (nobody, or everybody, in the sample produces
# them) can't have a logistic curve fit and get AoA = NA; same for items
# where glm() doesn't converge, returns a non-positive age slope (i.e.
# production doesn't increase with age -- shouldn't happen for real
# vocabulary but guards against a degenerate fit), or -- the main source of
# NAs in practice -- where the fitted 50%-crossing point falls outside the
# sample's observed age range (see the coverage note above: extrapolating
# there would be unreliable rather than merely uncertain).
compute_item_aoa <- function(d_wide, d_demo, d_items) {
  age_range <- range(d_demo$age, na.rm = TRUE)

  d_long <- d_wide %>%
    left_join(d_demo %>% select(response_id, age), by = "response_id") %>%
    pivot_longer(cols = -c(response_id, age), names_to = "definition", values_to = "produces")

  items <- unique(d_long$definition)
  aoa_vals <- map_dbl(items, function(def) {
    d_i <- d_long %>% filter(definition == def)
    if (sum(d_i$produces) == 0 || sum(d_i$produces) == nrow(d_i)) return(NA_real_)
    fit <- tryCatch(glm(produces ~ age, data = d_i, family = binomial), error = function(e) NULL)
    if (is.null(fit) || !fit$converged) return(NA_real_)
    b0 <- unname(coef(fit)[1]); b1 <- unname(coef(fit)[2])
    if (b1 <= 0) return(NA_real_)
    aoa_est <- -b0 / b1
    if (aoa_est < age_range[1] || aoa_est > age_range[2]) return(NA_real_)
    aoa_est
  })

  tibble(definition = items, aoa = aoa_vals) %>%
    left_join(d_items %>% select(definition, item_kind, domain, item_definition), by = "definition")
}

aoa_en <- compute_item_aoa(en_data$d_wide, en_data$d_demo, en_data$d_items)
aoa_vi <- compute_item_aoa(vi_data$d_wide, vi_data$d_demo, vi_data$d_items)

cat("English: AoA estimated for", sum(!is.na(aoa_en$aoa)), "/", nrow(aoa_en), "items\n")
## English: AoA estimated for 146 / 681 items
cat("Vietnamese: AoA estimated for", sum(!is.na(aoa_vi$aoa)), "/", nrow(aoa_vi), "items\n")
## Vietnamese: AoA estimated for 380 / 687 items
aoa_en %>%
  filter(!is.na(aoa)) %>%
  arrange(aoa) %>%
  slice(c(1:15, (n() - 14):n())) %>%
  select(item_definition, item_kind, aoa) %>%
  kable(digits = 1, caption = "English: 15 earliest- and 15 latest-acquired words by empirical AoA (months).") %>%
  html_table_width(c(150, 100, 80))
English: 15 earliest- and 15 latest-acquired words by empirical AoA (months).
item_definition item_kind aoa
meow sounds 19.3
bye games_routines 19.5
hi games_routines 21.6
daddy people 22.0
apple food_drink 22.9
car vehicles 22.9
mommy people 23.9
banana food_drink 24.2
NA NA 24.7
eye body_parts 25.0
fish animals 26.1
dog animals 26.5
bus vehicles 27.1
duck animals 27.1
nose body_parts 27.2
ice food_drink 35.5
bottle household 35.5
pancake food_drink 35.6
wait action_words 35.6
green descriptive_words 35.6
you pronouns 35.7
alligator animals 35.7
NA NA 35.7
brown descriptive_words 35.7
cheese food_drink 35.8
box household 35.9
sleep action_words 35.9
moon outdoor 35.9
popsicle food_drink 36.0
tummy body_parts 36.0 
aoa_vi %>%
  filter(!is.na(aoa)) %>%
  arrange(aoa) %>%
  slice(c(1:15, (n() - 14):n())) %>%
  select(item_definition, item_kind, aoa) %>%
  kable(digits = 1, caption = "Vietnamese: 15 earliest- and 15 latest-acquired words by empirical AoA (months).") %>%
  html_table_width(c(150, 100, 80))
Vietnamese: 15 earliest- and 15 latest-acquired words by empirical AoA (months).
item_definition item_kind aoa
mũi body_parts 16.1
ông people 16.2
NA NA 17.3
mèo animals 17.7
mắt body_parts 17.8
chó animals 18.4
cá food_drink 18.7
đầu body_parts 18.7
tắm games_routines 18.8
chuối food_drink 19.1
tai body_parts 19.2
NA NA 19.4
chân body_parts 19.6
nước food_drink 19.6
tóc body_parts 19.8
ngoài places 35.5
NA NA 35.5
NA NA 35.6
NA NA 35.6
cao descriptive_words 35.6
che action_words 35.6
mềm descriptive_words 35.7
dây household 35.7
làm helping_verbs 35.8
bắt action_words 35.8
gõ action_words 35.9
NA NA 35.9
sói animals 35.9
hươu animals 36.0
nghĩ action_words 36.0

Does AoA agree with IRT difficulty?

Correlated separately per language (not pooled) – aoa (age, months) and b (IRT difficulty, logits) are on different, non-linearly-related scales (the logistic curve fit for AoA is itself a nonlinear function of age), so a rank-based Spearman correlation is used rather than Pearson. This is a within-language comparison (AoA vs. b for the same item, same language), unlike Part 3’s cross-language b-vs-b comparison below, which instead z-scores and uses Pearson since both quantities there are already on a (standardized) linear scale.

aoa_diff_en <- aoa_en %>% inner_join(diff_en %>% select(definition, b), by = "definition") %>% filter(!is.na(aoa))
aoa_diff_vi <- aoa_vi %>% inner_join(diff_vi %>% select(definition, b), by = "definition") %>% filter(!is.na(aoa))

cat("=== English: AoA vs. IRT difficulty (b) ===\n")
## === English: AoA vs. IRT difficulty (b) ===
print(cor.test(aoa_diff_en$aoa, aoa_diff_en$b, method = "spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  aoa_diff_en$aoa and aoa_diff_en$b
## S = 52284, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.8724905
cat("\n=== Vietnamese: AoA vs. IRT difficulty (b) ===\n")
## 
## === Vietnamese: AoA vs. IRT difficulty (b) ===
print(cor.test(aoa_diff_vi$aoa, aoa_diff_vi$b, method = "spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  aoa_diff_vi$aoa and aoa_diff_vi$b
## S = 56488, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.9780459
bind_rows(aoa_diff_en %>% mutate(language = "English"),
          aoa_diff_vi %>% mutate(language = "Vietnamese")) %>%
  ggplot(aes(x = aoa, y = b, color = domain)) +
  geom_point(alpha = 0.5, size = 1) +
  geom_smooth(method = "lm", se = TRUE, color = "grey30") +
  facet_wrap(~language, scales = "free") +
  theme_classic() +
  xlab("Empirical age of acquisition (months)") +
  ylab("IRT difficulty (b)")
## `geom_smooth()` using formula = 'y ~ x'

A strong positive correlation in both languages would confirm b is capturing the same underlying “when is this word learned” construct as the model-free AoA estimate; a weak or inconsistent correlation would suggest the IRT model’s difficulty ordering is picking up something else (e.g. discrimination/a absorbing some of what a single-parameter difficulty can’t, since this is the 2PL model, not Rasch).

Are the same categories easy/hard in both languages?

item_kind categories are shared between the English and Vietnamese fields files (unlike individual items, which need the translation crosswalk used in Part 3), so category-average difficulty can be compared between languages directly – no item-level translation matching required. As noted above, raw b isn’t on a shared scale between the two independently-calibrated models, so mean_b is z-scored within each language (across that language’s own categories) before comparing – the same within-language standardization used for ability_noun in 03_bifactor_noun_verb.Rmd, for the same reason (independently-calibrated models, arbitrary relative scale). Z-scoring (unlike converting to ranks) preserves relative spacing, not just order, so the comparison below uses a Pearson correlation on the z-scores rather than a rank-based Spearman correlation.

Results: Moderate positive correlation between English and Vietnamese.

category_wide <- category_difficulty %>%
  select(language, item_kind, domain, mean_b) %>%
  pivot_wider(names_from = language, values_from = mean_b) %>%
  filter(!is.na(English), !is.na(Vietnamese)) %>%
  mutate(
    en_z = as.numeric(scale(English)),
    vi_z = as.numeric(scale(Vietnamese))
  )

cor.test(category_wide$en_z, category_wide$vi_z, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  category_wide$en_z and category_wide$vi_z
## t = 3.4533, df = 20, p-value = 0.002512
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2553759 0.8211860
## sample estimates:
##       cor 
## 0.6111795
# geom_text_repel() needs the ggrepel package (not currently a project
# dependency); if unavailable, swap in plain geom_text(size = 3) -- labels
# will just overlap more.
ggplot(category_wide, aes(x = en_z, y = vi_z, label = item_kind, color = domain)) +
  geom_point() +
  ggrepel::geom_text_repel(size = 3, show.legend = FALSE) +
  geom_smooth(method = "lm", se = TRUE, color = "grey40") +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey70") +
  theme_classic() +
  xlab("English category difficulty (z-scored within language)") +
  ylab("Vietnamese category difficulty (z-scored within language)")
## `geom_smooth()` using formula = 'y ~ x'
## Warning: The following aesthetics were dropped during statistical transformation: label.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

Categories off the diagonal are inconsistent (easy in one language, hard in the other); categories on the diagonal are candidates for a genuinely conceptually-driven (language-independent) acquisition order.

*Results: sounds, pronouns, connecting words, clothing, time words relatively easier in Vietnamese; body parts, people, question words, action words, household relatively easier in English.

Empirical production-rate trajectories by age and category

A second, model-free view: for each item, the raw proportion of children who produce it, in 2-month age bins, aggregated up to item_kind – the standard Wordbank-style “vocabulary growth curve,” aggregated by category instead of shown item-by-item.

compute_category_trajectory <- function(d_wide, d_demo, d_items, age_bin_width = 2) {
  d_wide %>%
    left_join(d_demo %>% select(response_id, age), by = "response_id") %>%
    pivot_longer(cols = -c(response_id, age), names_to = "definition", values_to = "produces") %>%
    left_join(d_items %>% select(definition, item_kind, domain), by = "definition") %>%
    mutate(age_bin = floor(age / age_bin_width) * age_bin_width) %>%
    group_by(item_kind, domain, age_bin) %>%
    summarise(prop_producing = mean(produces), n_children = n_distinct(response_id), .groups = "drop")
}

traj_en <- compute_category_trajectory(en_data$d_wide, en_data$d_demo, en_data$d_items) %>%
  mutate(language = "English")
traj_vi <- compute_category_trajectory(vi_data$d_wide, vi_data$d_demo, vi_data$d_items) %>%
  mutate(language = "Vietnamese")
bind_rows(traj_en, traj_vi) %>%
  ggplot(aes(x = age_bin, y = prop_producing, color = item_kind)) +
  geom_line(alpha = 0.7) +
  facet_wrap(~language, ncol = 1) +
  theme_classic() +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(nrow = 3)) +
  xlab("Age (months, binned)") +
  ylab("Mean proportion of items in category produced")

Categories whose curves rise earliest (leftmost) are learned first; comparing panels shows whether the same categories rise early in both languages (conceptual/maturational account) or whether the rise order differs by language (input-driven account).

Quantifying category learning order: age at 50% production

To turn “which curve rises earliest” into a number, take each category’s age-binned prop_producing curve (just computed above) and find the age at which it crosses 50%, linearly interpolating between the two bins straddling 0.5 – the category-level analogue of the item-level AoA computed earlier, but using the already-aggregated category curve directly rather than needing an item-by-item logistic fit. As with item-level AoA, a category whose curve never reaches 50% within the observed age range (still climbing, or – rarer – already above 50% at the very first bin, i.e. crossed before the observed window) gets NA rather than an extrapolated guess.

Unlike the IRT b-based category comparison above, this crossing age is in months – a real, shared unit across both surveys (the same clock, not an independently-calibrated latent scale) – so it can be compared directly between languages without z-scoring; a raw difference in months is directly interpretable.

category_crossing_age <- function(traj_df, threshold = 0.25) {
  traj_df %>%
    group_by(item_kind, domain) %>%
    arrange(age_bin, .by_group = TRUE) %>%
    summarise(
      crossing_age = {
        idx <- which(prop_producing >= threshold)[1]
        if (is.na(idx) || idx == 1) {
          NA_real_  # never reaches threshold in-range, or already >= threshold at the first observed bin
        } else {
          age_bin[idx - 1] + (threshold - prop_producing[idx - 1]) /
            (prop_producing[idx] - prop_producing[idx - 1]) * (age_bin[idx] - age_bin[idx - 1])
        }
      },
      .groups = "drop"
    )
}

category_crossing_en <- category_crossing_age(traj_en) %>% mutate(language = "English")
category_crossing_vi <- category_crossing_age(traj_vi) %>% mutate(language = "Vietnamese")

category_crossing_en %>%
  arrange(crossing_age) %>%
  select(-language) %>%
  kable(digits = 1, caption = "English: categories ordered by age (months) at 50% production.") %>%
  html_table_width(c(140, 80, 90))
English: categories ordered by age (months) at 50% production.
item_kind domain crossing_age
sounds other 20.3
games_routines other 21.2
body_parts noun 21.5
vehicles noun 22.8
animals noun 23.3
toys noun 25.6
food_drink noun 28.4
people other 28.8
places other 29.0
outdoor noun 29.4
clothing noun 29.5
pronouns other 29.9
furniture_rooms noun 29.9
NA NA 30.1
question_words other 30.4
descriptive_words other 30.5
action_words verb 30.8
time_words other 34.6
connecting_words other 35.6
household noun 35.6
helping_verbs other NA
locations other NA
quantifiers other NA 
category_crossing_vi %>%
  arrange(crossing_age) %>%
  select(-language) %>%
  kable(digits = 1, caption = "Vietnamese: categories ordered by age (months) at 50% production.") %>%
  html_table_width(c(140, 80, 90))
Vietnamese: categories ordered by age (months) at 50% production.
item_kind domain crossing_age
sounds other 14.6
body_parts noun 17.6
clothing noun 19.3
food_drink noun 20.2
animals noun 20.3
toys noun 20.9
action_words verb 21.2
outdoor noun 21.2
furniture_rooms noun 21.2
household noun 22.1
quantifiers other 22.4
places other 22.5
time_words other 22.6
locations other 22.7
descriptive_words other 22.9
helping_verbs other 23.1
NA NA 23.8
pronouns other 24.0
connecting_words other 25.9
games_routines other NA
people other NA
question_words other NA
vehicles noun NA 
category_crossing_wide <- bind_rows(category_crossing_en, category_crossing_vi) %>%
  select(language, item_kind, domain, crossing_age) %>%
  pivot_wider(names_from = language, values_from = crossing_age) %>%
  filter(!is.na(English), !is.na(Vietnamese)) %>%
  mutate(diff_vi_minus_en = Vietnamese - English)

cor.test(category_crossing_wide$English, category_crossing_wide$Vietnamese, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  category_crossing_wide$English and category_crossing_wide$Vietnamese
## t = 4.9217, df = 14, p-value = 0.000225
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4961607 0.9262596
## sample estimates:
##       cor 
## 0.7960724
category_crossing_wide %>%
  arrange(desc(abs(diff_vi_minus_en))) %>%
  kable(digits = 1, caption = "Category 25%-production age (months), English vs. Vietnamese, sorted by |difference|.") %>%
  html_table_width(c(140, 60, 90, 90, 90))
Category 25%-production age (months), English vs. Vietnamese, sorted by |difference|.
item_kind domain English Vietnamese diff_vi_minus_en
household noun 35.6 22.1 -13.5
time_words other 34.6 22.6 -12.0
clothing noun 29.5 19.3 -10.2
connecting_words other 35.6 25.9 -9.7
action_words verb 30.8 21.2 -9.7
furniture_rooms noun 29.9 21.2 -8.7
food_drink noun 28.4 20.2 -8.2
outdoor noun 29.4 21.2 -8.2
descriptive_words other 30.5 22.9 -7.6
places other 29.0 22.5 -6.5
NA NA 30.1 23.8 -6.3
pronouns other 29.9 24.0 -5.9
sounds other 20.3 14.6 -5.7
toys noun 25.6 20.9 -4.6
body_parts noun 21.5 17.6 -3.9
animals noun 23.3 20.3 -3.0 
ggplot(category_crossing_wide, aes(x = English, y = Vietnamese, label = item_kind, color = domain)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey70") +
  geom_point() +
  ggrepel::geom_text_repel(size = 3, show.legend = FALSE) +
  theme_classic() +
  xlab("English: age (months) at 50% production") +
  ylab("Vietnamese: age (months) at 50% production")
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text_repel()`).

A dumbbell plot of the same numbers, ordered by how far apart the two languages are for that category – categories near the top show the largest EN-VI gap in learning order; the dashed connecting line points from whichever language is earlier to whichever is later.

category_crossing_wide %>%
  mutate(item_kind = fct_reorder(item_kind, abs(diff_vi_minus_en))) %>%
  ggplot(aes(y = item_kind)) +
  geom_segment(aes(x = English, xend = Vietnamese, y = item_kind, yend = item_kind),
               color = "grey70", linewidth = 0.6) +
  geom_point(aes(x = English, color = "English"), size = 2.5) +
  geom_point(aes(x = Vietnamese, color = "Vietnamese"), size = 2.5) +
  theme_classic() +
  xlab("Age (months) at 50% production") + ylab(NULL) +
  labs(color = NULL)

Part 2: Demographic bias (Differential Item Functioning) analysis

Following Kachergis, Francis, & Frank (2022, Estimating demographic bias on tests of children’s early vocabulary), we test items for DIF along three demographic dimensions, separately for English and Vietnamese:

Methods, following the paper

  1. Baseline Rasch (1PL) model (already fit in 02_IRT_analysis_en.Rmd/02_IRT_analysis_vi.Rmd, models/{en,vi}/mod_1pl.Rds), fit without regard to demographic group, gives each child a baseline ability theta. Regressing theta ~ age_c * group per demographic dimension characterizes the overall size of the demographic effect (paper’s “Demographic effects in the baseline Rasch model”). This step is just for the overall-effect-size context, not the item-level DIF test.
  2. Per-group Rasch models: rather than a multigroup mirt model, we fit two separate single-group Rasch models, one per demographic group. A single-group mirt() fit fixes that group’s latent mean at 0 by default, so fitting each group separately reproduces the paper’s identification assumption (“the mean language ability in each group is the same… thus pushing all between-group variation into the item [parameters]”). Rasch (1PL, a fixed at 1 for every item) is used throughout, matching the paper – it’s the simplest model and the most robust choice at our modest per-group sample sizes (~50-100 children, far smaller than the paper’s ~2000+ per group), and with Rasch a group difference is unambiguously a difficulty (b) shift (“uniform DIF”).
  3. For each item, b_diff = b_group1 - b_group2 is the difficulty difference between groups – the “uniform DIF” signal, and the only kind Rasch can detect (discrimination is fixed at 1 for every item, in every group). We summarize its distribution (median/mean/SD, count of items favoring each group) and visualize it as an ordered dot plot – a simplified, point-estimate-only version of the paper’s GLIMMER plots (which additionally draw 10,000 imputations from each item’s parameter covariance matrix to show uncertainty; that needs SE = TRUE when fitting mirt(), which is slower, so it’s omitted here – see the note in dif_dotplot() below).
  4. Pruning: at a range of |b_diff| exclusion thresholds (in SD units, 0.25 to 3), drop the most extreme items, refit a single-group Rasch model on the retained items (full sample, not split by group), and re-run the theta ~ age_c * group regression – replicating the paper’s Figure 3 (does trimming biased items shrink the estimated demographic effect?).
  5. Extreme items: items beyond a chosen threshold (paper settled on ~2.25 SD), listed with item_kind/ domain for qualitative inspection, as in the paper’s Figure 4.

Before running: things to verify

  • ppl_1_edu high/low split. Rather than the paper’s “any college” cutoff (which here would give a lopsided 130 high / 9 low split, too few “low” children to fit reliably), edu_group below is a median split on an 8-level ordinal scale (grade-8-or-less to master’s/doctoral, EN/VI label pairs collapsed to the same rank). Verified against both data files: identical distribution in the English and Vietnamese CDI (shared demographic field), median rank 6 (“Bachelor’s degree”), giving 109 high / 27 low. edu_na_levels (“Không rõ”/“Not sure”) are treated as missing, not “low”, since they’re non-response rather than a reported attainment level.
  • sex level names/casing. sex_group below assumes levels exactly "M"/"F" (confirmed against the raw data). If that’s ever wrong – e.g. after a data refresh changes the coding – factor(sex, levels = c("M","F")) would silently produce NA for every row and the whole sex-DIF analysis would silently run on zero data. The diagnostic table() calls in add_dif_groups() below are there specifically to catch this before it happens quietly – check their output before trusting anything downstream.
  • Row-order alignment. fit_group_irt() and prune_and_refit() subset d_mat by row position to match d_demo’s row order (there’s no join key back to d_mat’s rows). This matches how d_mat/d_demo are constructed in load_lang_data() above (same file, same read_csv() row order) and how the rest of this project’s Rmds rely on the same invariant – but it would silently break if d_demo were ever re-sorted without re-sorting d_mat to match.
  • ggrepel is used in Part 1’s plot only; not needed here.

Demographic group definitions

# Ordinal 8-level parent-1 education scale (EN/VI label pairs collapsed to
# the same rank), used for a median split into edu_group. edu_na_levels
# ("Không rõ"/"Not sure") are treated as missing, not "low", since they're
# non-response rather than a reported attainment level.
#
# Verified against both data files (identical distribution in the English
# and Vietnamese CDI, since ppl_1_edu is a shared demographic field): n =
# 136 non-missing / 6 missing, median rank = 6 ("Four-year college
# Bachelor's degree" / "Có bằng cử nhân bốn năm"), giving a 109 (high) / 27
# (low) split at the median -- far more balanced for per-group model fitting
# than a naive "any college" cutoff, which left only 6-9 children in "low".
edu_levels_ordered <- list(
  "1_grade8_or_less"      = c("Lớp 8 hoặc ít hơn"),
  "2_some_high_school"    = c("Vài năm trung học phổ thông (cấp 3)"),
  "3_high_school_grad"    = c("Hoàn thành trung học phổ thông (cấp 3)", "High school graduate or GED"),
  "4_some_college"        = c("One or more years of college, no degree", "Vài năm cao đẳng, không có bằng"),
  "5_two_year_college"    = c("Two-year college degree or vocational degree",
                               "Có bằng cao đẳng hoặc tốt nghiệp trường dạy nghề"),
  "6_bachelors"           = c("Four-year college Bachelor's degree", "Có bằng cử nhân bốn năm"),
  "7_some_grad_school"    = c("Vài năm cao học", "Some graduate education"),
  "8_masters_or_doctoral" = c("Master's or doctoral degree", "Có bằng thạc sĩ, tiến sĩ hoặc tương tự")
)
edu_na_levels    <- c("Không rõ", "Not sure")
edu_rank_lookup  <- setNames(rep(seq_along(edu_levels_ordered), lengths(edu_levels_ordered)),
                              unlist(edu_levels_ordered))
edu_median_rank  <- 6  # see note above; re-check if the sample changes

add_dif_groups <- function(d_demo, lang_label) {
  cat("--", lang_label, "sex distribution --\n")
  print(table(d_demo$sex, useNA = "always"))
  cat("\n--", lang_label, "ppl_1_edu distribution --\n")
  print(table(d_demo$ppl_1_edu, useNA = "always"))
  cat("\n--", lang_label, "n_languages_household distribution --\n")
  print(table(d_demo$n_languages_household, useNA = "always"))

  edu_rank <- if_else(d_demo$ppl_1_edu %in% edu_na_levels,
                       NA_integer_, unname(edu_rank_lookup[d_demo$ppl_1_edu]))

  d_demo %>%
    mutate(
      sex_group  = factor(sex, levels = c("M", "F")),
      edu_rank   = edu_rank,
      edu_group  = factor(if_else(edu_rank >= edu_median_rank, "high", "low"),
                           levels = c("low", "high")),
      lang_group = n_languages_group
    )
}

en_data$d_demo <- add_dif_groups(en_data$d_demo, "English")
## -- English sex distribution --
## 
##    F    M <NA> 
##   65   77    0 
## 
## -- English ppl_1_edu distribution --
## 
## Có bằng cao đẳng hoặc tốt nghiệp trường dạy nghề 
##                                                6 
##                          Có bằng cử nhân bốn năm 
##                                               38 
##           Có bằng thạc sĩ, tiến sĩ hoặc tương tự 
##                                               21 
##              Four-year college Bachelor's degree 
##                                               19 
##                      High school graduate or GED 
##                                                2 
##           Hoàn thành trung học phổ thông (cấp 3) 
##                                                2 
##                                         Không rõ 
##                                                2 
##                                Lớp 8 hoặc ít hơn 
##                                                1 
##                      Master's or doctoral degree 
##                                               28 
##                                         Not sure 
##                                                1 
##          One or more years of college, no degree 
##                                                6 
##                          Some graduate education 
##                                                1 
##     Two-year college degree or vocational degree 
##                                                2 
##                  Vài năm cao đẳng, không có bằng 
##                                                7 
##                                  Vài năm cao học 
##                                                2 
##              Vài năm trung học phổ thông (cấp 3) 
##                                                1 
##                                             <NA> 
##                                                3 
## 
## -- English n_languages_household distribution --
## 
##    1    2    3    4 <NA> 
##   11   89   34    8    0
vi_data$d_demo <- add_dif_groups(vi_data$d_demo, "Vietnamese")
## -- Vietnamese sex distribution --
## 
##    F    M <NA> 
##   65   77    0 
## 
## -- Vietnamese ppl_1_edu distribution --
## 
## Có bằng cao đẳng hoặc tốt nghiệp trường dạy nghề 
##                                                6 
##                          Có bằng cử nhân bốn năm 
##                                               38 
##           Có bằng thạc sĩ, tiến sĩ hoặc tương tự 
##                                               21 
##              Four-year college Bachelor's degree 
##                                               19 
##                      High school graduate or GED 
##                                                2 
##           Hoàn thành trung học phổ thông (cấp 3) 
##                                                2 
##                                         Không rõ 
##                                                2 
##                                Lớp 8 hoặc ít hơn 
##                                                1 
##                      Master's or doctoral degree 
##                                               28 
##                                         Not sure 
##                                                1 
##          One or more years of college, no degree 
##                                                6 
##                          Some graduate education 
##                                                1 
##     Two-year college degree or vocational degree 
##                                                2 
##                  Vài năm cao đẳng, không có bằng 
##                                                7 
##                                  Vài năm cao học 
##                                                2 
##              Vài năm trung học phổ thông (cấp 3) 
##                                                1 
##                                             <NA> 
##                                                3 
## 
## -- Vietnamese n_languages_household distribution --
## 
##    1    2    3    4 <NA> 
##   11   89   34    8    0

Helper functions

# Baseline (single-group, whole sample) Rasch ability regressed against age
# and a demographic group -- the paper's "Demographic effects in the
# baseline Rasch model" analysis. Always uses the baseline 1PL model
# (models/{en,vi}/mod_1pl.Rds), since this step is just about the overall
# demographic effect size, not item-level DIF.
# Assumes fscores(mod_1pl) returns one row per child in the same order as
# d_demo (true if mod_1pl.Rds was fit on the same file/row order as
# load_lang_data() reads here).
baseline_group_effect <- function(mod_1pl, d_demo, group_var) {
  theta <- fscores(mod_1pl, method = "MAP")[, 1]
  d <- d_demo %>%
    mutate(theta = theta, age_c = age - mean(age, na.rm = TRUE)) %>%
    filter(!is.na(.data[[group_var]]))
  fit <- lm(as.formula(paste0("theta ~ age_c * ", group_var)), data = d)
  list(data = d, fit = fit)
}

# Items with zero variance within a given row-subset (nobody, or everybody,
# in that subset endorses it) carry no information to estimate a difficulty
# from and make mirt() error out ("only one response category and cannot be
# estimated"). This happens routinely once the sample is split into small
# demographic subgroups (some rare/near-universal CDI words end up perfectly
# 0/1 within a ~27-100-child subgroup purely by chance, even though they
# have variance in the full ~140-child sample). Returns the column names of
# d_mat_subset that have at least one 0 and at least one 1.
usable_items <- function(d_mat_subset) {
  col_sums <- colSums(d_mat_subset)
  colnames(d_mat_subset)[col_sums > 0 & col_sums < nrow(d_mat_subset)]
}

# Fits a single-group Rasch model separately within each level of group_var
# (subsetting d_mat by row position via which(d_demo[[group_var]] == g)),
# returning one item-parameter tibble per group. Each fit has that
# subgroup's own latent mean fixed at 0, reproducing the paper's
# identification assumption. Items with zero variance *within that specific
# group* are dropped before fitting (see usable_items()); different groups
# can end up dropping different items, which is fine -- compute_dif() below
# only compares items present in both groups' output (via inner_join).
fit_group_irt <- function(d_mat, d_demo, group_var, ncycles = 2000, seed = 1234) {
  groups <- levels(d_demo[[group_var]])
  out <- list()
  for (g in groups) {
    idx <- which(d_demo[[group_var]] == g)
    d_mat_g <- d_mat[idx, , drop = FALSE]
    keep <- usable_items(d_mat_g)
    if (length(keep) < ncol(d_mat_g)) {
      message(sprintf("fit_group_irt(): dropping %d zero-variance item(s) for group '%s' (n = %d)",
                       ncol(d_mat_g) - length(keep), g, length(idx)))
    }
    d_mat_g <- d_mat_g[, keep, drop = FALSE]
    set.seed(seed)
    mod_g <- mirt(d_mat_g, 1, itemtype = "Rasch",
                   technical = list(NCYCLES = ncycles), verbose = FALSE)
    cf <- coef(mod_g, simplify = TRUE, IRTpars = TRUE)$items
    out[[g]] <- tibble(definition = rownames(cf), b = as.numeric(cf[, "b"]), n = length(idx))
  }
  out
}

# Item-level difficulty differences (group1 - group2) between two per-group
# Rasch fits from fit_group_irt(), with item_kind/domain attached and
# outlier flags at a given SD threshold. b_diff is the "uniform DIF" signal
# (positive b_diff means group1's b is higher/harder, i.e. the item is
# easier for group2) -- the only kind of DIF Rasch can detect, since
# discrimination is fixed at 1 for every item, in every group.
compute_dif <- function(group_fits, group1, group2, d_items, sd_threshold = 2.25) {
  d <- group_fits[[group1]] %>%
    select(definition, b1 = b) %>%
    inner_join(group_fits[[group2]] %>% select(definition, b2 = b), by = "definition") %>%
    mutate(b_diff = b1 - b2) %>%
    left_join(d_items %>% select(definition, item_kind, domain, item_definition), by = "definition")

  m_b <- mean(d$b_diff); s_b <- sd(d$b_diff)
  d %>% mutate(
    b_diff_z = (b_diff - m_b) / s_b,
    extreme  = abs(b_diff_z) > sd_threshold,
    favors   = if_else(b_diff > 0, group2, group1)
  )
}

# Simplified GLIMMER-style plot: the top_n most extreme items ordered by
# |b_diff|, colored by which group finds them easier. The original GLIMMER
# additionally draws 10,000 imputations per item from the item parameter
# covariance matrix to show uncertainty around each b_diff -- that needs
# SE = TRUE when fitting mirt() in fit_group_irt() (slower) plus
# vcov(mod_g) per item; omitted here for a faster point-estimate-only
# version, but worth adding if the point estimates below suggest clusters
# worth confirming with uncertainty.
dif_dotplot <- function(dif_tab, group1, group2, top_n = 20) {
  dif_tab %>%
    slice_max(order_by = abs(b_diff), n = top_n) %>%
    mutate(item_definition = fct_reorder(item_definition, b_diff)) %>%
    ggplot(aes(x = b_diff, y = item_definition, color = favors)) +
    geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
    geom_point(size = 2) +
    theme_classic() +
    labs(x = paste0("Difficulty difference (", group1, " b - ", group2, " b)"),
         y = NULL, color = "Favors")
}

# Pruning sweep: at each SD threshold, drop items with |b_diff_z| > threshold,
# refit a single-group Rasch model on the RETAINED items using the full
# (both-groups-combined) sample, and re-run the theta ~ age_c * group
# regression -- replicating the paper's Figure 3.
prune_and_refit <- function(d_mat, d_demo, dif_tab, group_var, thresholds = seq(0.25, 3, by = 0.25),
                             ncycles = 2000, seed = 1234) {
  valid_idx    <- which(!is.na(d_demo[[group_var]]))
  d_demo_valid <- d_demo[valid_idx, ]
  d_mat_valid  <- d_mat[valid_idx, , drop = FALSE]

  usable <- usable_items(d_mat_valid)  # see fit_group_irt()'s note; unlikely but possible even in the combined sample

  map_dfr(thresholds, function(thr) {
    keep <- intersect(dif_tab$definition[abs(dif_tab$b_diff_z) <= thr], colnames(d_mat_valid))
    keep <- intersect(keep, usable)
    if (length(keep) < 10) {
      return(tibble(threshold = thr, n_excluded = ncol(d_mat_valid) - length(keep),
                     beta_group = NA_real_, p_value = NA_real_))
    }
    set.seed(seed)
    mod_thr <- mirt(d_mat_valid[, keep], 1, itemtype = "Rasch",
                     technical = list(NCYCLES = ncycles), verbose = FALSE)
    theta <- fscores(mod_thr, method = "MAP")[, 1]
    d <- d_demo_valid %>% mutate(theta = theta, age_c = age - mean(age, na.rm = TRUE))
    fit <- lm(as.formula(paste0("theta ~ age_c * ", group_var)), data = d)
    cf <- summary(fit)$coefficients
    group_row <- grep(paste0("^", group_var), rownames(cf), value = TRUE)[1]
    tibble(threshold = thr, n_excluded = ncol(d_mat_valid) - length(keep),
           beta_group = cf[group_row, "Estimate"], p_value = cf[group_row, "Pr(>|t|)"])
  })
}

Run: every (language x demographic) combination

Loops over all 6 combinations (2 languages x 3 demographics), caching each combination’s per-group Rasch fits and pruning sweep to models/{en,vi}/dif_<group_var>_Rasch_{fits,prune}.Rds (same on-disk caching convention as bootstrap_noun_ci() in 03_bifactor_noun_verb.Rmd – delete a cache file to force a recompute for that combination). Uses results='asis' so each iteration can print its own subheading, tables, and plots.

dif_specs <- tribble(
  ~language, ~group_var,   ~group1,       ~group2,
  "en",      "sex_group",  "M",           "F",
  "vi",      "sex_group",  "M",           "F",
  "en",      "edu_group",  "low",         "high",
  "vi",      "edu_group",  "low",         "high",
  "en",      "lang_group", "2_or_fewer",  "more_than_2",
  "vi",      "lang_group", "2_or_fewer",  "more_than_2"
)

lang_data  <- list(en = en_data, vi = vi_data)
lang_label <- c(en = "English", vi = "Vietnamese")

# mod_1pl.Rds was written via save(file = ..., "mod_1pl", "fscores_1pl", "coefs_1pl")
# in 02_IRT_analysis_en.Rmd/02_IRT_analysis_vi.Rmd -- a multi-object save()
# file, not a single-object saveRDS() one, so it needs load() (not
# readRDS(), which errors with "unknown input format" on this file) and
# each load() populates mod_1pl/fscores_1pl/coefs_1pl into the environment,
# overwriting the previous language's copy -- hence capturing mod_1pl into a
# renamed variable immediately after each load(), same as mod_2pl_en/
# mod_2pl_vi above.
load("models/en/mod_1pl.Rds"); mod_1pl_en <- mod_1pl; rm(mod_1pl)
load("models/vi/mod_1pl.Rds"); mod_1pl_vi <- mod_1pl; rm(mod_1pl)
rm(fscores_1pl, coefs_1pl)
mod_1pl <- list(en = mod_1pl_en, vi = mod_1pl_vi)

dif_results <- list()

for (i in seq_len(nrow(dif_specs))) {
  spec <- dif_specs[i, ]
  ld   <- lang_data[[spec$language]]
  key  <- paste(spec$language, spec$group_var, sep = "_")

  cat("\n\n### ", lang_label[[spec$language]], ": ", spec$group_var, "\n\n", sep = "")

  base <- baseline_group_effect(mod_1pl[[spec$language]], ld$d_demo, spec$group_var)
  cat("Baseline demographic effect (`theta ~ age_c * group`):\n\n")
  print(knitr::kable(broom::tidy(base$fit), digits = 3))

  fits_cache <- paste0("models/", spec$language, "/dif_", spec$group_var, "_Rasch_fits.Rds")
  if (!file.exists(fits_cache)) {
    fits <- fit_group_irt(ld$d_mat, ld$d_demo, spec$group_var)
    saveRDS(fits, fits_cache)
  } else {
    fits <- readRDS(fits_cache)
  }

  dif_tab <- compute_dif(fits, spec$group1, spec$group2, ld$d_items)
  dif_results[[key]] <- dif_tab

  cat("\n\nMedian b_diff (", spec$group1, " b - ", spec$group2, " b): ",
      round(median(dif_tab$b_diff), 2),
      " (mean ", round(mean(dif_tab$b_diff), 2), ", sd ", round(sd(dif_tab$b_diff), 2), ") | ",
      sum(dif_tab$b_diff > 0), "/", nrow(dif_tab), " items easier for ", spec$group2, "\n\n", sep = "")

  print(
    ggplot(dif_tab, aes(x = b_diff)) +
      geom_histogram(bins = 40) +
      geom_vline(xintercept = 0, linetype = "dashed") +
      theme_classic() +
      xlab(paste0("Difficulty difference (", spec$group1, " b - ", spec$group2, " b)"))
  )

  print(dif_dotplot(dif_tab, spec$group1, spec$group2))

  cat("\n\n**Items with extreme (>2.25 SD) difficulty difference:**\n\n")
  print(
    dif_tab %>%
      filter(extreme) %>%
      arrange(b_diff) %>%
      select(item_definition, item_kind, domain, b_diff, b_diff_z, favors) %>%
      knitr::kable(digits = 2)
  )

  prune_cache <- paste0("models/", spec$language, "/dif_", spec$group_var, "_Rasch_prune.Rds")
  if (!file.exists(prune_cache)) {
    prune_tab <- prune_and_refit(ld$d_mat, ld$d_demo, dif_tab, spec$group_var)
    saveRDS(prune_tab, prune_cache)
  } else {
    prune_tab <- readRDS(prune_cache)
  }

  print(
    ggplot(prune_tab, aes(x = threshold, y = beta_group)) +
      geom_point() + geom_line() +
      theme_classic() +
      xlab("Exclusion threshold (SD)") +
      ylab(paste0(spec$group_var, " coefficient (beta), ", lang_label[[spec$language]]))
  )
}

English: sex_group

Baseline demographic effect (theta ~ age_c * group):

term estimate std.error statistic p.value
(Intercept) -1.711 0.395 -4.330 0.000
age_c 0.116 0.070 1.663 0.099
sex_groupF 0.308 0.585 0.527 0.599
age_c:sex_groupF 0.075 0.099 0.758 0.450

Median b_diff (M b - F b): 0.68 (mean 0.58, sd 0.86) | 528/676 items easier for F

Items with extreme (>2.25 SD) difficulty difference:

item_definition item_kind domain b_diff b_diff_z favors
myself pronouns other -2.55 -3.64 M
they pronouns other -2.33 -3.37 M
tricycle vehicles noun -2.33 -3.37 M
bump action_words verb -1.84 -2.81 M
radio household noun -1.81 -2.78 M
beside locations other -1.81 -2.78 M
on top of locations other -1.81 -2.78 M
our pronouns other -1.81 -2.78 M
gum food_drink noun -1.52 -2.44 M
could helping_verbs other -1.52 -2.44 M
about locations other -1.52 -2.44 M
circus places other -1.52 -2.44 M
downtown places other -1.52 -2.44 M
their pronouns other -1.52 -2.44 M
these pronouns other -1.52 -2.44 M
yourself pronouns other -1.52 -2.44 M
tonight time_words other -1.43 -2.34 M
coke food_drink noun -1.43 -2.34 M
TV furniture_rooms noun 2.61 2.35 F
owie/ boo boo body_parts noun 2.61 2.36 F
swing action_words verb 2.69 2.45 F
allgone descriptive_words other 2.69 2.45 F
sister people other 2.69 2.45 F
sky outdoor noun 2.79 2.57 F
finish action_words verb 2.88 2.66 F

Vietnamese: sex_group

Baseline demographic effect (theta ~ age_c * group):

term estimate std.error statistic p.value
(Intercept) -0.551 0.241 -2.283 0.024
age_c 0.190 0.043 4.453 0.000
sex_groupF -0.109 0.357 -0.305 0.761
age_c:sex_groupF 0.013 0.060 0.217 0.829

Median b_diff (M b - F b): -0.75 (mean -0.76, sd 0.6) | 62/687 items easier for F

Items with extreme (>2.25 SD) difficulty difference:

item_definition item_kind domain b_diff b_diff_z favors
đài household noun -3.39 -4.37 M
cái của họ pronouns other -2.51 -2.91 M
để tôi helping_verbs other -2.51 -2.91 M
mái nhà outdoor noun -2.39 -2.71 M
cô ấy pronouns other -2.37 -2.67 M
dã ngoại places other -2.37 -2.67 M
găng tay hở ngón clothing noun -2.29 -2.53 M
trước time_words other -2.26 -2.49 M
vòi tưới outdoor noun -2.26 -2.49 M
quần áo mùa đông clothing noun -2.21 -2.40 M
bắp rang bơ food_drink noun -2.12 -2.26 M
ầm descriptive_words other 0.60 2.27 F
cá animals noun 0.63 2.31 F
tóe action_words verb 0.68 2.40 F
váy clothing noun 0.73 2.48 F
ghế dài furniture_rooms noun 0.75 2.51 F
mua sắm games_routines other 0.80 2.59 F
nước có ga food_drink noun 0.92 2.81 F
rết animals noun 1.02 2.96 F

English: edu_group

Baseline demographic effect (theta ~ age_c * group):

term estimate std.error statistic p.value
(Intercept) -1.541 0.691 -2.230 0.027
age_c 0.090 0.121 0.747 0.456
edu_grouphigh -0.006 0.766 -0.007 0.994
age_c:edu_grouphigh 0.070 0.133 0.529 0.598

Median b_diff (low b - high b): 1.51 (mean 1.5, sd 0.9) | 634/667 items easier for high

Items with extreme (>2.25 SD) difficulty difference:

item_definition item_kind domain b_diff b_diff_z favors
person people other -3.18 -5.21 low
on top of locations other -1.24 -3.05 low
porch furniture_rooms noun -1.11 -2.91 low
woods places other -0.93 -2.71 low
sweep action_words verb -0.85 -2.62 low
dump action_words verb -0.85 -2.62 low
his pronouns other -0.85 -2.62 low
farm places other -0.80 -2.56 low
country places other -0.71 -2.46 low
soda/pop food_drink noun -0.65 -2.39 low
every quantifiers other -0.65 -2.39 low
people people other -0.65 -2.39 low
him pronouns other -0.65 -2.39 low
cereal/granola food_drink noun 3.65 2.40 high
peanut butter food_drink noun 3.65 2.40 high
swing action_words verb 4.21 3.02 high
cheese food_drink noun 4.27 3.09 high

Vietnamese: edu_group

Baseline demographic effect (theta ~ age_c * group):

term estimate std.error statistic p.value
(Intercept) -0.850 0.425 -2.002 0.047
age_c 0.105 0.074 1.420 0.158
edu_grouphigh 0.317 0.470 0.674 0.502
age_c:edu_grouphigh 0.118 0.082 1.451 0.149

Median b_diff (low b - high b): 0.4 (mean 0.46, sd 0.72) | 506/683 items easier for high

Items with extreme (>2.25 SD) difficulty difference:

item_definition item_kind domain b_diff b_diff_z favors
gà tây animals noun -1.70 -3.01 low
nhà places other -1.37 -2.56 low
cái của họ pronouns other -1.24 -2.37 low
ghế bành furniture_rooms noun -1.24 -2.37 low
mọi cái quantifiers other -1.24 -2.37 low
không phải quantifiers other -1.22 -2.35 low
rết animals noun -1.21 -2.34 low
nai animals noun -1.21 -2.33 low
phòng furniture_rooms noun -1.18 -2.28 low
bữa tiệc places other -1.17 -2.27 low
be be sounds other -1.16 -2.26 low
đồng xu household noun 2.13 2.31 high
khăn tắm household noun 2.15 2.34 high
giày thể thao clothing noun 2.15 2.35 high
thìa household noun 2.21 2.42 high
bơ đậu phộng food_drink noun 2.25 2.49 high
ầm descriptive_words other 2.25 2.49 high
chậm descriptive_words other 2.32 2.59 high
người tuyết outdoor noun 2.32 2.59 high
cốc household noun 2.35 2.63 high
máy cắt cỏ outdoor noun 2.38 2.66 high
chổi lau sàn household noun 2.81 3.27 high
máy kéo vehicles noun 2.91 3.41 high
băng dính household noun 3.01 3.54 high

English: lang_group

Baseline demographic effect (theta ~ age_c * group):

term estimate std.error statistic p.value
(Intercept) -2.006 0.338 -5.932 0.000
age_c 0.146 0.054 2.715 0.007
lang_groupmore_than_2 1.417 0.627 2.259 0.025
age_c:lang_groupmore_than_2 0.069 0.121 0.568 0.571

Median b_diff (2_or_fewer b - more_than_2 b): 0.53 (mean 0.48, sd 0.85) | 505/677 items easier for more_than_2

Items with extreme (>2.25 SD) difficulty difference:

item_definition item_kind domain b_diff b_diff_z favors
for locations other -2.34 -3.31 2_or_fewer
at locations other -2.34 -3.31 2_or_fewer
church places other -2.34 -3.31 2_or_fewer
the quantifiers other -2.34 -3.31 2_or_fewer
think action_words verb -2.11 -3.04 2_or_fewer
not quantifiers other -1.79 -2.66 2_or_fewer
are helping_verbs other -1.61 -2.45 2_or_fewer
is helping_verbs other -1.61 -2.45 2_or_fewer
tear action_words verb -1.57 -2.40 2_or_fewer
tiny descriptive_words other -1.57 -2.40 2_or_fewer
with locations other -1.57 -2.40 2_or_fewer
him pronouns other -1.57 -2.40 2_or_fewer
yourself pronouns other -1.57 -2.40 2_or_fewer
their pronouns other -1.57 -2.40 2_or_fewer
any quantifiers other -1.57 -2.40 2_or_fewer
same quantifiers other -1.57 -2.40 2_or_fewer
don’t helping_verbs other -1.52 -2.35 2_or_fewer
raisin food_drink noun 2.42 2.27 more_than_2
swing action_words verb 2.60 2.48 more_than_2
this little piggy games_routines other 2.87 2.80 more_than_2
vagina body_parts noun 2.93 2.87 more_than_2
patty cake games_routines other 3.22 3.20 more_than_2
nail household noun 3.52 3.56 more_than_2

Vietnamese: lang_group

Baseline demographic effect (theta ~ age_c * group):

term estimate std.error statistic p.value
(Intercept) -0.742 0.209 -3.558 0.001
age_c 0.195 0.033 5.872 0.000
lang_groupmore_than_2 0.487 0.387 1.257 0.211
age_c:lang_groupmore_than_2 0.034 0.075 0.453 0.651

Median b_diff (2_or_fewer b - more_than_2 b): 1.13 (mean 1.05, sd 0.75) | 634/687 items easier for more_than_2

Items with extreme (>2.25 SD) difficulty difference:

item_definition item_kind domain b_diff b_diff_z favors
nghèo descriptive_words other -2.46 -4.66 2_or_fewer
tóe action_words verb -1.78 -3.75 2_or_fewer
trống rỗng descriptive_words other -1.78 -3.75 2_or_fewer
kinh tởm descriptive_words other -1.61 -3.53 2_or_fewer
thời gian time_words other -1.61 -3.53 2_or_fewer
sân cát outdoor noun -1.43 -3.29 2_or_fewer
cũng thế quantifiers other -1.43 -3.29 2_or_fewer
phà outdoor noun -1.43 -3.29 2_or_fewer
chúng tôi pronouns other -1.43 -3.29 2_or_fewer
một vài quantifiers other -1.43 -3.29 2_or_fewer
khoai tây chiên food_drink noun -1.37 -3.21 2_or_fewer
giấu action_words verb -1.29 -3.09 2_or_fewer
thế và connecting_words other -1.24 -3.03 2_or_fewer
ghế bành furniture_rooms noun -1.24 -3.03 2_or_fewer
gạt action_words verb -1.09 -2.84 2_or_fewer
vợt toys noun -1.09 -2.84 2_or_fewer
cách xa locations other -0.97 -2.68 2_or_fewer
hạt trang trí clothing noun -0.77 -2.41 2_or_fewer
bắt cứ cái gì quantifiers other -0.77 -2.41 2_or_fewer
nếm action_words verb -0.74 -2.37 2_or_fewer
chú hề people other -0.71 -2.33 2_or_fewer
người đưa thư people other -0.71 -2.33 2_or_fewer
của họ pronouns other -0.71 -2.33 2_or_fewer
ngay cạnh locations other -0.71 -2.33 2_or_fewer
bánh quy giòn food_drink noun 3.34 3.04 more_than_2

Cross-language comparison of DIF patterns

Once the loop above has run, dif_results holds one item-level DIF table per (language, demographic) combination (keyed "en_sex_group", "vi_sex_group", etc.). A natural follow-up, paralleling Part 1’s category-level cross-language check: are the same categories of items biased for the same demographic groups in both languages (e.g. are stereotypically-gendered nouns sex-biased in both English and Vietnamese, as the paper found for English alone)?

summarize_dif_by_category <- function(dif_tab, language, group_var) {
  dif_tab %>%
    group_by(item_kind, domain) %>%
    summarise(n_items = n(), mean_b_diff = mean(b_diff), pct_extreme = mean(extreme), .groups = "drop") %>%
    mutate(language = language, group_var = group_var)
}

dif_by_category <- bind_rows(
  summarize_dif_by_category(dif_results[["en_sex_group"]],  "English",    "sex_group"),
  summarize_dif_by_category(dif_results[["vi_sex_group"]],  "Vietnamese", "sex_group"),
  summarize_dif_by_category(dif_results[["en_edu_group"]],  "English",    "edu_group"),
  summarize_dif_by_category(dif_results[["vi_edu_group"]],  "Vietnamese", "edu_group"),
  summarize_dif_by_category(dif_results[["en_lang_group"]], "English",    "lang_group"),
  summarize_dif_by_category(dif_results[["vi_lang_group"]], "Vietnamese", "lang_group")
)

dif_by_category %>%
  kable(digits = 2, caption = "Mean item-level DIF and % extreme items by category, language, and demographic.") %>%
  html_table_width(c(120, 60, 60, 90, 80, 90, 90))
Mean item-level DIF and % extreme items by category, language, and demographic.
item_kind domain n_items mean_b_diff pct_extreme language group_var
action_words verb 103 0.81 0.03 English sex_group
animals noun 43 0.70 0.00 English sex_group
body_parts noun 26 1.05 0.04 English sex_group
clothing noun 27 0.33 0.00 English sex_group
connecting_words other 6 0.64 0.00 English sex_group
descriptive_words other 63 0.75 0.02 English sex_group
food_drink noun 68 0.60 0.03 English sex_group
furniture_rooms noun 34 0.74 0.03 English sex_group
games_routines other 25 1.15 0.00 English sex_group
helping_verbs other 20 0.41 0.05 English sex_group
household noun 50 0.64 0.02 English sex_group
locations other 26 0.16 0.12 English sex_group
outdoor noun 31 0.57 0.03 English sex_group
people other 28 0.59 0.04 English sex_group
places other 21 -0.12 0.10 English sex_group
pronouns other 25 -0.03 0.24 English sex_group
quantifiers other 17 0.61 0.00 English sex_group
question_words other 7 0.87 0.00 English sex_group
sounds other 12 0.40 0.00 English sex_group
time_words other 12 -0.45 0.08 English sex_group
toys noun 18 0.57 0.00 English sex_group
vehicles noun 14 -0.23 0.07 English sex_group
action_words verb 103 -0.71 0.01 Vietnamese sex_group
animals noun 46 -0.51 0.04 Vietnamese sex_group
body_parts noun 29 -0.73 0.00 Vietnamese sex_group
clothing noun 28 -0.66 0.11 Vietnamese sex_group
connecting_words other 6 -0.61 0.00 Vietnamese sex_group
descriptive_words other 62 -0.98 0.02 Vietnamese sex_group
food_drink noun 69 -0.69 0.03 Vietnamese sex_group
furniture_rooms noun 33 -0.70 0.03 Vietnamese sex_group
games_routines other 25 -0.70 0.04 Vietnamese sex_group
helping_verbs other 12 -0.57 0.08 Vietnamese sex_group
household noun 58 -0.77 0.02 Vietnamese sex_group
locations other 24 -0.83 0.00 Vietnamese sex_group
outdoor noun 35 -0.96 0.06 Vietnamese sex_group
people other 29 -0.60 0.00 Vietnamese sex_group
places other 23 -0.91 0.04 Vietnamese sex_group
pronouns other 24 -1.05 0.08 Vietnamese sex_group
quantifiers other 16 -0.70 0.00 Vietnamese sex_group
question_words other 7 -0.55 0.00 Vietnamese sex_group
sounds other 12 -0.62 0.00 Vietnamese sex_group
time_words other 12 -1.37 0.08 Vietnamese sex_group
toys noun 19 -0.52 0.00 Vietnamese sex_group
vehicles noun 15 -1.21 0.00 Vietnamese sex_group
action_words verb 103 1.27 0.03 English edu_group
animals noun 42 2.02 0.00 English edu_group
body_parts noun 27 1.87 0.00 English edu_group
clothing noun 26 1.77 0.00 English edu_group
connecting_words other 6 1.21 0.00 English edu_group
descriptive_words other 63 1.42 0.00 English edu_group
food_drink noun 67 2.03 0.06 English edu_group
furniture_rooms noun 34 1.36 0.03 English edu_group
games_routines other 25 1.61 0.00 English edu_group
helping_verbs other 16 0.99 0.00 English edu_group
household noun 48 1.73 0.00 English edu_group
locations other 25 0.74 0.04 English edu_group
outdoor noun 31 1.43 0.00 English edu_group
people other 27 1.14 0.07 English edu_group
places other 22 0.82 0.14 English edu_group
pronouns other 25 0.92 0.08 English edu_group
quantifiers other 17 1.37 0.06 English edu_group
question_words other 7 1.56 0.00 English edu_group
sounds other 12 2.20 0.00 English edu_group
time_words other 12 1.40 0.00 English edu_group
toys noun 18 1.88 0.00 English edu_group
vehicles noun 14 1.59 0.00 English edu_group
action_words verb 103 0.56 0.00 Vietnamese edu_group
animals noun 46 0.31 0.07 Vietnamese edu_group
body_parts noun 29 0.83 0.00 Vietnamese edu_group
clothing noun 27 0.62 0.04 Vietnamese edu_group
connecting_words other 6 -0.02 0.00 Vietnamese edu_group
descriptive_words other 62 0.48 0.03 Vietnamese edu_group
food_drink noun 67 0.51 0.01 Vietnamese edu_group
furniture_rooms noun 33 0.42 0.06 Vietnamese edu_group
games_routines other 25 0.30 0.00 Vietnamese edu_group
helping_verbs other 12 -0.04 0.00 Vietnamese edu_group
household noun 58 0.84 0.10 Vietnamese edu_group
locations other 24 0.53 0.00 Vietnamese edu_group
outdoor noun 35 0.77 0.06 Vietnamese edu_group
people other 29 0.20 0.00 Vietnamese edu_group
places other 22 0.14 0.09 Vietnamese edu_group
pronouns other 24 -0.14 0.04 Vietnamese edu_group
quantifiers other 16 0.17 0.12 Vietnamese edu_group
question_words other 7 -0.14 0.00 Vietnamese edu_group
sounds other 12 0.74 0.08 Vietnamese edu_group
time_words other 12 0.28 0.00 Vietnamese edu_group
toys noun 19 0.02 0.00 Vietnamese edu_group
vehicles noun 15 0.78 0.07 Vietnamese edu_group
action_words verb 103 0.29 0.03 English lang_group
animals noun 43 0.60 0.00 English lang_group
body_parts noun 27 0.82 0.04 English lang_group
clothing noun 27 0.80 0.00 English lang_group
connecting_words other 6 0.49 0.00 English lang_group
descriptive_words other 63 0.19 0.02 English lang_group
food_drink noun 68 0.92 0.01 English lang_group
furniture_rooms noun 33 0.87 0.00 English lang_group
games_routines other 25 1.13 0.08 English lang_group
helping_verbs other 20 0.02 0.15 English lang_group
household noun 50 0.63 0.02 English lang_group
locations other 25 -0.15 0.12 English lang_group
outdoor noun 31 0.51 0.00 English lang_group
people other 29 0.42 0.00 English lang_group
places other 22 0.42 0.05 English lang_group
pronouns other 25 -0.40 0.12 English lang_group
quantifiers other 17 -0.44 0.24 English lang_group
question_words other 7 0.61 0.00 English lang_group
sounds other 12 1.19 0.00 English lang_group
time_words other 12 0.76 0.00 English lang_group
toys noun 18 0.55 0.00 English lang_group
vehicles noun 14 0.52 0.00 English lang_group
action_words verb 103 0.80 0.04 Vietnamese lang_group
animals noun 46 1.40 0.00 Vietnamese lang_group
body_parts noun 29 0.97 0.00 Vietnamese lang_group
clothing noun 28 1.22 0.04 Vietnamese lang_group
connecting_words other 6 0.38 0.17 Vietnamese lang_group
descriptive_words other 62 0.81 0.05 Vietnamese lang_group
food_drink noun 69 1.41 0.03 Vietnamese lang_group
furniture_rooms noun 33 1.21 0.03 Vietnamese lang_group
games_routines other 25 1.27 0.00 Vietnamese lang_group
helping_verbs other 12 1.03 0.00 Vietnamese lang_group
household noun 58 1.24 0.00 Vietnamese lang_group
locations other 24 0.99 0.08 Vietnamese lang_group
outdoor noun 35 0.96 0.06 Vietnamese lang_group
people other 29 0.77 0.07 Vietnamese lang_group
places other 23 1.12 0.00 Vietnamese lang_group
pronouns other 24 0.53 0.08 Vietnamese lang_group
quantifiers other 16 0.75 0.19 Vietnamese lang_group
question_words other 7 1.29 0.00 Vietnamese lang_group
sounds other 12 1.39 0.00 Vietnamese lang_group
time_words other 12 1.02 0.08 Vietnamese lang_group
toys noun 19 0.94 0.05 Vietnamese lang_group
vehicles noun 15 1.24 0.00 Vietnamese lang_group

Same visual logic as Part 1’s category-level EN-VI comparison (part1-category-crosslang-plot): each point is one item_kind, with its mean b_diff (favors group1 when negative, group2 when positive – see dif_specs above for which demographic level is group1/group2 in each panel) plotted for English (x) against Vietnamese (y). Categories in the upper-right or lower-left quadrants are biased toward the same group in both languages; categories in the upper-left/lower-right are biased toward opposite groups depending on language.

dif_by_category_wide <- dif_by_category %>%
  select(group_var, item_kind, domain, language, mean_b_diff) %>%
  pivot_wider(names_from = language, values_from = mean_b_diff) %>%
  filter(!is.na(English), !is.na(Vietnamese))

ggplot(dif_by_category_wide, aes(x = English, y = Vietnamese, color = domain, label = item_kind)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey70") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey70") +
  geom_point() +
  ggrepel::geom_text_repel(size = 2.8, show.legend = FALSE, max.overlaps = 15) +
  facet_wrap(~group_var, nrow = 1) +
  theme_classic() +
  xlab("English: mean b_diff (group1 - group2)") +
  ylab("Vietnamese: mean b_diff (group1 - group2)")

Results: Most words show the same demographic

Future extensions

  • GLIMMER-style uncertainty bands: refit fit_group_irt() with SE = TRUE and draw imputations from each item’s parameter covariance matrix (see Stenhaug, Frank, & Domingue 2021, cited in Kachergis et al. 2022) to distinguish genuine DIF clusters from noise, rather than relying on point estimates alone.
  • Race/ethnicity, birth order: not available in this dataset’s demographics as loaded here, but the same fit_group_irt()/compute_dif()/prune_and_refit() machinery generalizes to any binary (or map-able multi-level) grouping variable added to d_demo.
  • Child-directed speech frequency: the paper’s final analysis relates sex-based DIF to word frequency in child-directed speech by target-child sex (CHILDES). No equivalent Vietnamese corpus is used here; if one becomes available, the same log-odds-of-frequency-by-group approach could be correlated against b_diff.

Part 3: Cross-language item-level difficulty comparison

Are item difficulties consistent in English and Vietnamese – i.e. do words that are hard for English-Vietnamese bilingual children to produce in English also tend to be hard for them to produce in Vietnamese? This extracts each item’s IRT difficulty (b) from the same mod_2pl_en/mod_2pl_vi models used in Part 1, z-scores it within each language (each item’s b relative to that language’s own full item-set mean/SD – the same within-language standardization used for category_wide in Part 1 and for ability_noun in 03_bifactor_noun_verb.Rmd, for the same reason: the two models are independently calibrated, so raw b isn’t on a shared scale), and correlates the standardized difficulties (Pearson) for translation-equivalent item pairs (e.g. “dog” / “con chó”) – unlike Parts 1-2, this needs an item-level (not just category-level) EN-VI mapping.

A translation-equivalence crosswalk is required and does not yet exist in this repo

Matching items across languages needs an item-level EN-VI translation crosswalk (which English item_definition corresponds to which Vietnamese item_definition). The two fields files share the same item_kind categories (used in Parts 1-2 above) but list words in each language’s own order (Vietnamese items appear alphabetized by the Vietnamese word), so there’s no positional correspondence to exploit – the mapping has to come from the original CDI adaptation/translation documentation, or be built by hand.

The chunk below looks for data/en_vi_item_crosswalk.csv (columns: item_kind, en_definition, vi_definition). If it isn’t there, it writes out a template (one row per English word-production item, vi_definition left blank) for someone to fill in by hand or from the adaptation documentation, and the rest of this Part is skipped until that file exists.

crosswalk_path <- "data/en_vi_item_crosswalk.csv"

if (!file.exists(crosswalk_path)) {
  template <- en_data$d_items %>%
    filter(item_kind %in% intersect(en_data$d_items$item_kind, vi_data$d_items$item_kind)) %>%
    transmute(item_kind, en_definition = item_definition, vi_definition = NA_character_)
  write_csv(template, crosswalk_path)
  warning(
    "No translation crosswalk found. Wrote a template with ", nrow(template),
    " English items (vi_definition left blank) to '", crosswalk_path, "' -- ",
    "fill in the Vietnamese translation equivalent for each row (or source it from the ",
    "CDI adaptation documentation) before re-running this analysis."
  )
}
crosswalk_ready <- file.exists(crosswalk_path) && !all(is.na(read_csv(crosswalk_path, show_col_types = FALSE)$vi_definition))

The remaining chunks in this Part are set to only evaluate once data/en_vi_item_crosswalk.csv exists and is filled in (eval=crosswalk_ready); until then they’re a no-op.

crosswalk <- read_csv(crosswalk_path, show_col_types = FALSE) %>% filter(!is.na(vi_definition))

# z-scored against each language's FULL item set (diff_en/diff_vi, all
# ~680/~690 word-production items), not just the matched subset below --
# same convention as category_wide in Part 1 and cross_lang_join in
# 03_bifactor_noun_verb.Rmd, so a matched pair's z-scores reflect where
# each word sits in its own language's overall difficulty distribution.
en_mean_b <- mean(diff_en$b); en_sd_b <- sd(diff_en$b)
vi_mean_b <- mean(diff_vi$b); vi_sd_b <- sd(diff_vi$b)

matched <- crosswalk %>%
  inner_join(diff_en %>% rename(en_b = b, en_item_kind = item_kind) %>% select(-difficulty_rank),
             by = c("en_definition" = "item_definition")) %>%
  inner_join(diff_vi %>% rename(vi_b = b, vi_item_kind = item_kind) %>% select(-difficulty_rank),
             by = c("vi_definition" = "item_definition")) %>%
  mutate(
    en_b_z = (en_b - en_mean_b) / en_sd_b,
    vi_b_z = (vi_b - vi_mean_b) / vi_sd_b
  )
## Warning in inner_join(., diff_en %>% rename(en_b = b, en_item_kind = item_kind) %>% : Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 11 of `x` matches multiple rows in `y`.
## ℹ Row 262 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
## Warning in inner_join(., diff_vi %>% rename(vi_b = b, vi_item_kind = item_kind) %>% : Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 36 of `x` matches multiple rows in `y`.
## ℹ Row 324 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
cat(nrow(matched), "translation-equivalent item pairs matched",
    "(out of", nrow(diff_en), "English and", nrow(diff_vi), "Vietnamese word-production items).")
## 728 translation-equivalent item pairs matched (out of 681 English and 687 Vietnamese word-production items).

Overall correlation (standardized difficulty)

cor.test(matched$en_b_z, matched$vi_b_z, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  matched$en_b_z and matched$vi_b_z
## t = 19.912, df = 726, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5452106 0.6393792
## sample estimates:
##       cor 
## 0.5943283
ggplot(matched, aes(x = en_b_z, y = vi_b_z)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey70") +
  theme_classic() +
  xlab("English difficulty (z-scored within language)") +
  ylab("Vietnamese difficulty (z-scored within language)")
## `geom_smooth()` using formula = 'y ~ x'

Is the correlation mediated by item category (noun/verb/other)?

matched <- matched %>%
  mutate(domain = case_when(
    en_item_kind %in% verb_kinds ~ "verb",
    en_item_kind %in% noun_kinds ~ "noun",
    TRUE ~ "other"
  ))

matched %>%
  group_by(domain) %>%
  summarise(
    n_items = n(),
    pearson_r = cor(en_b_z, vi_b_z, method = "pearson"),
    .groups = "drop"
  ) %>%
  kable(digits = 2, caption = "Item-difficulty correlation (standardized b) between English and Vietnamese, by domain.") %>%
  html_table_width(c(100, 90, 100))
Item-difficulty correlation (standardized b) between English and Vietnamese, by domain.
domain n_items pearson_r
noun 342 0.64
other 270 0.52
verb 116 0.54 
ggplot(matched, aes(x = en_b_z, y = vi_b_z, color = domain)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = FALSE) +
  facet_wrap(~domain) +
  theme_classic() +
  theme(legend.position = "none") +
  xlab("English difficulty (z-scored within language)") +
  ylab("Vietnamese difficulty (z-scored within language)")
## `geom_smooth()` using formula = 'y ~ x'

Future extensions

  • Concreteness / socialness / other item-level norms: not currently in this repo; would need an external norms dataset (e.g. Brysbaert et al. concreteness norms) joined onto matched by en_definition, then tested as a moderator of the correlation (e.g. lm(vi_b_z ~ en_b_z * concreteness, data = matched)).
  • Age-of-acquisition (AoA) norms: AoA is hard to estimate directly for bilingual children (interacts with age of exposure onset), which is part of why Part 1 above uses item-difficulty rank / production-rate trajectory as a cleaner developmental-order proxy. If AoA norms for both languages become available, the natural check is the same approach above, run between item difficulty (z-scored) and AoA within each language.