Motivation & Approach

Specific questions: Do E-V bilingual children show a noun/verb bias? If so, are they consistent across languages? - Does noun/verb dominance also change with age / other demographic characteristics? Possible big question: How should we characterize noun-verb dominance? Is noun-verb dominance driven by conceptual underpinnings (e.g., Mcnamara 1972) vs. language input features (Tardif 1996; also Gentner 1982)?

Alternatives

H1: Cross-linguistic consistency. Noun/verb bias is consistent across the two languages of E–V bilingual children. This pattern would be consistent with noun/verb dominance being underpinned primarily by language-independent conceptual differences.

H2: Cross-linguistic inconsistency. Noun/verb bias differs across the two languages of E–V bilingual children. This pattern would be consistent with noun/verb dominance being shaped primarily by language-specific input features. Importantly, cross-linguistic inconsistency could take at least two forms:

H2.1: Monolingual-like biases: Bilingual children show a noun/verb bias in each language that is comparable to the bias observed in monolingual children acquiring that language. This would suggest that each language’s bias is largely shaped by its own input.

H2.2: Language-specific biases with a lexical “accent”: Bilingual children show different noun/verb biases across their two languages, but each language’s bias is shifted relative to the corresponding monolingual pattern. This would suggest that language-specific input interacts with the bilingual child’s cross-linguistic lexical experience.

Note: We would need additional data from monolingual speakers to differentiate between H2.1 and H2.2.

Methodological Approach

  1. Fit bifactor IRT models that estimate kids’ noun ability & verb ability, in both English and Vietnamese. Check whether these models fit the data better than the corresponding unidimensional IRT model.
  2. Use estimated noun ability to bootstrap 95% CI noun ability for each kid, in each language.
  3. Check whether 95% CI EN-VI overlap, how much, does this change with age / language exposure for each language?

Methods Overview

The existing 2PL models (02_IRT_analysis_en.Rmd, 02_IRT_analysis_vi_0716.Rmd) estimate a single, unidimensional production-vocabulary ability per child. These were the best fitting models for English and Vietnamese vocabulary respectively (compared with unidimensional Rasch/3PL/4PL). Here we additionally estimate noun-specific and verb-specific ability using a confirmatory bifactor model: every item loads on a general factor (G), and items belonging to noun or verb categories additionally load on a noun-specific (S1) or verb-specific (S2) factor. Each child therefore gets three ability estimates per language:

  • ability_general (G): overall production vocabulary, net of noun/verb group membership
  • ability_noun (S1): noun ability beyond what is explained by general vocabulary ability
  • ability_verb (S2): verb ability beyond what is explained by general vocabulary ability

Item-to-domain mapping is based on item_kind (the CDI category embedded in each column name, e.g. action_words, animals, …):

  • verb = action_words
  • noun = animals, body_parts, clothing, food_drink, furniture_rooms, household, outdoor, people, places, toys, vehicles
  • all other item_kinds (descriptive_words, pronouns, quantifiers, helping_verbs, question_words, connecting_words, locations, games_routines, sounds, time_words) load on the general factor only.

Modeling notes

  • A plain 3-dimensional mirt() fit using rectangular quadrature does not scale to ~680 items and is extremely slow to converge. method = "QMCEM" (quasi-Monte Carlo integration) is mirt’s documented fix for 3+ dimensional models.
  • Without priors, item discriminations on the specific (noun/verb) factors blow up into Heywood cases (|a| > 50, NaN log-likelihood) given the current sample size (~100 usable children per language). Log-normal priors on a1/a2/a3 keep discriminations positive and bounded, which resolves this.
  • We use mirt() directly (not the bfactor() convenience wrapper) because bfactor()’s fast dimension-reduction EM algorithm does not respect the PRIOR statements needed to stabilize the fit at this sample size.
noun_kinds <- c("animals", "body_parts", "clothing", "food_drink", "furniture_rooms",
                "household", "outdoor", "people", "places", "toys", "vehicles")
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_",
    out_dir     = "models/en"
  ),
  vi = list(
    data_file   = "data/VietnameseWS_Bui_data.csv",
    fields_file = "data/VietnameseWS_Bui_fields.csv",
    prefix      = "cat_vie_",
    out_dir     = "models/vi"
  )
)
# For a demo data frame with language_1..4 and lang_1_pct_active..lang_4_pct_active
# (lang_i_pct_active is the % of the child's *active* language use taken up by
# whichever language sits in slot i -- these sum to ~100 across the 4 slots),
# returns each child's % active use of `lang_name`, or 0 if `lang_name` is not
# among their reported languages.
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)
}
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, language_1, ppl_1_edu,
           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"),
      pct_active_vietnamese = get_pct_active_lang(., "Vietnamese")
    )

  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, "^.*?___"),
    ) |>
    rename(item_id = field) |>
    select(-c(group, type))

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

en_data <- load_lang_data("en")
d_wide_en  <- en_data$d_wide
d_demo_en  <- en_data$d_demo
d_items_en <- en_data$d_items

vi_data <- load_lang_data("vi")
d_wide_vi  <- vi_data$d_wide
d_demo_vi  <- vi_data$d_demo
d_items_vi <- vi_data$d_items

We are currently not filtering any data, but this function is written here in case we want to filter kids who do not produce any English words etc.

# to_remove would drop children too young to be producing words (<12mo),
# who produce zero words at >=12mo, or who are at ceiling (produce every
# item) -- but that exclusion is currently disabled below. In particular,
# the 42 English zero-producers still produce substantial Vietnamese
# vocabulary (median ~244 words), so they're kept rather than treated as
# non-responders. to_remove is still computed for reference/diagnostics.
filter_d_mat <- function(d_wide, d_demo, d_items) {
  d_mat <- d_wide %>% data.frame %>% select(-response_id) %>% data.matrix
  stopifnot(all(colnames(d_mat) == d_items$definition))

  d_demo$production <- rowSums(d_mat)
  too_young <- which(d_demo$age < 12)
  no_words  <- which(d_demo$production == 0 & d_demo$age >= 12)
  ceiling   <- which(d_demo$production == ncol(d_mat))
  to_remove <- unique(c(too_young, no_words, ceiling))
  # if (length(to_remove) > 0) {
  #   d_mat  <- d_mat[-to_remove, ]
  #   d_demo <- d_demo[-to_remove, ]
  # }
  list(d_mat = d_mat, d_demo = d_demo)
}

en_filtered <- filter_d_mat(d_wide_en, d_demo_en, d_items_en)
d_mat_en  <- en_filtered$d_mat
d_demo_en <- en_filtered$d_demo

vi_filtered <- filter_d_mat(d_wide_vi, d_demo_vi, d_items_vi)
d_mat_vi  <- vi_filtered$d_mat
d_demo_vi <- vi_filtered$d_demo
# Fits the same G/S1(noun)/S2(verb) bifactor structure at a given itemtype
# (1PL/2PL/3PL/4PL), all treated the same way, so the four variants can be
# compared the way the original 02_IRT_analysis Rmds compare unidimensional
# 1PL/2PL/3PL/4PL.
#
# itemtype = "Rasch" (1PL) fixes all discriminations to 1, so there is
# nothing to put a prior on for a1/a2/a3 (only "2PL"/"3PL"/"4PL" freely
# estimate discriminations, and need lnorm priors to avoid Heywood cases --
# see the modeling notes above).

fit_bifactor_variant <- function(d_mat, d_items, out_dir, itemtype, ncycles = 4000, seed = 1234) {
  noun_idx <- which(d_items$item_kind %in% noun_kinds)
  verb_idx <- which(d_items$item_kind %in% verb_kinds)
  n_items  <- ncol(d_mat)

  model_lines <- c(
    paste0("G = 1-", n_items),
    paste0("S1 = ", paste(noun_idx, collapse = ",")),
    paste0("S2 = ", paste(verb_idx, collapse = ","))
  )
  if (itemtype != "Rasch") {
    model_lines <- c(model_lines, paste0(
      "PRIOR = (1-", n_items, ", a1, lnorm, 0, 0.5), (",
      paste(noun_idx, collapse = ","), ", a2, lnorm, 0, 0.5), (",
      paste(verb_idx, collapse = ","), ", a3, lnorm, 0, 0.5)"
    ))
  }
  mod_def <- mirt.model(paste(model_lines, collapse = "\n"))

  set.seed(seed)
  mod <- mirt(d_mat, mod_def, itemtype = itemtype, method = "QMCEM",
              verbose = TRUE, technical = list(NCYCLES = ncycles))

  dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)
  label <- c(Rasch = "1pl", `2PL` = "2pl", `3PL` = "3pl", `4PL` = "4pl")[[itemtype]]
  saveRDS(mod, file = file.path(out_dir, paste0("mod_bifactor_", label, ".Rds")))
  mod
}

# builds one fit-statistics row per model (AIC/SABIC/HQ/BIC/logLik) so the
# four bifactor variants can be compared side by side, in addition to the
# pairwise nested LRTs from get_anova_table() in IRT_helpers.R
get_fit_table <- function(models, model_names) {
  tab <- dplyr::bind_rows(lapply(models, function(m) as.data.frame(anova(m))))
  tab <- data.frame(Model = model_names, tab, row.names = NULL)
  tab
}

# item coefficients (with domain labels) and per-child G/S1(noun)/S2(verb)
# ability estimates (+ SE) for a fitted bifactor model -- used on whichever
# variant (typically 2PL) is chosen as the primary model for the rest of
# this report.
extract_bifactor_summary <- function(mod, d_demo, d_items) {
  coefs <- as_tibble(coef(mod, simplify = TRUE)$items) %>%
    mutate(definition = rownames(coef(mod, simplify = TRUE)$items)) %>%
    left_join(d_items, by = "definition") %>%
    mutate(domain = case_when(
      item_kind %in% verb_kinds ~ "verb",
      item_kind %in% noun_kinds ~ "noun",
      TRUE ~ "other"
    ))

  fs <- fscores(mod, method = "MAP", full.scores.SE = TRUE)
  fscores_tab <- tibble(
    row_id          = d_demo$row_id,
    response_id     = d_demo$response_id,
    ability_general = fs[, "G"],
    ability_noun    = fs[, "S1"],
    ability_verb    = fs[, "S2"],
    ability_general_SE = fs[, "SE_G"],
    ability_noun_SE    = fs[, "SE_S1"],
    ability_verb_SE    = fs[, "SE_S2"]
  )

  list(coefs = coefs, fscores = fscores_tab)
}

# Nonparametric item-resampling bootstrap CI for ability_noun (S1), holding the
# already-fitted item parameters fixed. Refitting the full 3-dimensional
# bifactor model per replicate isn't feasible -- a single fit already takes
# 40-90+ minutes -- so instead, for each child and each of B replicates, items
# are resampled with replacement (equivalently: re-weighted by how many times
# each item is drawn) and that child's (G, S1, S2) MAP is re-estimated under
# the model's fixed, already-estimated item parameters via a small weighted
# penalized-logistic optimization; the 2.5th/97.5th percentiles of the
# resulting S1 draws become the per-child CI.
# cache_file, if given, is read (skipping the bootstrap entirely) when it
# already exists, and written after computing otherwise -- same on-disk
# caching convention as fit_bifactor_variant()'s saveRDS().
#
# IMPORTANT: the prior variances on (G, S1, S2) are NOT fixed at 1 in general
# -- e.g. the Rasch/1PL bifactor variant fixes all item discriminations to 1
# for identification, which forces the *latent variance* to freely absorb the
# item-difficulty spread instead (checked: English's fitted var(G) = 24.55,
# var(S1) = 0.87, var(S2) = 0.58; Vietnamese's differ again). An earlier
# version of this function assumed independent N(0,1) priors, which silently
# used the wrong (much tighter) prior on G in particular, over-shrinking G
# and biasing S1 to compensate. Pulling the actual fitted variances from
# coef(mod, simplify=TRUE)$cov (confirmed diagonal / orthogonal for these
# models) fixes this.
bootstrap_noun_ci <- function(mod, d_mat, d_demo, B = 1000, seed = 1234, cache_file = NULL) {
  if (!is.null(cache_file) && file.exists(cache_file)) {
    return(readRDS(cache_file))
  }

  cf <- coef(mod, simplify = TRUE)$items
  a1 <- cf[, "a1"]; a2 <- cf[, "a2"]; a3 <- cf[, "a3"]
  d_par <- cf[, "d"]; g <- cf[, "g"]; u <- cf[, "u"]
  n_items <- nrow(cf)

  prior_var <- diag(coef(mod, simplify = TRUE)$cov)[c("G", "S1", "S2")]

  neg_log_post <- function(theta, y, w) {
    z <- a1 * theta[1] + a2 * theta[2] + a3 * theta[3] + d_par
    p <- g + (u - g) / (1 + exp(-z))
    p <- pmin(pmax(p, 1e-10), 1 - 1e-10)
    -(sum(w * (y * log(p) + (1 - y) * log(1 - p))) +
        sum(dnorm(theta, mean = 0, sd = sqrt(prior_var), log = TRUE)))
  }
  gr_neg_log_post <- function(theta, y, w) {
    z <- a1 * theta[1] + a2 * theta[2] + a3 * theta[3] + d_par
    ez <- 1 / (1 + exp(-z))
    p  <- g + (u - g) * ez
    p  <- pmin(pmax(p, 1e-10), 1 - 1e-10)
    dp_dz  <- (u - g) * ez * (1 - ez)
    dll_dz <- w * (y / p - (1 - y) / (1 - p)) * dp_dz
    -c(sum(dll_dz * a1) - theta[1] / prior_var[1],
       sum(dll_dz * a2) - theta[2] / prior_var[2],
       sum(dll_dz * a3) - theta[3] / prior_var[3])
  }

  fs <- fscores(mod, method = "MAP")
  n_people <- nrow(d_mat)
  set.seed(seed)
  boot_S1 <- matrix(NA_real_, n_people, B)
  for (i in seq_len(n_people)) {
    y <- d_mat[i, ]
    start <- fs[i, ]
    for (b in seq_len(B)) {
      w <- tabulate(sample.int(n_items, n_items, replace = TRUE), nbins = n_items)
      boot_S1[i, b] <- optim(start, neg_log_post, gr_neg_log_post, y = y, w = w, method = "BFGS")$par[2]
    }
  }

  out <- tibble(
    row_id                 = d_demo$row_id,
    response_id            = d_demo$response_id,
    ability_noun_boot_mean = rowMeans(boot_S1),
    ability_noun_boot_SD   = apply(boot_S1, 1, sd),
    ability_noun_CI_lower  = apply(boot_S1, 1, quantile, probs = 0.025),
    ability_noun_CI_upper  = apply(boot_S1, 1, quantile, probs = 0.975)
  )

  if (!is.null(cache_file)) saveRDS(out, cache_file)
  out
}

English

en_bifactor_1pl <- fit_bifactor_variant(d_mat_en, d_items_en, out_dir = "models/en", itemtype = "Rasch")
en_bifactor_2pl <- fit_bifactor_variant(d_mat_en, d_items_en, out_dir = "models/en", itemtype = "2PL")
en_bifactor_3pl <- fit_bifactor_variant(d_mat_en, d_items_en, out_dir = "models/en", itemtype = "3PL")
en_bifactor_4pl <- fit_bifactor_variant(d_mat_en, d_items_en, out_dir = "models/en", itemtype = "4PL")
load("models/en/mod_2pl.Rds")  # unidimensional baseline
mod_bifactor_1pl <- readRDS("models/en/mod_bifactor_1pl.Rds")
mod_bifactor_2pl <- readRDS("models/en/mod_bifactor_2pl.Rds")
mod_bifactor_3pl <- readRDS("models/en/mod_bifactor_3pl.Rds")
mod_bifactor_4pl <- readRDS("models/en/mod_bifactor_4pl.Rds")

Model comparison: Bifactor 1PL/2PL/3PL/4PL

Compares the bifactor structure (G + noun + verb factors) fit at each itemtype, plus the plain unidimensional 2PL (mod_2pl) as a baseline.

mc1_en <- get_anova_table(mod_bifactor_1pl, mod_bifactor_2pl, c("Bifactor 1PL", "Bifactor 2PL"))
mc2_en <- get_anova_table(mod_bifactor_1pl, mod_bifactor_3pl, c("Bifactor 1PL", "Bifactor 3PL"))
mc3_en <- get_anova_table(mod_bifactor_1pl, mod_bifactor_4pl, c("Bifactor 1PL", "Bifactor 4PL"))
mc4_en <- get_anova_table(mod_2pl, mod_bifactor_1pl, c("Unidim 2PL", "Bifactor 1PL"))

kable(rbind(mc1_en, mc2_en, mc3_en, mc4_en), digits = 2,
      caption = "English: pairwise nested-model comparisons.") %>%
  html_table_width(c(90, 90, 90, 90, 50))
English: pairwise nested-model comparisons.
Model AIC BIC logLik df
Bifactor 1PL 39314.19 41335.98 -18973.10 NA
Bifactor 2PL 40909.21 46315.41 -18625.60 1145
Bifactor 1PL 39314.19 41335.98 -18973.10 NA
Bifactor 3PL 42201.92 49621.05 -18590.96 1826
Bifactor 1PL 39314.19 41335.98 -18973.10 NA
Bifactor 4PL 43372.01 52804.06 -18495.01 2507
Unidim 2PL 40365.38 44391.22 -18820.69 NA
Bifactor 1PL 39314.19 41335.98 -18973.10 -678 
bifactor_variants_en <- list(`1PL` = mod_bifactor_1pl, `2PL` = mod_bifactor_2pl, `3PL` = mod_bifactor_3pl, `4PL` = mod_bifactor_4pl, `Unidim2PL` = mod_2pl)
fit_table_en <- get_fit_table(bifactor_variants_en, names(bifactor_variants_en))
kable(fit_table_en, digits = 2,
      caption = "English: fit statistics across the four bifactor itemtypes.") %>%
  html_table_width(c(90, 90, 90, 90, 90))
English: fit statistics across the four bifactor itemtypes.
Model AIC SABIC HQ BIC logLik logPost
1PL 39314.19 39171.76 40135.76 41335.98 -18973.10 NA
2PL 40909.21 40528.34 43106.07 46315.41 -18625.60 -19303.36
3PL 42201.92 41679.24 45216.75 49621.05 -18590.96 -19291.79
4PL 43372.01 42707.52 47204.81 52804.06 -18495.01 -19183.33
Unidim2PL 40365.38 40081.76 42001.32 44391.22 -18820.69 NA

Model selection

Selects the best-fitting bifactor variant by BIC. All of the inferential analysis below (noun/verb dominance, age effects, etc.) is run on this selected model, referred to from here on as mod_bifactor.

The best model is the Bifactor 1PL.

best_variant_en <- names(bifactor_variants_en)[which.min(fit_table_en$BIC)]
mod_bifactor <- bifactor_variants_en[[best_variant_en]]
cat("Best-fitting English bifactor variant (by BIC):", best_variant_en)
## Best-fitting English bifactor variant (by BIC): 1PL
en_bifactor_summary <- extract_bifactor_summary(mod_bifactor, d_demo_en, d_items_en)
en_coefs <- en_bifactor_summary$coefs
en_fscores <- en_bifactor_summary$fscores
cat("Converged:", extract.mirt(mod_bifactor, "converged"),
    "| Iterations:", extract.mirt(mod_bifactor, "iterations"),
    "| logLik:", round(logLik(mod_bifactor), 1))
## Converged: TRUE | Iterations: 115 | logLik: -18973.1

Investigate ability_noun and ability_verb for spread. Note that these estimates are not on the same scale.

summary(en_fscores$ability_noun)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1.7579 -0.4934 -0.1396 -0.1017  0.1870  4.4349
summary(en_fscores$ability_verb)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -1.573131 -0.279244 -0.007532 -0.018926  0.068660  2.350589
sd(en_fscores$ability_noun)
## [1] 0.7762871
sd(en_fscores$ability_verb)
## [1] 0.5171733

Step 2: 95% CI on noun ability (item-resampling bootstrap)

bootstrap_noun_ci() (defined above) gives each child a nonparametric 95% bootstrap CI on ability_noun, built by resampling items with replacement (B = 1000 replicates) and re-estimating that child’s (ability-general, ability-noun, ability-verb) MAP under the fixed, already-fitted item parameters – see the function definition for why a full model-refit bootstrap isn’t used. Results (the point estimates and CI bounds, not the raw B draws) are cached to models/en/noun_boot_ci_<variant>.Rds. We also attach demographics (age, % active English/Vietnamese exposure) here for use in the cross-language comparison below.

en_noun_boot <- bootstrap_noun_ci(mod_bifactor, d_mat_en, d_demo_en, B = 1000,
                                   cache_file = file.path("models/en", paste0("noun_boot_ci_", best_variant_en, ".Rds")))
en_fscores <- en_fscores %>% left_join(en_noun_boot, by = c("row_id", "response_id"))

noun_idx_en <- which(d_items_en$item_kind %in% noun_kinds)
noun_raw_en <- rowSums(d_mat_en[, noun_idx_en])

en_fscores_demog <- en_fscores %>%
  left_join(d_demo_en %>% select(row_id, response_id, age, pct_active_english, pct_active_vietnamese, production),
            by = c("row_id", "response_id")) %>%
  mutate(noun_raw = noun_raw_en, floor_noun = noun_raw == 0)

summary(en_fscores$ability_noun_boot_SD)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.01364 0.01463 0.22220 0.20442 0.30295 0.45878
summary(en_fscores$ability_noun_CI_upper - en_fscores$ability_noun_CI_lower)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.05252 0.05732 0.86037 0.79391 1.18746 1.85418

ability_noun and ability_verb share a nominal prior scale (the bifactor model fixes all three factor variances to 1), but their empirical spread differs. Each specific factor is measured by different items, and there are more noun items than verb items, so MAP shrinkage compresses them to different degrees (e.g. sd(ability_verb) is meaningfully smaller than sd(ability_noun) in this sample).

ggplot(aes(x = ability_noun, y = ability_verb),
       data = en_fscores) +
  geom_point() +
  geom_abline(yintercept = 0, slope = 1, linetype = "dashed")
## Warning in geom_abline(yintercept = 0, slope = 1, linetype = "dashed"):
## Ignoring unknown parameters: `yintercept`

Double-checking that nouns only load onto nouns + general, verbs only load onto verbs + general, others only load onto general.

en_coefs %>%
  group_by(domain) %>%
  summarise(n_items = n(), mean_a1 = mean(a1), mean_a2 = mean(a2), mean_a3 = mean(a3)) %>%
  kable(digits = 2, caption = "English: mean discrimination on G/S1(noun)/S2(verb) by domain.") %>%
  html_table_width(rep(120, 5))
English: mean discrimination on G/S1(noun)/S2(verb) by domain.
domain n_items mean_a1 mean_a2 mean_a3
noun 364 1 1 0
other 214 1 0 0
verb 103 1 0 1 
en_fscores %>%
  select(response_id, ability_general, ability_noun, ability_verb) %>%
  pivot_longer(-response_id, names_to = "dimension", values_to = "theta") %>%
  ggplot(aes(x = theta, fill = dimension)) +
  geom_histogram(alpha = 0.7, bins = 30, position = "identity") +
  facet_wrap(~dimension) +
  theme_classic() +
  theme(legend.position = "none") +
  xlab("Ability (theta)")

Vietnamese

vi_bifactor_1pl <- fit_bifactor_variant(d_mat_vi, d_items_vi, out_dir = "models/vi", itemtype = "Rasch")
vi_bifactor_2pl <- fit_bifactor_variant(d_mat_vi, d_items_vi, out_dir = "models/vi", itemtype = "2PL")
vi_bifactor_3pl <- fit_bifactor_variant(d_mat_vi, d_items_vi, out_dir = "models/vi", itemtype = "3PL")
vi_bifactor_4pl <- fit_bifactor_variant(d_mat_vi, d_items_vi, out_dir = "models/vi", itemtype = "4PL")
load("models/vi/mod_2pl.Rds")  # unidimensional baseline
mod_bifactor_1pl <- readRDS("models/vi/mod_bifactor_1pl.Rds")
mod_bifactor_2pl <- readRDS("models/vi/mod_bifactor_2pl.Rds")
mod_bifactor_3pl <- readRDS("models/vi/mod_bifactor_3pl.Rds")
mod_bifactor_4pl <- readRDS("models/vi/mod_bifactor_4pl.Rds")

Model comparison: Bifactor 1PL/2PL/3PL/4PL

Compares the bifactor structure (G + noun + verb factors) fit at each itemtype, plus the plain unidimensional 2PL (mod_2pl) as a baseline.

mc1_vi <- get_anova_table(mod_bifactor_1pl, mod_bifactor_2pl, c("Bifactor 1PL", "Bifactor 2PL"))
mc2_vi <- get_anova_table(mod_bifactor_1pl, mod_bifactor_3pl, c("Bifactor 1PL", "Bifactor 3PL"))
mc3_vi <- get_anova_table(mod_bifactor_1pl, mod_bifactor_4pl, c("Bifactor 1PL", "Bifactor 4PL"))
mc4_vi <- get_anova_table(mod_2pl, mod_bifactor_1pl, c("Unidim 2PL", "Bifactor 1PL"))

kable(rbind(mc1_vi, mc2_vi, mc3_vi, mc4_vi), digits = 2,
      caption = "Vietnamese: pairwise nested-model comparisons.") %>%
  html_table_width(c(90, 90, 90, 90, 50))
Vietnamese: pairwise nested-model comparisons.
Model AIC BIC logLik df
Bifactor 1PL 62866.50 64906.02 -30743.25 NA
Bifactor 2PL 62597.43 68098.22 -29437.71 1171
Bifactor 1PL 62866.50 64906.02 -30743.25 NA
Bifactor 3PL 63875.66 71407.11 -29389.83 1858
Bifactor 1PL 62866.50 64906.02 -30743.25 NA
Bifactor 4PL 65000.88 74562.98 -29265.44 2545
Unidim 2PL 64254.78 68316.09 -30753.39 NA
Bifactor 1PL 62866.50 64906.02 -30743.25 -684 
bifactor_variants_vi <- list(`1PL` = mod_bifactor_1pl, `2PL` = mod_bifactor_2pl,
                             `3PL` = mod_bifactor_3pl, `4PL` = mod_bifactor_4pl,
                             `Unidim2PL` = mod_2pl
                             )
fit_table_vi <- get_fit_table(bifactor_variants_vi, names(bifactor_variants_vi))
kable(fit_table_vi, digits = 2,
      caption = "Vietnamese: fit statistics across the four bifactor itemtypes.") %>%
  html_table_width(c(90, 90, 90, 90, 90))
Vietnamese: fit statistics across the four bifactor itemtypes.
Model AIC SABIC HQ BIC logLik logPost
1PL 62866.50 62722.81 63695.28 64906.02 -30743.25 NA
2PL 62597.43 62209.89 64832.73 68098.22 -29437.71 -30081.22
3PL 63875.66 63345.07 66936.14 71407.11 -29389.83 -30021.71
4PL 65000.88 64327.22 68886.53 74562.98 -29265.44 -29905.58
Unidim2PL 64254.78 63968.66 65905.13 68316.09 -30753.39 NA

Model selection

Selects the best-fitting bifactor variant by BIC (see the English section above for the rationale). All of the inferential analysis below is run on this selected model, referred to from here on as mod_bifactor.

The best model is Bifactor 1PL (same as English).

best_variant_vi <- names(bifactor_variants_vi)[which.min(fit_table_vi$BIC)]
mod_bifactor <- bifactor_variants_vi[[best_variant_vi]]
cat("Best-fitting Vietnamese bifactor variant (by BIC):", best_variant_vi)
## Best-fitting Vietnamese bifactor variant (by BIC): 1PL
vi_bifactor_summary <- extract_bifactor_summary(mod_bifactor, d_demo_vi, d_items_vi)
vi_coefs <- vi_bifactor_summary$coefs
vi_fscores <- vi_bifactor_summary$fscores
cat("Converged:", extract.mirt(mod_bifactor, "converged"),
    "| Iterations:", extract.mirt(mod_bifactor, "iterations"),
    "| logLik:", round(logLik(mod_bifactor), 1))
## Converged: TRUE | Iterations: 62 | logLik: -30743.2

Investigate ability_noun and ability_verb for spread.

summary(vi_fscores$ability_noun)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -2.59437 -0.51041 -0.07837 -0.04878  0.35845  2.06732
summary(vi_fscores$ability_verb)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -2.56118 -0.74693 -0.22793 -0.02973  0.55725  4.02561
sd(vi_fscores$ability_noun)
## [1] 0.7596996
sd(vi_fscores$ability_verb)
## [1] 1.053916

Step 2: 95% CI on noun ability (item-resampling bootstrap)

Same procedure as the English section above (bootstrap_noun_ci(), defined earlier): a nonparametric 95% bootstrap CI on ability_noun per child, from B = 1000 item-resampling replicates under the fixed, already-fitted item parameters.

vi_noun_boot <- bootstrap_noun_ci(mod_bifactor, d_mat_vi, d_demo_vi, B = 1000,
                                   cache_file = file.path("models/vi", paste0("noun_boot_ci_", best_variant_vi, ".Rds")))
vi_fscores <- vi_fscores %>% left_join(vi_noun_boot, by = c("row_id", "response_id"))

noun_idx_vi <- which(d_items_vi$item_kind %in% noun_kinds)
noun_raw_vi <- rowSums(d_mat_vi[, noun_idx_vi])

vi_fscores_demog <- vi_fscores %>%
  left_join(d_demo_vi %>% select(row_id, response_id, age, pct_active_english, pct_active_vietnamese, production),
            by = c("row_id", "response_id")) %>%
  mutate(noun_raw = noun_raw_vi, floor_noun = noun_raw == 0)

summary(vi_fscores$ability_noun_boot_SD)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1660  0.2034  0.2399  0.2599  0.2982  0.3972
summary(vi_fscores$ability_noun_CI_upper - vi_fscores$ability_noun_CI_lower)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.6578  0.7995  0.9237  1.0141  1.1816  1.5812

ability_noun and ability_verb share a nominal prior scale but differ in empirical spread due to differential shrinkage (see the English section above for why).

ggplot(aes(x = ability_noun, y = ability_verb),
       data = vi_fscores) +
  geom_point() +
  geom_abline(yintercept = 0, slope = 1)
## Warning in geom_abline(yintercept = 0, slope = 1): Ignoring unknown parameters:
## `yintercept`

vi_coefs %>%
  group_by(domain) %>%
  summarise(n_items = n(), mean_a1 = mean(a1), mean_a2 = mean(a2), mean_a3 = mean(a3)) %>%
  kable(digits = 2, caption = "Vietnamese: mean discrimination on G/S1(noun)/S2(verb) by domain.") %>%
  html_table_width(rep(120, 5))
Vietnamese: mean discrimination on G/S1(noun)/S2(verb) by domain.
domain n_items mean_a1 mean_a2 mean_a3
noun 384 1 1 0
other 200 1 0 0
verb 103 1 0 1 
vi_fscores %>%
  select(response_id, ability_general, ability_noun, ability_verb) %>%
  pivot_longer(-response_id, names_to = "dimension", values_to = "theta") %>%
  ggplot(aes(x = theta, fill = dimension)) +
  geom_histogram(alpha = 0.7, bins = 30, position = "identity") +
  facet_wrap(~dimension) +
  theme_classic() +
  theme(legend.position = "none") +
  xlab("Ability (theta)")

Cross-language comparison

row_id is a child-level identifier shared by the English and Vietnamese surveys (unlike response_id, which is per-survey-response), confirmed to fully overlap across the two raw data files, and age is identical across the two surveys for the same child. It’s used here to join each child’s independently-derived English and Vietnamese noun-ability bootstrap CIs (Step 2 above), to ask whether a child’s noun ability is consistent across their two languages.

Step 3: Do the English and Vietnamese noun-ability 95% CIs overlap?

The English and Vietnamese bifactor models are fit independently, so their ability_noun scales aren’t directly comparable and have different spreads. So before computing overlap, ability_noun and its CI bounds are z-scored within each language, using that language’s own full-sample mean/SD.

For each child with both an English and a Vietnamese noun-ability estimate, overlap_amount is the signed length of the intersection of their two (standardized) 95% CIs: min(upper_en, upper_vi) - max(lower_en, lower_vi). A positive value is the width of CI overlap; a negative value is the size of the gap between two non-overlapping intervals. overlap is whether overlap_amount > 0.

en_mean_noun <- mean(en_fscores$ability_noun)
en_sd_noun   <- sd(en_fscores$ability_noun)
vi_mean_noun <- mean(vi_fscores$ability_noun)
vi_sd_noun   <- sd(vi_fscores$ability_noun)

cross_lang <- inner_join(
  en_fscores_demog %>%
    select(row_id, ability_noun_en = ability_noun,
           CI_lower_en = ability_noun_CI_lower, CI_upper_en = ability_noun_CI_upper,
           floor_noun_en = floor_noun, production_en = production,
           age, pct_active_english, pct_active_vietnamese),
  vi_fscores_demog %>%
    select(row_id, ability_noun_vi = ability_noun,
           CI_lower_vi = ability_noun_CI_lower, CI_upper_vi = ability_noun_CI_upper,
           floor_noun_vi = floor_noun, production_vi = production),
  by = "row_id"
) %>%
  mutate(
    ability_noun_en = (ability_noun_en - en_mean_noun) / en_sd_noun,
    CI_lower_en     = (CI_lower_en - en_mean_noun) / en_sd_noun,
    CI_upper_en     = (CI_upper_en - en_mean_noun) / en_sd_noun,
    ability_noun_vi = (ability_noun_vi - vi_mean_noun) / vi_sd_noun,
    CI_lower_vi     = (CI_lower_vi - vi_mean_noun) / vi_sd_noun,
    CI_upper_vi     = (CI_upper_vi - vi_mean_noun) / vi_sd_noun,
    overlap_amount  = pmin(CI_upper_en, CI_upper_vi) - pmax(CI_lower_en, CI_lower_vi),
    overlap         = overlap_amount > 0,
    gap_direction   = case_when(
      overlap_amount >= 0        ~ "Overlap",
      ability_noun_en > ability_noun_vi ~ "Gap: EN > VI",
      TRUE                        ~ "Gap: VI > EN"
    )
  )
cat(round(100 * mean(cross_lang$overlap), 1), "% of children have overlapping EN/VI 95% CIs on ability_noun\n")
## 67.6 % of children have overlapping EN/VI 95% CIs on ability_noun
summary(cross_lang$overlap_amount)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -4.3944 -0.1779  0.3336  0.2588  0.8983  1.9681

Gap cases (overlap_amount < 0, non-overlapping CIs) are colored by which language is higher, overlaps are greys.

ggplot(cross_lang, aes(x = overlap_amount, fill = gap_direction)) +
  geom_histogram(bins = 30) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  theme_classic() +
  scale_fill_manual(values = c(Overlap = "grey70", `Gap: EN > VI` = "#F8766D", `Gap: VI > EN` = "#00BFC4")) +
  labs(x = "EN-VI noun-ability CI overlap (ability z-scored)", fill = NULL)

42 of 142 children produced zero total words in at least one language (all English non-producers). The estimates for these children might not be reliable. Results when they are excluded below.

cross_lang_nonzero <- cross_lang %>% filter(production_en > 0, production_vi > 0)

cat(round(100 * mean(cross_lang_nonzero$overlap), 1),
    "% of children have overlapping EN/VI 95% CIs on ability_noun (excluding zero-producers)\n")
## 81 % of children have overlapping EN/VI 95% CIs on ability_noun (excluding zero-producers)
summary(cross_lang_nonzero$overlap_amount)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -4.3944  0.2810  0.7227  0.5464  1.0067  1.9681
ggplot(cross_lang_nonzero, aes(x = overlap_amount, fill = gap_direction)) +
  geom_histogram(bins = 30) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  theme_classic() +
  scale_fill_manual(values = c(Overlap = "grey70", `Gap: EN > VI` = "#F8766D", `Gap: VI > EN` = "#00BFC4")) +
  labs(x = "EN-VI noun-ability CI overlap (ability z-scored)", fill = NULL,
       title = "Excluding children who produced 0 words in either language")

One row per child, showing both language-specific 95% CIs on ability_noun side by side. Children are sorted by overlap_amount (ascending), so the least overlapping/most-gapped pairs are at the bottom and the most-overlapping pairs are at the top. Floor-score children (0 of that language’s noun items produced) are marked with a star instead of a dot. These are the cases where the estimates (and bootstrap CIs) are not trustworthy. These kids also make up most of the cases where CIs do not include the estimates, because for these kids most of the sampled replicates would

cross_lang_plot_df <- cross_lang %>%
  arrange(overlap_amount) %>%
  mutate(child_rank = row_number()) %>%
  select(child_rank, row_id, ability_noun_en, CI_lower_en, CI_upper_en, floor_noun_en,
         ability_noun_vi, CI_lower_vi, CI_upper_vi, floor_noun_vi) %>%
  pivot_longer(
    cols = -c(child_rank, row_id),
    names_to = c(".value", "lang"),
    names_pattern = "(ability_noun|CI_lower|CI_upper|floor_noun)_(en|vi)"
  ) %>%
  mutate(
    lang  = recode(lang, en = "English", vi = "Vietnamese"),
    y_pos = child_rank + ifelse(lang == "English", 0.2, -0.2),
    floor_noun = factor(floor_noun, levels = c(FALSE, TRUE), labels = c("Not floor", "Floor score (0 nouns)"))
  )

ggplot(cross_lang_plot_df, aes(y = y_pos, x = ability_noun, color = lang)) +
  geom_errorbar(aes(xmin = CI_lower, xmax = CI_upper), orientation = "y", width = 0, linewidth = 0.3) +
  geom_point(aes(shape = floor_noun, size = floor_noun)) +
  scale_shape_manual(values = c(`Not floor` = 16, `Floor score (0 nouns)` = 8)) +
  scale_size_manual(values = c(`Not floor` = 0.6, `Floor score (0 nouns)` = 1.6)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
  theme_classic() +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
  labs(x = "z(ability_noun), within-language (95% bootstrap CI)",
       y = "Children (sorted by EN-VI CI overlap, ascending)",
       color = "Language", shape = "Score type", size = "Score type",
       title = "Per-child 95% CI on noun ability: English vs. Vietnamese")

Does CI overlap change with age or language exposure?

pct_active_english and pct_active_vietnamese are substantially negatively correlated for these bilingual children (most report only these two active languages, so more of one tends to mean less of the other) – see the correlation below. To avoid collinearity in one regression, we lead with pct_active_vietnamese (the heritage-language exposure) as the primary predictor; swap in pct_active_english to check the mirror-image model.

cat("cor(pct_active_english, pct_active_vietnamese):",
    round(cor(cross_lang$pct_active_english, cross_lang$pct_active_vietnamese, use = "complete.obs"), 3))
## cor(pct_active_english, pct_active_vietnamese): -0.639

Older kids are more likely to show divergence between noun dominance in the two languages. I added pct_active_vietnamese as a way to “control” for language exposure, but not sure if that’s the correct approach.

summary(lm(overlap_amount ~ age + pct_active_vietnamese, data = cross_lang))
## 
## Call:
## lm(formula = overlap_amount ~ age + pct_active_vietnamese, data = cross_lang)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4624 -0.5035  0.1288  0.6184  1.5552 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            1.354982   0.366589   3.696 0.000314 ***
## age                   -0.039500   0.012903  -3.061 0.002646 ** 
## pct_active_vietnamese -0.002079   0.002318  -0.897 0.371207    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9032 on 139 degrees of freedom
## Multiple R-squared:  0.06523,    Adjusted R-squared:  0.05178 
## F-statistic:  4.85 on 2 and 139 DF,  p-value: 0.009206
summary(lm(overlap_amount ~ age + pct_active_vietnamese, 
           data = cross_lang |> filter(production_en > 0, production_vi > 0)))
## 
## Call:
## lm(formula = overlap_amount ~ age + pct_active_vietnamese, data = filter(cross_lang, 
##     production_en > 0, production_vi > 0))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.6242 -0.3224  0.1345  0.5428  1.3331 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            1.930303   0.418965   4.607 1.24e-05 ***
## age                   -0.053725   0.014085  -3.814  0.00024 ***
## pct_active_vietnamese -0.001166   0.003188  -0.366  0.71542    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8393 on 97 degrees of freedom
## Multiple R-squared:  0.1311, Adjusted R-squared:  0.1132 
## F-statistic:  7.32 on 2 and 97 DF,  p-value: 0.001095
ggplot(cross_lang, aes(x = age, y = overlap_amount)) +
  geom_point(aes(color = gap_direction), alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_classic() +
  scale_color_manual(values = c(Overlap = "grey70", `Gap: EN > VI` = "#F8766D", `Gap: VI > EN` = "#00BFC4")) +
  labs(x = "Age (months)", y = "EN-VI noun-ability CI overlap", color = NULL)
## `geom_smooth()` using formula = 'y ~ x'

ggplot(cross_lang |> filter(production_en > 0, production_vi > 0),
       aes(x = age, y = overlap_amount)) +
  geom_point(aes(color = gap_direction), alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_classic() +
  scale_color_manual(values = c(Overlap = "grey70", `Gap: EN > VI` = "#F8766D", `Gap: VI > EN` = "#00BFC4")) +
  labs(x = "Age (months)", y = "EN-VI noun-ability CI overlap (SD units)", color = NULL,
       title = "Excluding children who produced 0 words in either language")
## `geom_smooth()` using formula = 'y ~ x'

Noun ability in Vietnamese tracks with % active in Vietnamese (makes sense!) But, kids who have an even split seem to have both overlap 95% CI and gaps. Some kids whose parents reported high Viet use also show EN > VI gap.

ggplot(cross_lang, aes(x = pct_active_vietnamese, y = ability_noun_vi)) +
  geom_point(aes(color = gap_direction), alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE) +
  theme_classic() +
  scale_color_manual(values = c(Overlap = "grey70", `Gap: EN > VI` = "#F8766D", `Gap: VI > EN` = "#00BFC4")) +
  labs(x = "% active Vietnamese exposure", y = "Ability Noun VI", color = NULL)
## `geom_smooth()` using formula = 'y ~ x'

ggplot(cross_lang |> filter(production_en > 0, production_vi > 0),
       aes(x = pct_active_vietnamese, y = ability_noun_vi)) +
  geom_point(aes(color = gap_direction), alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE) +
  theme_classic() +
  scale_color_manual(values = c(Overlap = "grey70", `Gap: EN > VI` = "#F8766D", `Gap: VI > EN` = "#00BFC4")) +
  labs(x = "% active Vietnamese exposure", y = "Ability Noun VI", color = NULL,
       title = "Excluding children who produced 0 words in either language")
## `geom_smooth()` using formula = 'y ~ x'

ggplot(cross_lang,
       aes(x = pct_active_vietnamese, y = overlap_amount)) +
  geom_point(aes(color = gap_direction), alpha = 0.6) +
  geom_smooth(se = TRUE) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_classic() +
  scale_color_manual(values = c(Overlap = "grey70", `Gap: EN > VI` = "#F8766D", `Gap: VI > EN` = "#00BFC4")) +
  labs(x = "% active Vietnamese exposure", y = "EN-VI noun-ability CI overlap", color = NULL)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

ggplot(cross_lang |> filter(production_en > 0, production_vi > 0),
       aes(x = pct_active_vietnamese, y = overlap_amount)) +
  geom_point(aes(color = gap_direction), alpha = 0.6) +
  geom_smooth(se = TRUE) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_classic() +
  scale_color_manual(values = c(Overlap = "grey70", `Gap: EN > VI` = "#F8766D", `Gap: VI > EN` = "#00BFC4")) +
  labs(x = "% active Vietnamese exposure", y = "EN-VI noun-ability CI overlap", color = NULL,
       title = "Excluding children who produced 0 words in either language")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Conclusions / Questions

The majority of kids (67.6-81%, depending on whether we count the kids who do not produce words in one language) show “noun-dominance” that are consistent between English and Vietnamese. This would be consistent with the hypothesis that noun/verb dominance are underpinned by conceptual differences.

Another possibility is that kids tend to get exposed to translation equivalents. We would need to establish that their language input in the languages are actually different, but I think this is unlikely…

Questions:

  1. How do we test whether language exposure determine overlap or not?

  2. For kids who are more “balanced” in language exposure, what determines whether they have better ability in English vs. Vietnamese?