Motivation

Among the 142 children with both an English and a Vietnamese CDI, data/participant_info.csv provides a researcher-coded “bi”/“mono” classification per child (group column) – this is the authoritative grouping used throughout this analysis, rather than a production-based proxy. An earlier version of this analysis defined “Monolingual” as “produced zero English words on this survey” (42/142 children), since the actual group label wasn’t yet joined in; that production-based proxy is kept below purely as a cross-check against the authoritative label, and any disagreements between the two are reported explicitly rather than papered over.

Question: are there specific Vietnamese words that are differentially likely to be produced by one group vs. the other – i.e. words that are “representative” of (more diagnostic of) monolingual vs. bilingual production, beyond what’s just explained by chance or by age?

Caveat on what “monolingual” means here: we don’t have participant_info.csv’s exact coding criterion in hand (e.g. whether “mono” reflects parent-reported language exposure, household language, or something else) – it’s used here as the intended ground-truth group label precisely because it should capture more than this one survey’s production snapshot, but the precise operational definition behind it should be confirmed against the original study documentation before this label is reported in a paper. We check below whether this grouping is confounded with age (it isn’t, much) and control for age in a secondary check regardless.

Setup: load data and define groups

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_")
)

# Mirrors load_lang_data() in 03_bifactor_noun_verb.Rmd / 05_item_difficulty_crosslang.Rmd:
# d_mat's column names are implicitly sanitized to match d_items$definition via
# data.frame()'s default check.names = TRUE (e.g. "french fries" -> "french.fries"),
# which is why we build d_mat this way rather than pivoting the raw tibble directly
# (see 05_item_difficulty_crosslang.Rmd's compute_item_aoa()/compute_category_trajectory()
# fix for what goes wrong if you don't).
load_lang_data <- function(lang) {
  cfg <- lang_config[[lang]]
  raw <- read_csv(cfg$data_file, show_col_types = FALSE)

  d_demo <- raw %>% select(row_id, response_id, age)
  d_wide <- raw %>%
    select(response_id, starts_with(cfg$prefix)) %>%
    mutate(across(everything(), ~replace_na(.x, 0)))
  d_mat <- d_wide %>% data.frame() %>% select(-response_id) %>% data.matrix()

  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, "^.*?___")
    ) %>%
    select(definition, item_kind, item_definition)
  stopifnot(all(colnames(d_mat) == d_items$definition))

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

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

row_id is the child-level identifier shared by the English and Vietnamese surveys (established in 03_bifactor_noun_verb.Rmd), confirmed below to fully overlap (142 children in both) with identical age across surveys, so we can attach each child’s English production status onto their Vietnamese responses positionally. participant_info.csv’s own child identifier is ResponseId, which we verified separately matches response_id in both the English and Vietnamese data files for all 142 children (not row_id, which doesn’t appear in participant_info.csv at all) – so that join uses response_id instead.

production_en <- tibble(row_id = en_data$d_demo$row_id, production_en = rowSums(en_data$d_mat))

participant_info <- read_csv("data/participant_info.csv", show_col_types = FALSE) %>%
  transmute(
    response_id = ResponseId,
    group_reported = factor(recode(group, bi = "Bilingual", mono = "Monolingual"),
                             levels = c("Bilingual", "Monolingual"))
  )
## New names:
## • `` -> `...1`
vi_demo <- vi_data$d_demo %>%
  left_join(production_en, by = "row_id") %>%
  left_join(participant_info, by = "response_id")

stopifnot(!anyNA(vi_demo$production_en))     # every Vietnamese respondent also has an English response
stopifnot(!anyNA(vi_demo$group_reported))    # every Vietnamese respondent also has a participant_info row
stopifnot(nrow(vi_demo) == nrow(vi_data$d_demo))  # neither left_join duplicated/dropped any rows

vi_demo <- vi_demo %>%
  mutate(
    group_production = factor(if_else(production_en == 0, "Monolingual", "Bilingual"),
                               levels = c("Bilingual", "Monolingual")),
    group = group_reported,  # authoritative grouping used throughout the rest of this analysis
    age_c = age - mean(age, na.rm = TRUE)
  )

table(vi_demo$group)
## 
##   Bilingual Monolingual 
##         101          41

Does the authoritative grouping agree with the production-based proxy?

table(Reported = vi_demo$group_reported, `Production-based` = vi_demo$group_production)
##              Production-based
## Reported      Bilingual Monolingual
##   Bilingual         100           1
##   Monolingual         0          41
contradictions <- vi_demo %>%
  filter(group_reported != group_production) %>%
  select(row_id, age, production_en, group_reported, group_production)

cat(nrow(contradictions), "of", nrow(vi_demo), "children disagree between the two definitions.\n")
## 1 of 142 children disagree between the two definitions.
contradictions %>%
  kable(digits = 1, caption = "Children where the reported group and the production-based proxy disagree.") %>%
  html_table_width(c(70, 60, 90, 100, 110))
Children where the reported group and the production-based proxy disagree.
row_id age production_en group_reported group_production
72 26 0 Bilingual Monolingual

Only one child (of 142) disagrees between the two definitions: reported as Bilingual, but producing zero English words on this particular survey. This is the single case the earlier caveat anticipated – a child whose reported group reflects something broader than this one administration’s production snapshot (e.g. exposure, or ordinary developmental timing) can land on zero English words here without contradicting that broader label. With agreement at 141/142, the reported and production-based definitions are overwhelmingly consistent; the switch to participant_info.csv changes the Monolingual group from 42 to 41 children.

Checking for an age confound: the two groups turn out to be closely age-matched (means ~24.6 vs. ~25.0 months, similar spread), so raw group differences below are unlikely to just be proxying for age – but a secondary, age-adjusted check is still reported alongside for each word.

vi_demo %>%
  group_by(group) %>%
  summarise(n = n(), mean_age = mean(age), sd_age = sd(age), min_age = min(age), max_age = max(age), .groups = "drop") %>%
  kable(digits = 1, caption = "Age (months) by group.") %>%
  html_table_width(rep(90, 5))
Age (months) by group.
group n mean_age sd_age min_age max_age
Bilingual 101 24.6 6.1 14 36
Monolingual 41 25.0 5.7 15 35 
ggplot(vi_demo, aes(x = group, y = age)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.15, alpha = 0.5) +
  theme_classic() +
  xlab(NULL) + ylab("Age (months)")

Item-level group comparison

For each Vietnamese word, we compare production rate between the two groups two ways:

  1. Primary: a 2x2 (group x produces) contingency table with a Haldane-Anscombe continuity correction (+0.5 to every cell) to get a stable log-odds ratio even for words at or near 0%/100% in one group (plain MLE log-odds would be +-Inf for such words, which are exactly the most “distinctive” ones we want to be able to rank) – and a Fisher’s exact test p-value, FDR-corrected (p.adjust(method = "BH")) across all ~680 items tested. log_or > 0 means the word is more likely to be produced by Monolingual (Vietnamese-only) children; log_or < 0 means more likely for Bilingual children.
  2. Secondary, age-adjusted robustness check: glm(produces ~ group + age_c, family = binomial) per word, reported alongside so a word whose raw association disappears (or flips sign) once age is controlled for can be flagged as likely age-driven rather than genuinely group-distinctive. This can fail to converge or blow up (quasi- complete separation) for rare words even with age included – those are set to NA rather than trusted, the same convention used for non-convergent per-item fits elsewhere in this project (e.g. compute_item_aoa() in 05_item_difficulty_crosslang.Rmd).
compare_mono_bilingual <- function(d_mat, group, age_c, d_items) {
  n_items <- ncol(d_mat)
  defs <- colnames(d_mat)

  map_dfr(seq_len(n_items), function(j) {
    y <- d_mat[, j]

    tab <- table(group, y)
    if (!("0" %in% colnames(tab))) tab <- cbind(`0` = c(0, 0), tab)
    if (!("1" %in% colnames(tab))) tab <- cbind(tab, `1` = c(0, 0))
    tab <- tab[, c("0", "1")]

    a <- tab["Bilingual", "1"] + 0.5;    b <- tab["Bilingual", "0"] + 0.5
    cc <- tab["Monolingual", "1"] + 0.5; d <- tab["Monolingual", "0"] + 0.5
    log_or    <- log((cc / d) / (a / b))
    se_log_or <- sqrt(1 / a + 1 / b + 1 / cc + 1 / d)

    fisher_p <- tryCatch(fisher.test(tab)$p.value, error = function(e) NA_real_)

    n_bi <- sum(group == "Bilingual"); n_mono <- sum(group == "Monolingual")
    p_bi <- mean(y[group == "Bilingual"]); p_mono <- mean(y[group == "Monolingual"])

    glm_or <- NA_real_; glm_se <- NA_real_
    if (length(unique(y)) > 1) {
      fit <- tryCatch(glm(y ~ group + age_c, family = binomial), error = function(e) NULL)
      if (!is.null(fit) && fit$converged) {
        cf <- summary(fit)$coefficients
        if ("groupMonolingual" %in% rownames(cf)) {
          est <- cf["groupMonolingual", "Estimate"]; se <- cf["groupMonolingual", "Std. Error"]
          if (se <= 5) { glm_or <- est; glm_se <- se }  # else: likely quasi-separation, leave NA
        }
      }
    }

    tibble(definition = defs[j], n_bilingual = n_bi, n_monolingual = n_mono,
           p_bilingual = p_bi, p_monolingual = p_mono, diff = p_mono - p_bi,
           log_or = log_or, se_log_or = se_log_or, fisher_p = fisher_p,
           glm_or_age_adj = glm_or, glm_se_age_adj = glm_se)
  }) %>%
    mutate(p_fdr = p.adjust(fisher_p, method = "BH")) %>%
    left_join(d_items, by = "definition")
}

vi_word_group_diff <- compare_mono_bilingual(vi_data$d_mat, vi_demo$group, vi_demo$age_c, vi_data$d_items)

Words representative of Vietnamese-only (Monolingual) producers

Top 20 words by log-odds ratio (Haldane-Anscombe corrected): most over-represented among children who produce zero English words, relative to bilingual producers.

vi_word_group_diff %>%
  arrange(desc(log_or)) %>%
  slice_head(n = 20) %>%
  select(item_definition, item_kind, n_bilingual, n_monolingual, p_bilingual, p_monolingual,
         log_or, fisher_p, p_fdr, glm_or_age_adj) %>%
  kable(digits = c(0, 0, 0, 0, 2, 2, 2, 3, 3, 2),
        caption = "Top 20 words over-represented among Vietnamese-only (Monolingual) producers.") %>%
  html_table_width(c(140, 80, 60, 60, 70, 70, 60, 60, 60, 90))
Top 20 words over-represented among Vietnamese-only (Monolingual) producers.
item_definition item_kind n_bilingual n_monolingual p_bilingual p_monolingual log_or fisher_p p_fdr glm_or_age_adj
sân cát outdoor 101 41 0.02 0.17 2.16 0.002 0.064 2.55
chúng tôi pronouns 101 41 0.02 0.17 2.16 0.002 0.064 2.55
cách xa locations 101 41 0.05 0.27 1.89 0.001 0.047 1.98
của anh ấy pronouns 101 41 0.07 0.29 1.68 0.001 0.047 1.80
nệm household 101 41 0.14 0.46 1.65 0.000 0.047 1.94
động vật animals 101 41 0.03 0.15 1.64 0.017 0.117 1.84
cũng thế quantifiers 101 41 0.03 0.15 1.64 0.017 0.117 1.97
xẻng outdoor 101 41 0.06 0.24 1.59 0.003 0.072 1.67
hạt trang trí clothing 101 41 0.02 0.10 1.56 0.058 0.203 1.67
tóe action_words 101 41 0.04 0.17 1.55 0.014 0.110 1.61
hoa outdoor 101 41 0.43 0.78 1.53 0.000 0.047 1.69
của chúng ta pronouns 101 41 0.05 0.20 1.49 0.011 0.104 1.60
tất đầu gối clothing 101 41 0.03 0.12 1.44 0.045 0.175 1.55
thế và connecting_words 101 41 0.03 0.12 1.44 0.045 0.175 1.54
bà people 101 41 0.79 0.95 1.44 0.023 0.135 1.63
tạm biệt games_routines 101 41 0.39 0.73 1.43 0.000 0.047 1.54
hươu animals 101 41 0.16 0.44 1.41 0.001 0.047 1.54
gà trống animals 101 41 0.18 0.46 1.36 0.001 0.047 1.45
cổ body_parts 101 41 0.38 0.71 1.36 0.000 0.047 1.50
gấu animals 101 41 0.41 0.73 1.35 0.001 0.047 1.53

Words representative of Bilingual producers

Top 20 words by log-odds ratio, in the opposite direction: most over-represented among children who also produce English words.

vi_word_group_diff %>%
  arrange(log_or) %>%
  slice_head(n = 20) %>%
  select(item_definition, item_kind, n_bilingual, n_monolingual, p_bilingual, p_monolingual,
         log_or, fisher_p, p_fdr, glm_or_age_adj) %>%
  kable(digits = c(0, 0, 0, 0, 2, 2, 2, 3, 3, 2),
        caption = "Top 20 words over-represented among Bilingual producers.") %>%
  html_table_width(c(140, 80, 60, 60, 70, 70, 60, 60, 60, 90))
Top 20 words over-represented among Bilingual producers.
item_definition item_kind n_bilingual n_monolingual p_bilingual p_monolingual log_or fisher_p p_fdr glm_or_age_adj
bắp rang bơ food_drink 101 41 0.15 0.05 -1.04 0.152 0.313 -1.27
quần áo mùa đông clothing 101 41 0.11 0.05 -0.70 0.348 0.513 -0.91
thạch Jello food_drink 101 41 0.14 0.07 -0.60 0.395 0.557 -0.75
pizza food_drink 101 41 0.47 0.34 -0.50 0.195 0.365 -0.65
đậu Hà Lan food_drink 101 41 0.09 0.05 -0.48 0.511 0.653 -0.67
quần đùi clothing 101 41 0.12 0.07 -0.43 0.554 0.680 -0.56
áo ngủ clothing 101 41 0.24 0.17 -0.37 0.502 0.648 -0.48
bánh mì gối food_drink 101 41 0.08 0.05 -0.36 0.724 0.830 -0.53
be be sounds 101 41 0.48 0.39 -0.34 0.457 0.610 -0.43
người đưa thư people 101 41 0.11 0.07 -0.33 0.757 0.853 -0.49
khối toys 101 41 0.11 0.07 -0.33 0.757 0.853 -0.48
phố places 101 41 0.14 0.10 -0.32 0.589 0.713 -0.41
bây giờ time_words 101 41 0.23 0.17 -0.32 0.505 0.649 -0.41
ngày time_words 101 41 0.17 0.12 -0.32 0.613 0.737 -0.41
xe vehicles 101 41 0.72 0.66 -0.31 0.543 0.673 -0.33
quạc quạc sounds 101 41 0.74 0.68 -0.30 0.535 0.668 -0.41
meo sounds 101 41 0.82 0.78 -0.28 0.638 0.757 -0.35
một cái quantifiers 101 41 0.25 0.20 -0.27 0.662 0.777 -0.34
áo len clothing 101 41 0.22 0.17 -0.26 0.648 0.765 -0.35
sân sau outdoor 101 41 0.13 0.10 -0.24 0.778 0.870 -0.33

Volcano plot: effect size vs. significance

ggplot(vi_word_group_diff, aes(x = log_or, y = -log10(fisher_p), color = p_fdr < 0.05)) +
  geom_point(alpha = 0.6) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
  theme_classic() +
  scale_color_manual(values = c(`TRUE` = "#D55E00", `FALSE` = "grey70")) +
  xlab("Log-odds ratio (Monolingual vs. Bilingual, Haldane-Anscombe corrected)") +
  ylab(expression(-log[10]("Fisher's exact p-value"))) +
  labs(color = "FDR < .05",
       title = "Words distinguishing Vietnamese-only vs. bilingual producers")

How many words survive FDR correction at all, in each direction:

vi_word_group_diff %>%
  filter(p_fdr < 0.05) %>%
  mutate(direction = if_else(log_or > 0, "More likely: Monolingual", "More likely: Bilingual")) %>%
  count(direction) %>%
  kable(caption = "Number of words significant at FDR < .05, by direction.") %>%
  html_table_width(c(200, 80))
Number of words significant at FDR < .05, by direction.
direction n
More likely: Monolingual 13

Since all 13 FDR-significant words favor the Monolingual direction, the raw production-rate gap for each is shown directly below (a dumbbell plot: Bilingual rate vs. Monolingual rate per word, ordered by gap size) – this is the same 13 words as the top rows of the “words representative of Vietnamese-only producers” table above, but only the FDR-significant subset, and as a plot rather than a table.

sig_words_plot_df <- vi_word_group_diff %>%
  filter(p_fdr < 0.05, log_or > 0) %>%
  mutate(item_definition = fct_reorder(item_definition, diff))

ggplot(sig_words_plot_df, aes(y = item_definition)) +
  geom_segment(aes(x = p_bilingual, xend = p_monolingual, y = item_definition, yend = item_definition),
               color = "grey70", linewidth = 0.6) +
  geom_point(aes(x = p_bilingual, color = "Bilingual"), size = 2.5) +
  geom_point(aes(x = p_monolingual, color = "Monolingual"), size = 2.5) +
  theme_classic() +
  xlab("Proportion of children producing") + ylab(NULL) +
  labs(color = NULL, title = "Production rate: Bilingual vs. Vietnamese-only, FDR-significant words")

Are distinguishing words concentrated in particular categories?

vi_word_group_diff %>%
  filter(p_fdr < 0.05) %>%
  mutate(direction = if_else(log_or > 0, "More likely: Monolingual", "More likely: Bilingual")) %>%
  count(item_kind, direction) %>%
  pivot_wider(names_from = direction, values_from = n, values_fill = 0) %>%
  kable(caption = "Count of FDR-significant words by category and direction.") %>%
  html_table_width(c(140, 100, 100))
Count of FDR-significant words by category and direction.
item_kind More likely: Monolingual
action_words 1
animals 4
body_parts 1
food_drink 1
games_routines 1
household 2
locations 1
outdoor 1
pronouns 1 
vi_word_group_diff %>%
  filter(p_fdr < 0.05, log_or > 0) %>%
  count(item_kind) %>%
  mutate(item_kind = fct_reorder(item_kind, n)) %>%
  ggplot(aes(x = n, y = item_kind)) +
  geom_col(fill = "#D55E00") +
  theme_classic() +
  xlab("Number of FDR-significant words") + ylab(NULL) +
  labs(title = "Where the 13 Monolingual-favoring words concentrate, by category")

Interpreting the asymmetry

At FDR < .05, 13 words are significantly over-represented among Vietnamese-only producers, and zero words are significantly over-represented among bilingual producers (the strongest bilingual-leaning candidate, “bắp rang bơ”/popcorn, only reaches p = .15, q = .31). This asymmetry is itself worth noting rather than a null result to skip past.

A look at the 13 Monolingual-favoring words is informative: many are ordinary, non-culturally-specific concepts a bilingual 2-3-year-old plausibly also knows – animals (gấu/bear, hươu/deer, gà trống/rooster, khỉ/monkey), household items (nệm/mattress), food (thịt/meat), a body part (cổ/neck). These don’t look like words that are specifically absent from a bilingual household’s Vietnamese input. A more likely explanation: the CDI only records production in the surveyed language. A bilingual child who knows the concept “bear” but has started saying “bear” in English rather than “gấu” in Vietnamese would show up as a Vietnamese non-producer for that item, even though they’re not missing the underlying concept – only the Vietnamese-specific expression of it. Under this reading, these 13 words could be less “words distinctive of monolingual family environments” and more “words where bilingual children have already shifted productive labeling to English” – a substantively different (and, if true, more interesting) claim than the original “representative vocabulary” framing. The check below tests this directly.

Testing the language-shift account

For each of the 13 significant Monolingual-favoring words, using data/EngVieConceptMatch.xlsx (a researcher-coded EN-VI concept-match table, more complete than the earlier hand-filled crosswalk): among the bilingual children who do not produce that Vietnamese word, what fraction do produce its English translation equivalent? A high fraction supports the language-shift account (they know it, just in English); a low fraction supports a genuine vocabulary gap (they don’t produce it in either language).

# EngVieConceptMatch.xlsx's `viet` column already matches vi_data$d_items$definition
# exactly (both are make.names()-sanitized the same way -- verified directly). Its
# `MB.eng.matched` column is the researcher-verified English match (as opposed to
# `eng.unmatched_1/2`, rougher candidate guesses) but has two quirks worth handling
# rather than ignoring: (1) "no match found" is coded as the *literal string* "NA",
# not a true missing value, so it must be converted explicitly; (2) a handful of
# Vietnamese words have multiple plausible English matches recorded as a single
# semicolon-separated string (e.g. "lấy" -> "take;get") -- these are kept as a list of
# candidates, and a child is counted as producing "the English equivalent" if they
# produce ANY one of them.
concept_match <- read_excel("data/EngVieConceptMatch.xlsx", sheet = "Sheet1") %>%
  transmute(
    vi_definition = viet,
    en_candidates = str_split(na_if(MB.eng.matched, "NA"), ";")
  )
## New names:
## • `` -> `...1`
## • `` -> `...6`
## • `` -> `...7`
sig_mono_words <- vi_word_group_diff %>%
  filter(p_fdr < 0.05, log_or > 0) %>%
  select(definition, item_definition, item_kind, log_or, p_fdr)

# Only candidate strings that actually exist as English items survive (2 of 673
# candidate strings across the full 687-item concept-match table are typos that don't
# match any real column, e.g. a missing "___" separator -- checked directly against
# en_data$d_items$definition); a word whose only candidate(s) fail this check is
# treated the same as having no usable match at all.
sig_mono_crosswalk <- sig_mono_words %>%
  left_join(concept_match, by = c("definition" = "vi_definition")) %>%
  mutate(en_candidates = map(en_candidates, ~ intersect(.x, en_data$d_items$definition)),
         n_en_candidates = lengths(en_candidates))

cat(nrow(sig_mono_words), "significant Monolingual-favoring words;",
    sum(sig_mono_crosswalk$n_en_candidates > 0), "have a usable translation match in EngVieConceptMatch.xlsx.\n")
## 13 significant Monolingual-favoring words; 11 have a usable translation match in EngVieConceptMatch.xlsx.
# Row-position alignment (vi_demo/vi_data$d_mat, en_data$d_demo/en_data$d_mat) follows
# the same invariant used throughout this project's Rmds (see load_lang_data() above).
vi_mat_by_row <- as_tibble(vi_data$d_mat) %>% mutate(row_id = vi_demo$row_id, group = vi_demo$group)
en_mat_by_row <- as_tibble(en_data$d_mat) %>% mutate(row_id = en_data$d_demo$row_id)

# A child counts as producing "the English equivalent" if they produce any one of
# (possibly several) candidate English items for that Vietnamese word.
check_language_shift <- function(vi_col, en_cols) {
  bi_non_producers <- vi_mat_by_row %>% filter(group == "Bilingual", .data[[vi_col]] == 0) %>% pull(row_id)
  en_sub <- en_mat_by_row %>% filter(row_id %in% bi_non_producers) %>% select(all_of(en_cols))
  produced_any <- if (ncol(en_sub) == 1) en_sub[[1]] == 1 else apply(en_sub == 1, 1, any)
  tibble(
    n_bilingual_non_producers   = length(bi_non_producers),
    n_produce_en_equivalent     = sum(produced_any),
    pct_produce_en_equivalent   = round(100 * mean(produced_any), 1)
  )
}

language_shift_results <- sig_mono_crosswalk %>%
  filter(n_en_candidates > 0) %>%
  mutate(
    en_definition = map_chr(en_candidates, ~ paste(str_remove(.x, "^.*?___"), collapse = " / ")),
    check = map2(definition, en_candidates, check_language_shift)
  ) %>%
  unnest(check) %>%
  mutate(account = if_else(pct_produce_en_equivalent >= 50, "Language shift", "Vocabulary gap"))

language_shift_results %>%
  select(item_definition, en_definition, item_kind, n_bilingual_non_producers,
         n_produce_en_equivalent, pct_produce_en_equivalent, account) %>%
  arrange(desc(pct_produce_en_equivalent)) %>%
  kable(digits = 1, caption = paste0(
    "Among bilingual non-producers of each Vietnamese word, % who produce its English equivalent ",
    "(>=50% = language-shift account favored).")) %>%
  html_table_width(c(120, 130, 80, 110, 110, 110, 100)) %>%
  row_spec(which(language_shift_results$account == "Language shift"), background = "#FFF3CD")
Among bilingual non-producers of each Vietnamese word, % who produce its English equivalent (>=50% = language-shift account favored).
item_definition en_definition item_kind n_bilingual_non_producers n_produce_en_equivalent pct_produce_en_equivalent account
tạm biệt bye games_routines 62 49 79.0 Language shift
hoa flower outdoor 58 21 36.2 Vocabulary gap
khỉ monkey animals 62 21 33.9 Vocabulary gap
gấu bear animals 60 20 33.3 Vocabulary gap
hươu deer animals 85 16 18.8 Vocabulary gap
dậy wake action_words 74 12 16.2 Vocabulary gap
gà trống rooster animals 83 12 14.5 Vocabulary gap
tiền money household 75 9 12.0 Vocabulary gap
của anh ấy his pronouns 94 9 9.6 Vocabulary gap
thịt meat food_drink 65 6 9.2 Vocabulary gap
cách xa away locations 96 8 8.3 Vocabulary gap 
language_shift_results %>%
  mutate(item_definition = fct_reorder(item_definition, pct_produce_en_equivalent)) %>%
  ggplot(aes(x = pct_produce_en_equivalent, y = item_definition, fill = account)) +
  geom_col() +
  geom_vline(xintercept = 50, linetype = "dashed", color = "grey30") +
  theme_classic() +
  scale_fill_manual(values = c(`Language shift` = "#E69F00", `Vocabulary gap` = "grey60")) +
  xlab("% of bilingual non-producers who DO produce the English equivalent") + ylab(NULL) +
  labs(fill = NULL, title = "Testing the language-shift account, word by word")

language_shift_results %>%
  count(account) %>%
  kable(caption = "Language-shift vs. vocabulary-gap account, tallied across the 13 significant words.") %>%
  html_table_width(c(150, 80))
Language-shift vs. vocabulary-gap account, tallied across the 13 significant words.
account n
Language shift 1
Vocabulary gap 10

Result: the data mostly contradicts the language-shift hypothesis. Only 1 of the 11 matched words (“tạm biệt” /bye, 79.0%) shows the pattern the language-shift account predicts. The other 10 – including every animal word in the list (khỉ/monkey 33.9%, gấu/bear 33.3%, hươu/deer 18.8%, gà trống/rooster 14.5%) and most of the rest (hoa/flower 36.2%, dậy/wake 16.2%, tiền/money 12.0%, của anh ấy/his 9.6%, thịt/meat 9.2%, cách xa/away 8.3%) – show that most bilingual non-producers of the Vietnamese word also don’t produce its English equivalent. So for most of these words, bilingual children in this sample aren’t quietly saying “bear” instead of “gấu” – they mostly don’t productively say either. That favors a genuine production gap specifically for these items (disproportionately animal vocabulary) among bilingual children, over the language-shift account, and reopens the original “representative vocabulary” framing rather than explaining it away: whatever is driving the Vietnamese-only group’s advantage on these specific words, it looks like it’s not just a reporting artifact of which language a bilingual child happens to answer in. (Two significant words, nệm/mattress and cổ/neck, had no usable English match in EngVieConceptMatch.xlsx and so aren’t part of this check.)

Acquisition trajectories for these 13 words

A third, developmental view of the same question: for these 13 words, how does production change with age in each of three group x language combinations – Bilingual children’s Vietnamese production, Bilingual children’s English production (of the crosswalk-matched equivalent, where one exists), and Monolingual children’s Vietnamese production? If this is a genuine production gap rather than a language-of-expression shift, the Bilingual-English line should stay low across the whole observed age range, not just lag behind and catch up later.

sig_words_tbl <- vi_word_group_diff %>%
  filter(p_fdr < 0.05, log_or > 0) %>%
  select(definition, item_definition)

# Vietnamese production of each of the 13 words, by child, with that child's group and age.
vi_traj_long <- as_tibble(vi_data$d_mat[, sig_words_tbl$definition, drop = FALSE]) %>%
  mutate(row_id = vi_demo$row_id, age = vi_demo$age, group = vi_demo$group) %>%
  pivot_longer(cols = all_of(sig_words_tbl$definition), names_to = "definition", values_to = "produces") %>%
  left_join(sig_words_tbl, by = "definition") %>%
  mutate(language = "Vietnamese")

# English production of each word's crosswalk-matched equivalent (where one exists),
# for Bilingual children only -- reuses sig_mono_crosswalk's candidate-column lookup
# from the language-shift check above (a child counts as producing the "equivalent"
# if they produce ANY of that word's candidate English items).
en_traj_long <- sig_mono_crosswalk %>%
  filter(n_en_candidates > 0) %>%
  select(definition, item_definition, en_candidates) %>%
  pmap_dfr(function(definition, item_definition, en_candidates) {
    sub <- en_data$d_mat[, en_candidates, drop = FALSE]
    produces <- if (ncol(sub) == 1) as.integer(sub[, 1]) else as.integer(apply(sub == 1, 1, any))
    tibble(definition = definition, item_definition = item_definition,
           row_id = en_data$d_demo$row_id, age = en_data$d_demo$age, produces = produces)
  }) %>%
  left_join(vi_demo %>% select(row_id, group), by = "row_id") %>%
  filter(group == "Bilingual") %>%
  mutate(language = "English")

traj_combined <- bind_rows(vi_traj_long, en_traj_long) %>%
  mutate(line = paste0(group, " – ", language))

cat("Words with no English trajectory line (no crosswalk match):",
    paste(setdiff(sig_words_tbl$item_definition, unique(en_traj_long$item_definition)), collapse = ", "), "\n")
## Words with no English trajectory line (no crosswalk match): cổ, nệm

Aggregate trajectory (averaged across all 13 words)

age_bin_width <- 3  # wider than 05_item_difficulty_crosslang.Rmd's category trajectories (2 months),
                     # since the Monolingual group here is only 41 children spread across 13 words

traj_combined %>%
  mutate(age_bin = floor(age / age_bin_width) * age_bin_width) %>%
  group_by(group, language, line, age_bin) %>%
  summarise(prop_producing = mean(produces), n_obs = n(), .groups = "drop") %>%
  ggplot(aes(x = age_bin, y = prop_producing, color = group, linetype = language)) +
  geom_line(linewidth = 0.8) +
  geom_point(size = 1.8) +
  theme_classic() +
  scale_color_manual(values = c(Bilingual = "#1F78B4", Monolingual = "#E31A1C")) +
  xlab("Age (months, binned)") + ylab("Mean proportion producing (across 13 words)") +
  labs(color = NULL, linetype = NULL,
       title = "Acquisition trajectory of the 13 Monolingual-favoring words")

Per-word trajectories

traj_combined %>%
  mutate(age_bin = floor(age / age_bin_width) * age_bin_width) %>%
  group_by(item_definition, group, language, line, age_bin) %>%
  summarise(prop_producing = mean(produces), n_obs = n(), .groups = "drop") %>%
  ggplot(aes(x = age_bin, y = prop_producing, color = group, linetype = language)) +
  geom_line(linewidth = 0.6) +
  geom_point(size = 1) +
  facet_wrap(~item_definition, ncol = 4) +
  theme_classic() +
  theme(legend.position = "bottom") +
  scale_color_manual(values = c(Bilingual = "#1F78B4", Monolingual = "#E31A1C")) +
  xlab("Age (months, binned)") + ylab("Proportion producing") +
  labs(color = NULL, linetype = NULL)

Reading note: per-word bins can be thin (as few as 1-3 Monolingual children in an age bin), so individual-word curves are noisy by construction – the aggregate plot above is the more reliable summary; per-word panels are for spotting words that depart from the aggregate pattern (e.g. a word already near-ceiling in both groups by the youngest ages sampled, vs. one where the gap only opens up later).

Does Vietnamese exposure explain the gap?

The leading account for why these particular 13 words differ (worked out above: they’re moderate-difficulty, already-common-for-Monolinguals words, with animal vocabulary over-represented ~4.6x relative to its base rate) is input-dose dilution: Bilingual children split their waking hours across two languages, so for any given word they receive roughly half the cumulative Vietnamese exposure a Vietnamese-dominant peer gets – which mostly only matters for this moderate-difficulty band (the earliest/highest-frequency words reach ceiling under either amount of exposure; the hardest words are below floor for both groups regardless). That account makes a specific, testable prediction: among Bilingual children, more reported Vietnamese exposure should predict producing more of these 13 words – i.e. the gap should be exposure-graded, not a flat Bilingual-vs-Monolingual step function.

# Mirrors get_pct_active_lang() in 03_bifactor_noun_verb.Rmd: lang_i_pct_active is
# the % of a child's *active* language use taken up by whichever language sits in
# slot i (language_1..4); summing across slots where that slot is "Vietnamese"
# gives each child's overall % active Vietnamese use.
get_pct_active_lang <- function(d_demo, lang_name) {
  lang_mat <- as.matrix(d_demo[, c("language_1", "language_2", "language_3", "language_4")])
  pct_mat  <- as.matrix(d_demo[, c("lang_1_pct_active", "lang_2_pct_active", "lang_3_pct_active", "lang_4_pct_active")])
  is_match <- lang_mat == lang_name
  is_match[is.na(is_match)] <- FALSE
  pct_mat[!is_match] <- 0
  rowSums(pct_mat, na.rm = TRUE)
}

vi_exposure <- read_csv("data/VietnameseWS_Bui_data.csv", show_col_types = FALSE) %>%
  select(response_id, language_1, language_2, language_3, language_4,
         lang_1_pct_active, lang_2_pct_active, lang_3_pct_active, lang_4_pct_active) %>%
  mutate(pct_active_vietnamese = get_pct_active_lang(., "Vietnamese")) %>%
  select(response_id, pct_active_vietnamese)

vi_demo <- vi_demo %>% left_join(vi_exposure, by = "response_id")

# Sanity check: Monolingual children should report ~100% active Vietnamese use,
# validating pct_active_vietnamese = 100 as the "full exposure" reference point
# the plot below compares Bilingual children against.
vi_demo %>% group_by(group) %>%
  summarise(mean_pct_active_vi = mean(pct_active_vietnamese, na.rm = TRUE),
            median_pct_active_vi = median(pct_active_vietnamese, na.rm = TRUE), .groups = "drop") %>%
  kable(digits = 1, caption = "% active Vietnamese use, by group.") %>%
  html_table_width(c(90, 110, 110))
% active Vietnamese use, by group.
group mean_pct_active_vi median_pct_active_vi
Bilingual 52.7 50
Monolingual 69.3 95 
n_sig_words_produced <- as_tibble(vi_data$d_mat[, sig_words_tbl$definition, drop = FALSE]) %>%
  mutate(row_id = vi_demo$row_id) %>%
  pivot_longer(-row_id, names_to = "definition", values_to = "produces") %>%
  group_by(row_id) %>%
  summarise(n_produced = sum(produces), prop_produced = mean(produces), .groups = "drop")

vi_demo_exposure <- vi_demo %>% left_join(n_sig_words_produced, by = "row_id")

monolingual_baseline <- vi_demo_exposure %>% filter(group == "Monolingual") %>% pull(prop_produced) %>% mean()
cat("Monolingual mean proportion of the 13 words produced (reference line below):",
    round(monolingual_baseline, 3), "\n")
## Monolingual mean proportion of the 13 words produced (reference line below): 0.57
bilingual_exposure <- vi_demo_exposure %>% filter(group == "Bilingual")

summary(lm(prop_produced ~ pct_active_vietnamese + age_c, data = bilingual_exposure))
## 
## Call:
## lm(formula = prop_produced ~ pct_active_vietnamese + age_c, data = bilingual_exposure)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.50196 -0.15545 -0.05066  0.16408  0.63192 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           0.0436849  0.0543616   0.804    0.424    
## pct_active_vietnamese 0.0042526  0.0009245   4.600 1.26e-05 ***
## age_c                 0.0209583  0.0040919   5.122 1.51e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.244 on 98 degrees of freedom
## Multiple R-squared:  0.2923, Adjusted R-squared:  0.2779 
## F-statistic: 20.24 on 2 and 98 DF,  p-value: 4.38e-08
ggplot(bilingual_exposure, aes(x = pct_active_vietnamese, y = prop_produced)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = TRUE, color = "#1F78B4") +
  geom_hline(yintercept = monolingual_baseline, linetype = "dashed", color = "#E31A1C") +
  annotate("text", x = 5, y = monolingual_baseline + 0.03, label = "Monolingual mean",
           color = "#E31A1C", hjust = 0, size = 3.5) +
  theme_classic() +
  xlab("% active Vietnamese use (Bilingual children)") +
  ylab("Proportion of the 13 words produced") +
  labs(title = "Does higher Vietnamese exposure close the gap?")
## `geom_smooth()` using formula = 'y ~ x'

If the input-dilution account is right, the regression line should rise with pct_active_vietnamese and approach the Monolingual reference line (dashed) at the high-exposure end; a flat line would instead suggest the gap is tied to group membership itself (or something else correlated with it) rather than graded by exposure amount.

Is this specific to these 13 words, or just a generic vocabulary-size effect?

Exposure predicting production of these 13 words isn’t informative on its own if exposure predicts everything equally well – overall Vietnamese vocabulary size should obviously track Vietnamese exposure too. The real question is whether these 13 words are more exposure-sensitive than a typical Vietnamese word, as the input-dilution + “moderate-difficulty band” account would predict, or whether they just ride along with a generic effect that’s no stronger here than anywhere else in the vocabulary. Three outcomes, same model (~ pct_active_vietnamese + age_c), same Bilingual children, compared on a common scale (standardized exposure effect, i.e. the correlation between pct_active_vietnamese and each outcome controlling for age):

non_sig_words <- setdiff(colnames(vi_data$d_mat), sig_words_tbl$definition)

bilingual_exposure <- bilingual_exposure %>%
  mutate(
    prop_produced_13       = prop_produced,  # the 13 significant words, already computed above
    prop_produced_all      = rowSums(vi_data$d_mat[vi_demo$group == "Bilingual", ]) / ncol(vi_data$d_mat),
    prop_produced_nonsig13 = rowSums(vi_data$d_mat[vi_demo$group == "Bilingual", non_sig_words]) / length(non_sig_words)
  )

fit_and_summarize <- function(outcome_var, label) {
  f <- as.formula(paste0(outcome_var, " ~ pct_active_vietnamese + age_c"))
  fit <- lm(f, data = bilingual_exposure)
  cf <- summary(fit)$coefficients
  # partial correlation of pct_active_vietnamese with the outcome, controlling for age_c
  partial_r <- sign(cf["pct_active_vietnamese", "Estimate"]) *
    sqrt(cf["pct_active_vietnamese", "t value"]^2 / (cf["pct_active_vietnamese", "t value"]^2 + fit$df.residual))
  tibble(outcome = label, n_items = NA_integer_,
         beta_exposure = cf["pct_active_vietnamese", "Estimate"],
         p_exposure = cf["pct_active_vietnamese", "Pr(>|t|)"],
         partial_r_exposure = partial_r,
         r_squared = summary(fit)$r.squared)
}

comparison <- bind_rows(
  fit_and_summarize("prop_produced_13", "The 13 significant words") %>% mutate(n_items = 13),
  fit_and_summarize("prop_produced_nonsig13", "All other Vietnamese words") %>% mutate(n_items = length(non_sig_words)),
  fit_and_summarize("prop_produced_all", "All 687 Vietnamese words (incl. the 13)") %>% mutate(n_items = ncol(vi_data$d_mat))
)

comparison %>%
  select(outcome, n_items, beta_exposure, p_exposure, partial_r_exposure, r_squared) %>%
  kable(digits = c(0, 0, 4, 4, 2, 2), caption = "Vietnamese-exposure effect on production, by word set (Bilingual children, controlling for age).") %>%
  html_table_width(c(220, 70, 100, 90, 110, 80))
Vietnamese-exposure effect on production, by word set (Bilingual children, controlling for age).
outcome n_items beta_exposure p_exposure partial_r_exposure r_squared
The 13 significant words 13 0.0043 0e+00 0.42 0.29
All other Vietnamese words 674 0.0029 4e-04 0.35 0.27
All 687 Vietnamese words (incl. the 13) 687 0.0030 4e-04 0.35 0.27 
bilingual_exposure %>%
  select(pct_active_vietnamese, prop_produced_13, prop_produced_nonsig13, prop_produced_all) %>%
  pivot_longer(-pct_active_vietnamese, names_to = "word_set", values_to = "prop_produced") %>%
  mutate(word_set = recode(word_set,
                            prop_produced_13 = "The 13 significant words",
                            prop_produced_nonsig13 = "All other Vietnamese words",
                            prop_produced_all = "All 687 words")) %>%
  ggplot(aes(x = pct_active_vietnamese, y = prop_produced, color = word_set)) +
  geom_smooth(method = "lm", se = TRUE) +
  theme_classic() +
  scale_color_manual(values = c("The 13 significant words" = "#E31A1C",
                                 "All other Vietnamese words" = "#1F78B4",
                                 "All 687 words" = "grey50")) +
  xlab("% active Vietnamese use (Bilingual children)") +
  ylab("Proportion of word set produced") +
  labs(color = NULL, title = "Exposure effect: the 13 words vs. the rest of the vocabulary")
## `geom_smooth()` using formula = 'y ~ x'

Lines are plotted on each outcome’s own proportion scale, so compare steepness relative to each line’s own range (or just the partial_r_exposure/r_squared columns above) rather than raw vertical position – the 13-word outcome and the 674-word outcome aren’t on directly comparable absolute scales (different difficulty composition), which is exactly why the table reports a standardized partial correlation instead of relying on the plot alone.

Permutation test: are these 13 words unusually exposure-sensitive?

The comparison above (13 words vs. the other 674) isn’t quite a fair test on its own: a 13-item average is noisier per child than a 674-item average purely from having fewer items to average over, and that extra noise should, if anything, attenuate a correlation with exposure (classic measurement-error dilution) – so it’s not obvious whether partial r = .42 is “high for a 13-word subset” or just typical of what any 13-word subset would show once you account for that smaller-sample noise. The direct fix: build a null distribution from many random same-sized (13-word) subsets of the Vietnamese vocabulary, and see where the real 13 words actually fall in it.

set.seed(2025)
n_perm <- 5000

bilingual_idx <- which(vi_demo$group == "Bilingual")
bilingual_mat <- vi_data$d_mat[bilingual_idx, ]
all_word_cols <- colnames(vi_data$d_mat)

# Same partial-correlation statistic as fit_and_summarize() above, computed for an
# arbitrary set of word columns against the fixed bilingual_exposure covariates.
compute_partial_r <- function(word_cols) {
  prop <- rowSums(bilingual_mat[, word_cols, drop = FALSE]) / length(word_cols)
  fit <- lm(prop ~ pct_active_vietnamese + age_c, data = bilingual_exposure)
  t <- summary(fit)$coefficients["pct_active_vietnamese", "t value"]
  sign(t) * sqrt(t^2 / (t^2 + fit$df.residual))
}

observed_r <- compute_partial_r(sig_words_tbl$definition)

null_r <- map_dbl(seq_len(n_perm), function(i) {
  compute_partial_r(sample(all_word_cols, length(sig_words_tbl$definition)))
})

p_perm <- mean(null_r >= observed_r)
cat("Observed partial r (the 13 significant words):", round(observed_r, 3), "\n")
## Observed partial r (the 13 significant words): 0.421
cat("Null distribution (", n_perm, "random 13-word subsets): mean =", round(mean(null_r), 3),
    ", 95th percentile =", round(quantile(null_r, .95), 3), "\n")
## Null distribution ( 5000 random 13-word subsets): mean = 0.327 , 95th percentile = 0.383
cat("One-sided permutation p-value:", round(p_perm, 4), "\n")
## One-sided permutation p-value: 0.0028
tibble(null_r = null_r) %>%
  ggplot(aes(x = null_r)) +
  geom_histogram(bins = 50, fill = "grey75", color = "white") +
  geom_vline(xintercept = observed_r, color = "#E31A1C", linewidth = 1) +
  annotate("text", x = observed_r, y = Inf, label = "The 13 words", color = "#E31A1C",
           hjust = -0.1, vjust = 1.5, size = 3.8) +
  theme_classic() +
  xlab("Partial r (exposure effect), random 13-word subset") +
  ylab("Count (of 5,000 random draws)") +
  labs(title = "Null distribution: exposure-sensitivity of random 13-word subsets")

If p_perm is small (conventionally < .05), that’s evidence these particular 13 words are more exposure-sensitive than a typical same-sized slice of the Vietnamese vocabulary – not just an artifact of averaging over fewer items. A large p_perm would mean the apparent specialness in the earlier comparison was mostly a small-subset-noise artifact, and the honest claim is the more modest one from the “all 674 vs. all 687” comparison above: exposure matters for Vietnamese vocabulary generally, without these 13 words being a distinct category.

Controlling for item difficulty: is this just a difficulty-band effect?

The permutation test above draws random 13-word subsets from the entire vocabulary, uniformly. But IRT item information is mechanically highest near an item’s own difficulty: a word that sits right at the ability level where children are actively transitioning from “doesn’t produce” to “produces” is, by construction, maximally sensitive to any small individual difference among children that nudges them one way or the other – including exposure – regardless of the word’s content. Since these 13 words skew toward easier-than-median difficulty (median b = 0.82 vs. 1.29 for the full vocabulary), part of their apparent exposure-sensitivity could be nothing more than “words in this difficulty range are sensitive to everything,” not anything specific to these words. The fix: instead of a free random draw, match each of the 13 real words to a pool of same-difficulty other words, and redraw from those pools only.

load("models/vi/mod_2pl.Rds")  # multi-object save() file -- see load() note in 05_item_difficulty_crosslang.Rmd
diff_vi_full <- as_tibble(coef(mod_2pl, simplify = TRUE, IRTpars = TRUE)$items) %>%
  mutate(definition = rownames(coef(mod_2pl, simplify = TRUE, IRTpars = TRUE)$items)) %>%
  select(definition, b)
rm(mod_2pl)

sig_difficulty <- diff_vi_full %>% filter(definition %in% sig_words_tbl$definition)

sig_difficulty %>%
  arrange(b) %>%
  left_join(sig_words_tbl, by = "definition") %>%
  select(item_definition, b) %>%
  kable(digits = 2, caption = "IRT difficulty (b) of the 13 significant words -- a wide spread, not a narrow band.") %>%
  html_table_width(c(120, 80))
IRT difficulty (b) of the 13 significant words – a wide spread, not a narrow band.
item_definition b
hoa 0.08
gấu 0.22
tạm biệt 0.23
khỉ 0.33
cổ 0.36
thịt 0.43
dậy 0.82
tiền 0.87
gà trống 1.26
nệm 1.28
hươu 1.34
của anh ấy 2.10
cách xa 2.41 
# For each of the 13 words, its 30 nearest neighbors by |difference in b|, excluding
# the 13 significant words themselves -- a k-nearest-neighbor pool rather than a fixed
# +-tolerance window, so it stays well-populated even for the hardest/easiest of the 13
# (checked: 28-166 words fall within +-0.15 to +-0.30 logits across this range, so k=30
# neighbors corresponds to a reasonably tight difficulty match throughout).
k_neighbors <- 30
match_pools <- map(sig_difficulty$b, function(b_target) {
  diff_vi_full %>%
    filter(!(definition %in% sig_words_tbl$definition)) %>%
    mutate(dist = abs(b - b_target)) %>%
    slice_min(dist, n = k_neighbors) %>%
    pull(definition)
})

cat("Pool size per word (all should be", k_neighbors, "):", paste(lengths(match_pools), collapse = ", "), "\n")
## Pool size per word (all should be 30 ): 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30
set.seed(2025)

draw_matched_set <- function(pools) {
  drawn <- character(0)
  for (pool in pools) {
    available <- setdiff(pool, drawn)
    if (length(available) == 0) available <- pool  # fallback; not expected to trigger with k=30
    drawn <- c(drawn, sample(available, 1))
  }
  drawn
}

null_r_matched <- map_dbl(seq_len(n_perm), function(i) {
  compute_partial_r(draw_matched_set(match_pools))
})

p_perm_matched <- mean(null_r_matched >= observed_r)
cat("Observed partial r:", round(observed_r, 3), "\n")
## Observed partial r: 0.421
cat("Difficulty-matched null (", n_perm, "draws): mean =", round(mean(null_r_matched), 3),
    ", 95th percentile =", round(quantile(null_r_matched, .95), 3), "\n")
## Difficulty-matched null ( 5000 draws): mean = 0.331 , 95th percentile = 0.378
cat("One-sided permutation p-value (difficulty-matched):", round(p_perm_matched, 4), "\n")
## One-sided permutation p-value (difficulty-matched): 6e-04
bind_rows(
  tibble(r = null_r, null_type = "Unmatched (any 13 words)"),
  tibble(r = null_r_matched, null_type = "Difficulty-matched")
) %>%
  ggplot(aes(x = r, fill = null_type)) +
  geom_histogram(bins = 50, alpha = 0.6, position = "identity") +
  geom_vline(xintercept = observed_r, color = "#E31A1C", linewidth = 1) +
  annotate("text", x = observed_r, y = Inf, label = "The 13 words", color = "#E31A1C",
           hjust = -0.1, vjust = 1.5, size = 3.8) +
  theme_classic() +
  scale_fill_manual(values = c("Unmatched (any 13 words)" = "grey60", "Difficulty-matched" = "#1F78B4")) +
  xlab("Partial r (exposure effect), random 13-word subset") +
  ylab("Count (of 5,000 random draws)") +
  labs(fill = NULL, title = "Does difficulty-matching absorb the effect?")

This directly answers both of the opening questions. If the difficulty-matched null shifts right (closer to the unmatched null’s right tail, or to observed_r itself) and p_perm_matched is no longer small, then yes, this is just a difficulty effect – any 687-item-vocabulary word near this difficulty range would show the same exposure- sensitivity, and there’s nothing else special about these 13. If p_perm_matched stays small even after matching, that rules out difficulty as the (sole) explanation: these words are more exposure-sensitive than other words at the same difficulty, which points back to something about their content (the animal-word enrichment noted earlier being the most concrete candidate) rather than a generic psychometric artifact of where they sit on the difficulty scale.

Which words are most exposure-sensitive, in each language?

Everything above started from the 13 words that happened to come out of the Monolingual-vs-Bilingual comparison. This section instead runs an unbiased sweep: for every word in both languages, among Bilingual children only, does production depend on how much of that specific language they’re reportedly exposed to? This both (a) checks whether the 13 words from before actually show up among the most exposure-sensitive words in an analysis that didn’t presuppose them, and (b) surfaces an English-side exposure-sensitive word list, which nothing above has looked at yet (everything so far was Vietnamese-only).

# pct_active_english, built the same way as pct_active_vietnamese earlier
# (get_pct_active_lang(), defined in the "Does Vietnamese exposure explain the gap?"
# section above) -- language_1..4 / lang_i_pct_active are demographic fields
# duplicated identically across the English and Vietnamese raw data files (verified
# directly: identical for every shared row_id), so either file would give the same answer.
pct_active_english_tbl <- read_csv("data/EnglishAmericanWS_Bui_data_redact.csv", show_col_types = FALSE) %>%
  select(response_id, language_1, language_2, language_3, language_4,
         lang_1_pct_active, lang_2_pct_active, lang_3_pct_active, lang_4_pct_active) %>%
  mutate(pct_active_english = get_pct_active_lang(., "English")) %>%
  select(response_id, pct_active_english)

vi_demo <- vi_demo %>% left_join(pct_active_english_tbl, by = "response_id")
bilingual_full <- vi_demo %>% filter(group == "Bilingual")

# en_data$d_mat and vi_data$d_mat share row order with vi_demo (same invariant used
# throughout this project: both raw CSVs were confirmed row-for-row aligned by
# response_id earlier in this analysis), so the same logical index subsets all three.
is_bilingual <- vi_demo$group == "Bilingual"
vi_bilingual_mat <- vi_data$d_mat[is_bilingual, ]
en_bilingual_mat <- en_data$d_mat[is_bilingual, ]
# Per-item logistic regression of production on that language's exposure measure,
# controlling for age -- the whole-vocabulary analogue of the single-item checks
# above. Items with zero variance, non-convergent fits, or quasi-separation
# (SE > 5, same guard used for compute_item_aoa()/DIF elsewhere in this project) are
# returned as NA rather than trusted.
compute_exposure_effect <- function(d_mat, d_items, exposure, age_c) {
  n_items <- ncol(d_mat)
  defs <- colnames(d_mat)
  map_dfr(seq_len(n_items), function(j) {
    y <- d_mat[, j]
    if (length(unique(y)) < 2) return(tibble(definition = defs[j], beta_exposure = NA_real_, p_exposure = NA_real_))
    fit <- tryCatch(glm(y ~ exposure + age_c, family = binomial), error = function(e) NULL)
    if (is.null(fit) || !fit$converged) return(tibble(definition = defs[j], beta_exposure = NA_real_, p_exposure = NA_real_))
    cf <- summary(fit)$coefficients
    if (!("exposure" %in% rownames(cf))) return(tibble(definition = defs[j], beta_exposure = NA_real_, p_exposure = NA_real_))
    est <- cf["exposure", "Estimate"]; se <- cf["exposure", "Std. Error"]
    if (!is.finite(se) || se > 5) return(tibble(definition = defs[j], beta_exposure = NA_real_, p_exposure = NA_real_))
    tibble(definition = defs[j], beta_exposure = est, p_exposure = cf["exposure", "Pr(>|z|)"])
  }) %>%
    mutate(p_fdr = p.adjust(p_exposure, method = "BH")) %>%
    left_join(d_items, by = "definition")
}

vi_exposure_effects <- compute_exposure_effect(vi_bilingual_mat, vi_data$d_items,
                                                bilingual_full$pct_active_vietnamese, bilingual_full$age_c)
en_exposure_effects <- compute_exposure_effect(en_bilingual_mat, en_data$d_items,
                                                bilingual_full$pct_active_english, bilingual_full$age_c)

cat("Vietnamese: ", sum(vi_exposure_effects$p_fdr < 0.05, na.rm = TRUE), "/", nrow(vi_exposure_effects),
    " words FDR-significant for exposure effect\n", sep = "")
## Vietnamese: 163/687 words FDR-significant for exposure effect
cat("English: ", sum(en_exposure_effects$p_fdr < 0.05, na.rm = TRUE), "/", nrow(en_exposure_effects),
    " words FDR-significant for exposure effect\n", sep = "")
## English: 345/681 words FDR-significant for exposure effect

Top 20 most exposure-sensitive words, Vietnamese

vi_exposure_effects %>%
  filter(!is.na(p_fdr)) %>%
  arrange(p_fdr) %>%
  slice_head(n = 20) %>%
  select(item_definition, item_kind, beta_exposure, p_exposure, p_fdr) %>%
  kable(digits = c(0, 0, 4, 4, 4), caption = "Vietnamese words most sensitive to % active Vietnamese exposure (Bilingual children).") %>%
  html_table_width(c(140, 90, 90, 90, 90))
Vietnamese words most sensitive to % active Vietnamese exposure (Bilingual children).
item_definition item_kind beta_exposure p_exposure p_fdr
kính household 0.0451 0.0000 0.0332
hổ animals 0.0491 0.0001 0.0344
chơi action_words 0.0274 0.0020 0.0350
cắn action_words 0.0292 0.0018 0.0350
giấu action_words 0.0361 0.0021 0.0350
nhặt action_words 0.0361 0.0014 0.0350
chuột animals 0.0356 0.0004 0.0350
cú animals 0.0369 0.0021 0.0350
cừu animals 0.0460 0.0002 0.0350
gấu animals 0.0355 0.0004 0.0350
gấu trúc animals 0.0384 0.0011 0.0350
hươu animals 0.0454 0.0012 0.0350
hươu cao cổ animals 0.0383 0.0005 0.0350
kiến animals 0.0350 0.0010 0.0350
rùa animals 0.0319 0.0021 0.0350
sóc animals 0.0374 0.0010 0.0350
sư tử animals 0.0346 0.0007 0.0350
rốn body_parts 0.0392 0.0004 0.0350
cứng descriptive_words 0.0421 0.0015 0.0350
nặng descriptive_words 0.0391 0.0014 0.0350

Top 20 most exposure-sensitive words, English

en_exposure_effects %>%
  filter(!is.na(p_fdr)) %>%
  arrange(p_fdr) %>%
  slice_head(n = 20) %>%
  select(item_definition, item_kind, beta_exposure, p_exposure, p_fdr) %>%
  kable(digits = c(0, 0, 4, 4, 4), caption = "English words most sensitive to % active English exposure (Bilingual children).") %>%
  html_table_width(c(140, 90, 90, 90, 90))
English words most sensitive to % active English exposure (Bilingual children).
item_definition item_kind beta_exposure p_exposure p_fdr
swing action_words 0.0446 4e-04 0.0166
owl animals 0.0335 6e-04 0.0166
cheek body_parts 0.0382 2e-04 0.0166
chin body_parts 0.0327 7e-04 0.0166
leg body_parts 0.0361 6e-04 0.0166
tongue body_parts 0.0467 3e-04 0.0166
tooth body_parts 0.0332 7e-04 0.0166
pants clothing 0.0342 6e-04 0.0166
shirt clothing 0.0357 4e-04 0.0166
shoe clothing 0.0399 1e-04 0.0166
sock clothing 0.0316 4e-04 0.0166
bad descriptive_words 0.0472 6e-04 0.0166
happy descriptive_words 0.0355 3e-04 0.0166
milk food_drink 0.0366 2e-04 0.0166
gonna get you! games_routines 0.0472 4e-04 0.0166
please games_routines 0.0346 3e-04 0.0166
basket household 0.0487 6e-04 0.0166
bottle household 0.0342 6e-04 0.0166
bowl household 0.0396 3e-04 0.0166
pillow household 0.0372 5e-04 0.0166

Do the original 13 words show up here?

vi_exposure_effects %>%
  filter(definition %in% sig_words_tbl$definition) %>%
  mutate(rank_out_of = nrow(vi_exposure_effects)) %>%
  left_join(vi_exposure_effects %>% arrange(p_fdr) %>% mutate(rank = row_number()) %>% select(definition, rank),
            by = "definition") %>%
  arrange(rank) %>%
  select(item_definition, item_kind, beta_exposure, p_fdr, rank, rank_out_of) %>%
  kable(digits = c(0, 0, 4, 4, 0, 0),
        caption = "Where the original 13 words rank in this unbiased, whole-vocabulary exposure sweep.") %>%
  html_table_width(c(140, 90, 90, 90, 70, 90))
Where the original 13 words rank in this unbiased, whole-vocabulary exposure sweep.
item_definition item_kind beta_exposure p_fdr rank rank_out_of
gấu animals 0.0355 0.0350 10 687
hươu animals 0.0454 0.0350 12 687
tạm biệt games_routines 0.0343 0.0350 32 687
nệm household 0.0440 0.0352 49 687
gà trống animals 0.0343 0.0388 67 687
khỉ animals 0.0237 0.0495 143 687
hoa outdoor 0.0215 0.0495 152 687
dậy action_words 0.0247 0.0496 156 687
cổ body_parts 0.0213 0.0543 184 687
thịt food_drink 0.0183 0.0842 288 687
tiền household 0.0193 0.0912 314 687
của anh ấy pronouns 0.0372 0.0923 320 687
cách xa locations 0.0194 0.3341 606 687

Volcano plots, both languages

bind_rows(
  vi_exposure_effects %>% mutate(language = "Vietnamese"),
  en_exposure_effects %>% mutate(language = "English")
) %>%
  filter(!is.na(p_exposure)) %>%
  ggplot(aes(x = beta_exposure, y = -log10(p_exposure), color = p_fdr < 0.05)) +
  geom_point(alpha = 0.5, size = 1) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
  facet_wrap(~language, scales = "free_x") +
  theme_classic() +
  scale_color_manual(values = c(`TRUE` = "#D55E00", `FALSE` = "grey70")) +
  xlab("Exposure effect (logistic regression coefficient)") +
  ylab(expression(-log[10]("p-value"))) +
  labs(color = "FDR < .05", title = "Exposure-sensitivity across the whole vocabulary, by language")

Does item difficulty moderate the exposure effect itself?

A sharper, whole-vocabulary version of the difficulty question from earlier: rather than asking whether the 13 words’ competition effect correlates with difficulty (RQ4, tested on the EN-only/VI-only subset), this asks whether each item’s raw exposure-sensitivity (beta_exposure, from the per-item regressions just computed, across all 687 Vietnamese and 681 English words) depends on its IRT difficulty – and specifically whether that dependence is an inverted U (peak sensitivity at moderate difficulty, low at both ceiling and floor), the shape predicted by IRT item information being maximized near an item’s own difficulty, rather than a simple linear trend.

# English difficulty, loaded the same way as diff_vi_full above (05_item_difficulty_crosslang.Rmd's convention).
load("models/en/mod_2pl.Rds")
diff_en_full <- as_tibble(coef(mod_2pl, simplify = TRUE, IRTpars = TRUE)$items) %>%
  mutate(definition = rownames(coef(mod_2pl, simplify = TRUE, IRTpars = TRUE)$items)) %>%
  select(definition, b)
rm(mod_2pl)

vi_moderation <- vi_exposure_effects %>%
  filter(!is.na(beta_exposure)) %>%
  left_join(diff_vi_full %>% select(definition, b), by = "definition") %>%
  mutate(language = "Vietnamese")
en_moderation <- en_exposure_effects %>%
  filter(!is.na(beta_exposure)) %>%
  left_join(diff_en_full %>% select(definition, b), by = "definition") %>%
  mutate(language = "English")
fit_moderation <- function(d, label) {
  lin <- lm(beta_exposure ~ b, data = d)
  quad <- lm(beta_exposure ~ b + I(b^2), data = d)
  comparison <- anova(lin, quad)
  vertex <- -coef(quad)["b"] / (2 * coef(quad)["I(b^2)"])
  cat("===", label, "===\n")
  cat("Linear b coefficient:", round(coef(lin)["b"], 4), ", p =", round(summary(lin)$coefficients["b","Pr(>|t|)"], 4), "\n")
  cat("Quadratic b^2 coefficient:", round(coef(quad)["I(b^2)"], 5),
      ", p =", round(summary(quad)$coefficients["I(b^2)","Pr(>|t|)"], 4), "\n")
  cat("Does adding the quadratic term improve fit? p =", round(comparison$`Pr(>F)`[2], 4), "\n")
  cat("Implied peak-sensitivity difficulty (vertex):", round(vertex, 2),
      "(negative b^2 coefficient means this is a maximum, not a minimum)\n\n")
}

fit_moderation(vi_moderation, "Vietnamese")
## === Vietnamese ===
## Linear b coefficient: 0.0034 , p = 0 
## Quadratic b^2 coefficient: -0.00095 , p = 0.002 
## Does adding the quadratic term improve fit? p = 0.002 
## Implied peak-sensitivity difficulty (vertex): 2.92 (negative b^2 coefficient means this is a maximum, not a minimum)
fit_moderation(en_moderation, "English")
## === English ===
## Linear b coefficient: 0.0042 , p = 0 
## Quadratic b^2 coefficient: -0.00179 , p = 0.0024 
## Does adding the quadratic term improve fit? p = 0.0024 
## Implied peak-sensitivity difficulty (vertex): 3.06 (negative b^2 coefficient means this is a maximum, not a minimum)
bind_rows(vi_moderation, en_moderation) %>%
  ggplot(aes(x = b, y = beta_exposure)) +
  geom_point(alpha = 0.3, size = 1) +
  geom_smooth(method = "lm", formula = y ~ x + I(x^2), se = TRUE, color = "#E31A1C") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
  facet_wrap(~language, scales = "free") +
  theme_classic() +
  xlab("IRT difficulty (b)") +
  ylab("Exposure-sensitivity (beta_exposure)") +
  labs(title = "Does exposure-sensitivity peak at moderate difficulty?")

If the quadratic term is significant and negative (and the model-comparison p-value is small), that’s direct evidence difficulty moderates the exposure effect in the inverted-U shape IRT information predicts. If it isn’t, exposure-sensitivity doesn’t simply track “how mechanically sensitive this difficulty level is” – reinforcing, from a different angle than the earlier matched-permutation test, that content (not psychometric position alone) is doing real work.

Language competition: does exposure balance determine which language a concept gets produced in?

Everything above treated each language’s vocabulary growth as if it were independent – more Vietnamese exposure predicts more Vietnamese words, more English exposure predicts more English words. A different, stronger claim is competition: for a single concept a child knows (e.g. “dog”), does which language they happen to produce it in depend on their relative exposure balance, as if the two languages’ words for the same concept compete for a single “production slot”? The clean test: restrict to children who produce a concept in exactly one of the two languages (not both, not neither), and check whether English-only producers of that concept skew toward English-dominant exposure while Vietnamese-only producers skew toward Vietnamese-dominant exposure.

# All usable EN-VI concept matches (not just the 13 words from earlier) -- same
# candidate-filtering logic as sig_mono_crosswalk above, applied to the full table.
concept_match_full <- concept_match %>%
  mutate(en_candidates = map(en_candidates, ~ intersect(.x, en_data$d_items$definition))) %>%
  filter(lengths(en_candidates) > 0) %>%
  left_join(vi_data$d_items %>% select(definition, item_definition), by = c("vi_definition" = "definition"))

cat(nrow(concept_match_full), "concept pairs have a usable match for this analysis.\n")
## 639 concept pairs have a usable match for this analysis.
relative_exposure <- bilingual_full$pct_active_english - bilingual_full$pct_active_vietnamese  # + = English-dominant

# For one concept: classifies each Bilingual child as Both/Neither/EN_only/VI_only,
# based on that concept's Vietnamese column and (possibly multiple) English candidate
# columns (produces the concept in English if they produce ANY candidate).
classify_concept <- function(vi_col, en_cols) {
  vi_produces <- vi_bilingual_mat[, vi_col]
  en_sub <- en_bilingual_mat[, en_cols, drop = FALSE]
  en_produces <- if (ncol(en_sub) == 1) as.integer(en_sub[, 1]) else as.integer(apply(en_sub == 1, 1, any))
  case_when(
    vi_produces == 1 & en_produces == 1 ~ "Both",
    vi_produces == 0 & en_produces == 0 ~ "Neither",
    vi_produces == 0 & en_produces == 1 ~ "EN_only",
    vi_produces == 1 & en_produces == 0 ~ "VI_only"
  )
}

Per-concept test: do EN-only and VI-only producers differ in exposure balance?

min_n_per_group <- 4  # below this, a two-group comparison isn't meaningfully powered

competition_results <- pmap_dfr(
  list(concept_match_full$vi_definition, concept_match_full$en_candidates, concept_match_full$item_definition),
  function(vi_col, en_cols, item_definition) {
    category <- classify_concept(vi_col, en_cols)
    n_en <- sum(category == "EN_only"); n_vi <- sum(category == "VI_only")
    if (n_en < min_n_per_group || n_vi < min_n_per_group) {
      return(tibble(vi_definition = vi_col, item_definition = item_definition,
                     n_en_only = n_en, n_vi_only = n_vi, effect = NA_real_, p_value = NA_real_))
    }
    exp_en <- relative_exposure[category == "EN_only"]
    exp_vi <- relative_exposure[category == "VI_only"]
    test <- wilcox.test(exp_en, exp_vi)
    tibble(vi_definition = vi_col, item_definition = item_definition,
           n_en_only = n_en, n_vi_only = n_vi,
           effect = mean(exp_en) - mean(exp_vi),  # + = EN-only producers are more English-dominant, as competition predicts
           p_value = test$p.value)
  }
) %>%
  left_join(vi_data$d_items %>% select(definition, item_kind), by = c("vi_definition" = "definition")) %>%
  mutate(p_fdr = p.adjust(p_value, method = "BH"))

cat(sum(!is.na(competition_results$p_value)), "/", nrow(competition_results),
    "concepts had enough EN-only and VI-only producers (n >=", min_n_per_group, "each) to test.\n")
## 587 / 639 concepts had enough EN-only and VI-only producers (n >= 4 each) to test.
cat(sum(competition_results$p_fdr < 0.05, na.rm = TRUE), "concepts show a significant competition effect at FDR < .05.\n")
## 551 concepts show a significant competition effect at FDR < .05.
competition_results %>%
  filter(!is.na(p_value)) %>%
  arrange(p_fdr) %>%
  slice_head(n = 20) %>%
  select(item_definition, item_kind, n_en_only, n_vi_only, effect, p_value, p_fdr) %>%
  kable(digits = c(0, 0, 0, 0, 1, 4, 4),
        caption = "Strongest language-competition words: EN-only producers' exposure balance vs. VI-only producers' (positive effect = competition-consistent).") %>%
  html_table_width(c(130, 90, 70, 70, 70, 80, 80))
Strongest language-competition words: EN-only producers’ exposure balance vs. VI-only producers’ (positive effect = competition-consistent).
item_definition item_kind n_en_only n_vi_only effect p_value p_fdr
giấu action_words 16 16 84.7 0 0e+00
yêu action_words 16 19 83.9 0 0e+00
chuột animals 15 14 101.4 0 0e+00
gấu animals 20 11 95.9 0 0e+00
ngựa animals 21 18 81.7 0 0e+00
nặng descriptive_words 13 14 107.7 0 0e+00
trứng food_drink 25 13 89.4 0 0e+00
bát household 19 18 83.0 0 0e+00
gối household 17 20 86.4 0 0e+00
kính household 13 20 100.1 0 0e+00
rổ household 13 16 106.3 0 0e+00
bò animals 16 15 94.4 0 0e+00
rốn body_parts 22 16 80.6 0 0e+00
búp bê toys 14 23 81.2 0 0e+00
trèo action_words 8 19 109.7 0 0e+00
ngoài places 17 12 89.4 0 0e+00
thuyền vehicles 19 9 101.9 0 0e+00
giày clothing 20 15 87.4 0 1e-04
lưỡi body_parts 13 25 87.0 0 1e-04
quần clothing 10 25 101.5 0 1e-04 
top_competition <- competition_results %>% filter(!is.na(p_value)) %>% slice_min(p_fdr, n = 10)

bind_rows(
  top_competition %>% transmute(item_definition, value = effect, type = "Observed effect"),
) %>%
  mutate(item_definition = fct_reorder(item_definition, value)) %>%
  ggplot(aes(x = value, y = item_definition)) +
  geom_col(fill = "#6A3D9A") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey40") +
  theme_classic() +
  xlab("Mean relative exposure, EN-only producers minus VI-only producers\n(+ = consistent with competition)") +
  ylab(NULL) +
  labs(title = "Top 10 language-competition words")

Result: this isn’t a handful of special words – it’s nearly the whole vocabulary

The honest summary looks different from a short “top competing words” list. Of 587 testable concepts, 551 (94%) show a significant competition effect, and – more strikingly – not a single one of the 587 goes in the wrong direction (every single effect estimate is positive: EN-only producers are always more English-dominant on average than VI-only producers for that same concept, never the reverse). Effect sizes are large throughout (median 76 percentage points of relative exposure, on a scale that tops out at 200). This is a far more pervasive, far stronger pattern than “a few words show competition” – it looks like which language a bilingual child happens to produce a given concept in is overwhelmingly determined by their overall language dominance, for nearly any concept they know in only one language.

competition_results %>%
  filter(!is.na(p_value)) %>%
  arrange(desc(p_fdr)) %>%
  slice_head(n = 15) %>%
  select(item_definition, item_kind, n_en_only, n_vi_only, effect, p_fdr) %>%
  kable(digits = c(0, 0, 0, 0, 1, 3),
        caption = "The 15 WEAKEST competition words -- still positive-direction, just smaller/less reliable than the rest.") %>%
  html_table_width(c(130, 110, 70, 70, 70, 80)
  )
The 15 WEAKEST competition words – still positive-direction, just smaller/less reliable than the rest.
item_definition item_kind n_en_only n_vi_only effect p_fdr
meo sounds 6 10 13.0 0.800
không games_routines 24 12 9.5 0.697
be be sounds 35 5 18.3 0.452
máy kéo vehicles 15 5 35.3 0.216
chào games_routines 31 15 23.5 0.169
cô ấy pronouns 11 7 39.0 0.146
em bé people 14 13 33.1 0.136
của locations 4 24 53.8 0.131
ngựa con animals 10 13 39.7 0.127
hơn quantifiers 26 13 32.3 0.118
nên connecting_words 4 5 77.5 0.113
của chúng ta pronouns 7 4 45.4 0.111
cúc cù cu cu sounds 11 7 50.0 0.107
cây household 5 25 48.9 0.107
và connecting_words 11 11 38.2 0.104

Given that, the more interesting question flips: not “which words show competition” (nearly all do) but “which words show it least.” The weakest cases above are disproportionately high-frequency routine words (không/no, chào/hello, sound effects like meo/meow, be be/baa) – exactly the category identified earlier as least exposure-sensitive in general (robust to dilution because they come up constantly regardless of total exposure). That’s a coherent story: routine words are frequent enough in either language’s input that even a non-dominant- language child picks them up and produces them somewhat independent of overall dominance, slightly blurring the EN-only/VI-only split for exactly these items – while content words with no such frequency floor sort almost perfectly by language dominance.

Pooled check (important caveat on independence)

pooled_data <- pmap_dfr(
  list(concept_match_full$vi_definition, concept_match_full$en_candidates),
  function(vi_col, en_cols) {
    category <- classify_concept(vi_col, en_cols)
    keep <- category %in% c("EN_only", "VI_only")
    if (sum(keep) == 0) return(NULL)
    tibble(is_en_only = as.integer(category[keep] == "EN_only"), relative_exposure = relative_exposure[keep])
  }
)

cat("Pooled rows (every child x concept where they produced exactly one language's word):", nrow(pooled_data), "\n")
## Pooled rows (every child x concept where they produced exactly one language's word): 16953
summary(glm(is_en_only ~ relative_exposure, family = binomial, data = pooled_data))
## 
## Call:
## glm(formula = is_en_only ~ relative_exposure, family = binomial, 
##     data = pooled_data)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -0.0436024  0.0194677   -2.24   0.0251 *  
## relative_exposure  0.0292787  0.0004358   67.18   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 23451  on 16952  degrees of freedom
## Residual deviance: 16304  on 16951  degrees of freedom
## AIC: 16308
## 
## Number of Fisher Scoring iterations: 4

Caveat this pooled model does not fix: each Bilingual child contributes one row per concept they’re a single- language producer for, so the same ~101 children appear many times each – these rows are not independent observations, which makes the pooled model’s standard errors (and its p-value) too optimistic. It’s included only as a large-sample sanity check that the direction and rough size of the effect matches the per-concept results above, not as a properly-powered significance test on its own; the per-concept Wilcoxon tests (one independent set of children per concept, though still reusing many of the same ~101 children across different concepts) are the more defensible unit of evidence, and even those should be read as suggestive given how many concepts are being screened at a fairly small per-concept n.

RQ4: Does the competition effect depend on how frequency-robust a word is?

The weakest competition words identified above (không/no, chào/hello, meo/meow, be be/baa) were only inspected qualitatively. H7 predicts this formally: high-frequency, functionally routine words (produced through obligatory daily-routine input that both languages supply regardless of overall dominance) should show a smaller competition effect than discretionary content words, whose production depends more on which language gets more dedicated engagement time. Two operationalizations of “frequency-robust,” tested against competition_results’ effect column (already confirmed earlier: positive for every one of the 587 testable concepts):

  1. Continuous: IRT difficulty (b, from the Vietnamese 2PL model, already loaded as diff_vi_full above) – easier/more frequent words should show a smaller competition effect, so H7 predicts a positive correlation between b and effect (harder/rarer words show bigger dominance-driven separation).
  2. Categorical: closed-class/functional categories (games_routines, sounds, connecting_words, helping_verbs, quantifiers, pronouns, time_words – chosen a priori as the routine/functional word classes, not post-hoc) vs. open-class content categories (animals, food/drink, clothing, vehicles, toys, outdoor, household, furniture/rooms, body parts, people, places, descriptive words, action words, locations). H7 predicts functional categories show smaller effect than content categories.
competition_with_difficulty <- competition_results %>%
  filter(!is.na(effect)) %>%
  left_join(diff_vi_full %>% select(definition, b), by = c("vi_definition" = "definition"))

routine_categories <- c("games_routines", "sounds", "connecting_words", "helping_verbs", "quantifiers",
                         "pronouns", "time_words")
competition_with_difficulty <- competition_with_difficulty %>%
  mutate(word_class = if_else(item_kind %in% routine_categories, "Functional/routine", "Content"))

table(competition_with_difficulty$word_class)
## 
##            Content Functional/routine 
##                497                 90

Continuous test: difficulty vs. competition effect size

cor.test(competition_with_difficulty$b, competition_with_difficulty$effect, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  competition_with_difficulty$b and competition_with_difficulty$effect
## t = 5.4855, df = 585, p-value = 6.143e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1428115 0.2967958
## sample estimates:
##       cor 
## 0.2211818
ggplot(competition_with_difficulty, aes(x = b, y = effect)) +
  geom_point(alpha = 0.4, size = 1) +
  geom_smooth(method = "lm", se = TRUE, color = "#E31A1C") +
  theme_classic() +
  xlab("IRT difficulty (b) -- higher = harder/rarer") +
  ylab("Competition effect (relative-exposure separation)") +
  labs(title = "Does a harder/rarer word show a bigger competition effect?")
## `geom_smooth()` using formula = 'y ~ x'

Categorical test: functional/routine vs. content words

competition_with_difficulty %>%
  group_by(word_class) %>%
  summarise(n = n(), mean_effect = mean(effect), median_effect = median(effect), .groups = "drop") %>%
  kable(digits = 1, caption = "Competition effect size by word class.") %>%
  html_table_width(c(150, 60, 100, 100))
Competition effect size by word class.
word_class n mean_effect median_effect
Content 497 77.3 77.1
Functional/routine 90 69.3 70.2 
wilcox.test(effect ~ word_class, data = competition_with_difficulty)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  effect by word_class
## W = 26758, p-value = 0.003011
## alternative hypothesis: true location shift is not equal to 0
ggplot(competition_with_difficulty, aes(x = word_class, y = effect, fill = word_class)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.15, alpha = 0.3) +
  theme_classic() +
  theme(legend.position = "none") +
  scale_fill_manual(values = c(`Functional/routine` = "#1F78B4", Content = "#E31A1C")) +
  xlab(NULL) + ylab("Competition effect (relative-exposure separation)") +
  labs(title = "Functional/routine vs. content words")

If both tests come back in the predicted direction (positive difficulty correlation; functional words showing a smaller effect than content words), H7 is supported as a formal, quantitative claim rather than an impression from eyeballing the weakest cases – and the paper can report the correlation coefficient and group means directly rather than naming four illustrative words.

Caveats

  • Small monolingual group (n = 42): many individual words won’t reach significance simply from limited power, especially rare items. Read the ranked tables above (raw effect size) as suggestive/exploratory, and treat only the FDR-corrected subset as a defensible list of “distinguishing” words.
  • Group definition is production-based, not exposure-based: a “Monolingual” label here means “produced zero English words on this survey,” not “has no English input.” Some of these children may simply be earlier in expressive language development for English specifically.
  • Age is not a strong confound empirically (checked above: nearly identical age distributions between groups), but the age-adjusted glm_or_age_adj column is included per word as a check regardless – a word whose sign flips or becomes NA there is less trustworthy as a genuine group effect.
  • Multiple comparisons: ~680 words tested; FDR (Benjamini-Hochberg) correction is applied, but with this much multiplicity and a modest smaller-group size, some false discoveries among the FDR-significant set are still expected.