Reproducibility note. Place all TIMSS 2023 Grade 8 student background files (bsg*m8.sav) in the same folder as this .Rmd file and knit. Cache is enabled (cache=TRUE); clear the cache folder before re-running from scratch.

1. Packages

required <- c(
  "haven", "dplyr", "tidyr", "purrr",
  "qgraph", "bootnet", "igraph",
  "nnet", "pROC", "caret",
  "ggplot2", "gridExtra", "ggrepel",
  "factoextra", "dendextend",
  "psych", "knitr", "kableExtra"
)
for (pkg in required) {
  if (!requireNamespace(pkg, quietly = TRUE))
    install.packages(pkg, repos = "https://cloud.r-project.org")
  library(pkg, character.only = TRUE)
}
## package 'ggrepel' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\KNUE\AppData\Local\Temp\RtmpEXjzFg\downloaded_packages
## package 'fracdiff' successfully unpacked and MD5 sums checked
## package 'lmtest' successfully unpacked and MD5 sums checked
## package 'urca' successfully unpacked and MD5 sums checked
## package 'zoo' successfully unpacked and MD5 sums checked
## package 'RcppArmadillo' successfully unpacked and MD5 sums checked
## package 'Deriv' successfully unpacked and MD5 sums checked
## package 'forecast' successfully unpacked and MD5 sums checked
## package 'modelr' successfully unpacked and MD5 sums checked
## package 'microbenchmark' successfully unpacked and MD5 sums checked
## package 'doBy' successfully unpacked and MD5 sums checked
## package 'SparseM' successfully unpacked and MD5 sums checked
## package 'MatrixModels' successfully unpacked and MD5 sums checked
## package 'lazyeval' successfully unpacked and MD5 sums checked
## package 'later' successfully unpacked and MD5 sums checked
## package 'otel' successfully unpacked and MD5 sums checked
## package 'carData' successfully unpacked and MD5 sums checked
## package 'pbkrtest' successfully unpacked and MD5 sums checked
## package 'quantreg' successfully unpacked and MD5 sums checked
## package 'crosstalk' successfully unpacked and MD5 sums checked
## package 'promises' successfully unpacked and MD5 sums checked
## package 'estimability' successfully unpacked and MD5 sums checked
## package 'car' successfully unpacked and MD5 sums checked
## package 'DT' successfully unpacked and MD5 sums checked
## package 'emmeans' successfully unpacked and MD5 sums checked
## package 'flashClust' successfully unpacked and MD5 sums checked
## package 'leaps' successfully unpacked and MD5 sums checked
## package 'multcompView' successfully unpacked and MD5 sums checked
## package 'scatterplot3d' successfully unpacked and MD5 sums checked
## package 'ggsci' successfully unpacked and MD5 sums checked
## package 'cowplot' successfully unpacked and MD5 sums checked
## package 'ggsignif' successfully unpacked and MD5 sums checked
## package 'rstatix' successfully unpacked and MD5 sums checked
## package 'FactoMineR' successfully unpacked and MD5 sums checked
## package 'ggpubr' successfully unpacked and MD5 sums checked
## package 'factoextra' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\KNUE\AppData\Local\Temp\RtmpEXjzFg\downloaded_packages

2. Item Structure

# ── Scale items ───────────────────────────────────────────────
slm_items <- paste0("BSBM19", LETTERS[1:9])   # 9 items
scm_items <- paste0("BSBM22", LETTERS[1:8])   # 8 items
svm_items <- paste0("BSBM23", LETTERS[1:9])   # 9 items
all_items <- c(slm_items, scm_items, svm_items)  # 26 total
label_cols <- c("BSDGSLM", "BSDGSCM", "BSDGSVM")

# Reverse-scored items
all_reverse <- c("BSBM19B","BSBM19C",
                 "BSBM22B","BSBM22C","BSBM22G","BSBM22H")

construct_vec <- c(rep("SLM",9), rep("SCM",8), rep("SVM",9))
names(construct_vec) <- all_items

item_content <- c(
  "Enjoy learning mathematics",
  "Wish I did not have to study mathematics (R)",
  "Mathematics is boring (R)",
  "Learn many interesting things in mathematics",
  "Like mathematics",
  "Like schoolwork involving numbers",
  "Like to solve mathematics problems",
  "Look forward to mathematics class",
  "Mathematics is one of my favorite subjects",
  "I usually do well in mathematics",
  "Mathematics is harder for me than for many classmates (R)",
  "Mathematics is not one of my strengths (R)",
  "Mathematics is easy for me",
  "I am good at working out difficult mathematics problems",
  "I am good at explaining mathematics to others",
  "Mathematics is harder for me than any other subject (R)",
  "Mathematics makes me confused (R)",
  "Math will help me in my daily life",
  "I need mathematics to learn other school subjects",
  "I need to do well to get into the university I want",
  "I need to do well in mathematics to get the job I want",
  "I would like a job that involves using mathematics",
  "It is important to learn mathematics to get ahead in the world",
  "Learning mathematics will give me more job opportunities",
  "My parents think it is important that I do well in mathematics",
  "It is important to do well in mathematics"
)
names(item_content) <- all_items

item_structure <- list(
  SLM = list(items = slm_items, label_var = "BSDGSLM",
             color = "#3498DB", label = "SLM: Students Like Learning Mathematics"),
  SCM = list(items = scm_items, label_var = "BSDGSCM",
             color = "#27AE60", label = "SCM: Students Confident in Mathematics"),
  SVM = list(items = svm_items, label_var = "BSDGSVM",
             color = "#9B59B6", label = "SVM: Students Value Mathematics")
)
cols_map <- setNames(sapply(item_structure, `[[`, "color"),
                     names(item_structure))
n_items <- length(all_items)  # 26

# Geographic region mapping
region_map <- c(
  KOR="East Asia",  JPN="East Asia",  TWN="East Asia",
  HKG="East Asia",  SGP="East Asia",
  ENG="Europe",     FIN="Europe",     FRA="Europe",
  SWE="Europe",     NOR="Europe",     IRL="Europe",
  ITA="Europe",     CZE="Europe",     HUN="Europe",
  LTU="Europe",     MLT="Europe",     CYP="Europe",
  NZL="Europe",     GEO="Europe",     PRT="Europe",
  ROM="Europe",     AUT="Europe",     AUS="Europe",
  ASH="Europe",
  JOR="Middle East/Africa",  SAU="Middle East/Africa",
  QAT="Middle East/Africa",  BHR="Middle East/Africa",
  OMN="Middle East/Africa",  MAR="Middle East/Africa",
  IRN="Middle East/Africa",  PSE="Middle East/Africa",
  ZAF="Middle East/Africa",  CIV="Middle East/Africa",
  KWT="Middle East/Africa",  ARE="Middle East/Africa",
  USA="Americas",   BRA="Americas",   CHL="Americas",
  ISR="Other",      KAZ="Other",      AZE="Other",
  UZB="Other",      TUR="Other",      MYS="Other"
)
region_colors <- c(
  "East Asia"          = "#E74C3C",
  "Europe"             = "#3498DB",
  "Middle East/Africa" = "#F39C12",
  "Americas"           = "#27AE60",
  "Other"              = "#9B59B6"
)
Table 1. TIMSS 2023 Grade 8 Mathematics Affective Items
Item Scale Content Reverse
BSBM19A BSBM19A SLM Enjoy learning mathematics FALSE
BSBM19B BSBM19B SLM Wish I did not have to study mathematics (R) TRUE
BSBM19C BSBM19C SLM Mathematics is boring (R) TRUE
BSBM19D BSBM19D SLM Learn many interesting things in mathematics FALSE
BSBM19E BSBM19E SLM Like mathematics FALSE
BSBM19F BSBM19F SLM Like schoolwork involving numbers FALSE
BSBM19G BSBM19G SLM Like to solve mathematics problems FALSE
BSBM19H BSBM19H SLM Look forward to mathematics class FALSE
BSBM19I BSBM19I SLM Mathematics is one of my favorite subjects FALSE
BSBM22A BSBM22A SCM I usually do well in mathematics FALSE
BSBM22B BSBM22B SCM Mathematics is harder for me than for many classmates (R) TRUE
BSBM22C BSBM22C SCM Mathematics is not one of my strengths (R) TRUE
BSBM22D BSBM22D SCM Mathematics is easy for me FALSE
BSBM22E BSBM22E SCM I am good at working out difficult mathematics problems FALSE
BSBM22F BSBM22F SCM I am good at explaining mathematics to others FALSE
BSBM22G BSBM22G SCM Mathematics is harder for me than any other subject (R) TRUE
BSBM22H BSBM22H SCM Mathematics makes me confused (R) TRUE
BSBM23A BSBM23A SVM Math will help me in my daily life FALSE
BSBM23B BSBM23B SVM I need mathematics to learn other school subjects FALSE
BSBM23C BSBM23C SVM I need to do well to get into the university I want FALSE
BSBM23D BSBM23D SVM I need to do well in mathematics to get the job I want FALSE
BSBM23E BSBM23E SVM I would like a job that involves using mathematics FALSE
BSBM23F BSBM23F SVM It is important to learn mathematics to get ahead in the world FALSE
BSBM23G BSBM23G SVM Learning mathematics will give me more job opportunities FALSE
BSBM23H BSBM23H SVM My parents think it is important that I do well in mathematics FALSE
BSBM23I BSBM23I SVM It is important to do well in mathematics FALSE

3. Data Loading

# SAV files must be in the same directory as this Rmd
DATA_DIR  <- "."
sav_files <- list.files(DATA_DIR, pattern = "^bsg.*m8\\.sav$",
                        full.names = TRUE)
cat(sprintf("Found %d country files\n", length(sav_files)))
## Found 45 country files
read_one <- function(fpath) {
  tryCatch({
    df <- haven::read_spss(fpath) %>%
      select(all_of(c(all_items, label_cols, "IDCNTRY", "CTY"))) %>%
      mutate(across(all_of(all_items),
                    ~ ifelse(. == 9, NA_real_, .))) %>%
      drop_na()
    if (nrow(df) < 100) return(NULL)
    df
  }, error = function(e) {
    message("Skipped: ", basename(fpath), " — ", conditionMessage(e))
    NULL
  })
}

data_raw <- sav_files %>% map(read_one) %>% compact() %>% bind_rows()

# Reverse scoring (4-point scale: 5 - x)
data_raw <- data_raw %>%
  mutate(across(all_of(all_reverse), ~ 5 - .))

# TIMSS official 3-level indices as ordered factor
# 1 = High, 2 = Intermediate, 3 = Low
for (lc in label_cols) {
  data_raw[[lc]] <- factor(data_raw[[lc]], levels = 1:3,
                           labels = c("High","Intermediate","Low"))
}

data <- data_raw
cat(sprintf("N = %s students from %d countries\n",
            format(nrow(data), big.mark = ","),
            n_distinct(data$IDCNTRY)))
## N = 232,970 students from 45 countries

4. Descriptive Statistics

cat(sprintf("Total students : %s\n", format(nrow(data), big.mark=",")))
## Total students : 232,970
cat(sprintf("Total countries: %d\n",  n_distinct(data$IDCNTRY)))
## Total countries: 45
cat(sprintf("Total items    : %d (SLM=%d, SCM=%d, SVM=%d)\n",
            n_items, 9L, 8L, 9L))
## Total items    : 26 (SLM=9, SCM=8, SVM=9)
country_n <- data %>%
  count(CTY, name = "N") %>%
  mutate(Region = unname(ifelse(
    as.character(CTY) %in% names(region_map),
    region_map[as.character(CTY)], "Other"))) %>%
  arrange(desc(N))

ggplot(country_n, aes(reorder(CTY, N), N, fill = Region)) +
  geom_col(alpha = 0.85) +
  geom_text(aes(label = format(N, big.mark=",")),
            hjust = -0.15, size = 2.8) +
  scale_fill_manual(values = region_colors) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.22))) +
  coord_flip() +
  labs(title = "TIMSS 2023 Grade 8 — Sample Size by Country",
       x = NULL, y = "N students") +
  theme_classic(base_size = 10) +
  theme(legend.position = "bottom",
        axis.text.y = element_text(size = 8))
Figure 1. Sample size by country, coloured by region.

Figure 1. Sample size by country, coloured by region. 

data %>%
  select(all_of(label_cols)) %>%
  pivot_longer(everything(), names_to = "Scale", values_to = "Level") %>%
  mutate(Scale = recode(Scale,
    BSDGSLM="SLM", BSDGSCM="SCM", BSDGSVM="SVM")) %>%
  count(Scale, Level) %>%
  group_by(Scale) %>%
  mutate(Pct = round(n / sum(n) * 100, 1)) %>%
  ggplot(aes(Scale, Pct, fill = Level)) +
  geom_col(position = "stack", width = 0.55, alpha = 0.88) +
  geom_text(aes(label = paste0(Pct, "%")),
            position = position_stack(vjust = 0.5),
            fontface = "bold", size = 4.5) +
  scale_fill_manual(values = c(High     = "#2ECC71",
                                Intermediate = "#F39C12",
                                Low      = "#E74C3C")) +
  labs(title = "TIMSS 3-Level Distribution by Scale (Overall)",
       x = NULL, y = "Percentage (%)") +
  theme_classic(base_size = 12) +
  theme(legend.position = "bottom")
Figure 2. TIMSS 3-level distribution by scale (overall).

Figure 2. TIMSS 3-level distribution by scale (overall). 

psych::describe(data[, all_items])[,
  c("n","mean","sd","min","max","skew","kurtosis")] %>%
  round(3) %>%
  kable(caption = "Table 2. Item Descriptive Statistics (after reverse scoring)") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE) %>%
  pack_rows("SLM", 1,  9) %>%
  pack_rows("SCM", 10, 17) %>%
  pack_rows("SVM", 18, 26)
Table 2. Item Descriptive Statistics (after reverse scoring)
n mean sd min max skew kurtosis
SLM
BSBM19A 232970 2.184 1.047 1 4 0.450 -0.991
BSBM19B 232970 2.468 1.135 1 4 0.028 -1.397
BSBM19C 232970 2.497 1.083 1 4 -0.017 -1.278
BSBM19D 232970 2.118 1.011 1 4 0.494 -0.874
BSBM19E 232970 2.289 1.100 1 4 0.308 -1.227
BSBM19F 232970 2.500 1.071 1 4 -0.012 -1.247
BSBM19G 232970 2.420 1.092 1 4 0.129 -1.282
BSBM19H 232970 2.628 1.093 1 4 -0.144 -1.288
BSBM19I 232970 2.586 1.183 1 4 -0.097 -1.495
SCM
BSBM22A 232970 2.056 0.972 1 4 0.587 -0.656
BSBM22B 232970 2.465 1.070 1 4 0.022 -1.246
BSBM22C 232970 2.606 1.110 1 4 -0.148 -1.320
BSBM22D 232970 2.426 1.033 1 4 0.120 -1.135
BSBM22E 232970 2.517 1.032 1 4 0.020 -1.151
BSBM22F 232970 2.527 1.050 1 4 0.001 -1.199
BSBM22G 232970 2.472 1.115 1 4 0.029 -1.353
BSBM22H 232970 2.603 1.072 1 4 -0.151 -1.227
SVM
BSBM23A 232970 1.942 0.970 1 4 0.750 -0.469
BSBM23B 232970 1.997 0.937 1 4 0.649 -0.478
BSBM23C 232970 1.745 0.917 1 4 1.066 0.176
BSBM23D 232970 1.871 0.972 1 4 0.830 -0.415
BSBM23E 232970 2.507 1.100 1 4 -0.006 -1.318
BSBM23F 232970 1.928 0.958 1 4 0.753 -0.438
BSBM23G 232970 1.777 0.899 1 4 1.012 0.189
BSBM23H 232970 1.633 0.848 1 4 1.268 0.844
BSBM23I 232970 1.658 0.861 1 4 1.247 0.793

5. GGM-EBIC (Overall)

set.seed(42)
network_overall <- bootnet::estimateNetwork(
  data[, all_items],
  default   = "EBICglasso",
  corMethod = "cor",
  tuning    = 0.5,
  missing   = "pairwise"
)
pc_overall <- network_overall$graph
rownames(pc_overall) <- colnames(pc_overall) <- all_items

n_edges <- sum(abs(pc_overall[upper.tri(pc_overall)]) > 0.03)
cat(sprintf("Non-zero edges (|r| > 0.03): %d\n", n_edges))
## Non-zero edges (|r| > 0.03): 102
# Cross-construct edges
cross_edges <- do.call(rbind, lapply(seq_along(all_items), function(i) {
  do.call(rbind, lapply(seq_along(all_items), function(j) {
    if (i < j &&
        abs(pc_overall[i,j]) > 0.03 &&
        construct_vec[all_items[i]] != construct_vec[all_items[j]]) {
      data.frame(
        Item1 = all_items[i], Item2 = all_items[j],
        Scale1 = construct_vec[all_items[i]],
        Scale2 = construct_vec[all_items[j]],
        PartialCorr = round(pc_overall[i,j], 3),
        stringsAsFactors = FALSE
      )
    }
  }))
})) %>% arrange(desc(abs(PartialCorr)))
group_list   <- lapply(item_structure, `[[`, "items")
names(group_list) <- sapply(item_structure, `[[`, "label")

qgraph(pc_overall,
       layout     = "spring",
       groups     = group_list,
       color      = unname(sapply(item_structure, `[[`, "color")),
       labels     = all_items,
       label.cex  = 0.85, vsize = 6,
       posCol     = "#1A5276",
       negCol     = "#922B21",
       edge.width = 1.5,
       title      = sprintf(
         "GGM — TIMSS 2023 Affective Scales (EBICglasso, N=%s, %d countries)",
         format(nrow(data), big.mark=","), n_distinct(data$IDCNTRY)),
       legend     = TRUE, legend.cex = 0.55)
Figure 3. Overall GGM network (EBICglasso). Node size proportional to betweenness centrality. Bold edges = cross-construct connections.

Figure 3. Overall GGM network (EBICglasso). Node size proportional to betweenness centrality. Bold edges = cross-construct connections. 

kable(head(cross_edges, 15),
      caption = "Table 3. Strongest Cross-construct Partial Correlations (top 15)") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE)
Table 3. Strongest Cross-construct Partial Correlations (top 15)
Item1 Item2 Scale1 Scale2 PartialCorr
BSBM19G BSBM19G BSBM22E SLM SCM 0.111
BSBM19D BSBM19D BSBM23A SLM SVM 0.107
BSBM19C2 BSBM19C BSBM22H SLM SCM 0.105
BSBM19A BSBM19A BSBM22A SLM SCM 0.093
BSBM19F BSBM19F BSBM23E SLM SVM 0.092
BSBM19I1 BSBM19I BSBM22D SLM SCM 0.086
BSBM22E BSBM22E BSBM23E SCM SVM 0.080
BSBM19I2 BSBM19I BSBM23E SLM SVM 0.078
BSBM19B1 BSBM19B BSBM22C SLM SCM 0.068
BSBM19B4 BSBM19B BSBM23A SLM SVM 0.060
BSBM19I BSBM19I BSBM22C SLM SCM 0.058
BSBM19B3 BSBM19B BSBM22H SLM SCM 0.055
BSBM19A1 BSBM19A BSBM23A SLM SVM 0.049
BSBM19E1 BSBM19E BSBM22D SLM SCM 0.048
BSBM22D BSBM22D BSBM23E SCM SVM 0.045

6. Bridge Centrality

ig_overall <- igraph::graph_from_adjacency_matrix(
  ifelse(abs(pc_overall) > 0.03, 1, 0),
  mode = "undirected", diag = FALSE)
igraph::V(ig_overall)$name      <- all_items
igraph::V(ig_overall)$construct <- construct_vec

between_overall <- igraph::betweenness(ig_overall, normalized = TRUE)

bridge_df_overall <- data.frame(
  Item        = all_items,
  Scale       = construct_vec,
  Content     = item_content,
  Betweenness = round(between_overall, 4),
  Degree      = igraph::degree(ig_overall),
  stringsAsFactors = FALSE
) %>% arrange(desc(Betweenness))

bridge_order_overall <- bridge_df_overall$Item

kable(bridge_df_overall,
      caption = "Table 4. Bridge Centrality — All 26 Items") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE) %>%
  row_spec(1:5, background = "#FDEBD0")
Table 4. Bridge Centrality — All 26 Items
Item Scale Content Betweenness Degree
BSBM23E BSBM23E SVM I would like a job that involves using mathematics 0.1736 12
BSBM23A BSBM23A SVM Math will help me in my daily life 0.0699 9
BSBM19D BSBM19D SLM Learn many interesting things in mathematics 0.0692 9
BSBM23H BSBM23H SVM My parents think it is important that I do well in mathematics 0.0630 6
BSBM22A BSBM22A SCM I usually do well in mathematics 0.0561 8
BSBM22B BSBM22B SCM Mathematics is harder for me than for many classmates (R) 0.0450 8
BSBM19A BSBM19A SLM Enjoy learning mathematics 0.0407 10
BSBM19E BSBM19E SLM Like mathematics 0.0359 10
BSBM19B BSBM19B SLM Wish I did not have to study mathematics (R) 0.0330 8
BSBM22H BSBM22H SCM Mathematics makes me confused (R) 0.0301 6
BSBM19I BSBM19I SLM Mathematics is one of my favorite subjects 0.0299 9
BSBM19C BSBM19C SLM Mathematics is boring (R) 0.0282 8
BSBM19H BSBM19H SLM Look forward to mathematics class 0.0274 9
BSBM23B BSBM23B SVM I need mathematics to learn other school subjects 0.0248 8
BSBM23I BSBM23I SVM It is important to do well in mathematics 0.0235 8
BSBM19G BSBM19G SLM Like to solve mathematics problems 0.0224 8
BSBM23G BSBM23G SVM Learning mathematics will give me more job opportunities 0.0218 8
BSBM22D BSBM22D SCM Mathematics is easy for me 0.0214 8
BSBM22C BSBM22C SCM Mathematics is not one of my strengths (R) 0.0197 8
BSBM23F BSBM23F SVM It is important to learn mathematics to get ahead in the world 0.0196 8
BSBM22F BSBM22F SCM I am good at explaining mathematics to others 0.0182 6
BSBM22E BSBM22E SCM I am good at working out difficult mathematics problems 0.0151 6
BSBM22G BSBM22G SCM Mathematics is harder for me than any other subject (R) 0.0123 6
BSBM23D BSBM23D SVM I need to do well in mathematics to get the job I want 0.0095 6
BSBM19F BSBM19F SLM Like schoolwork involving numbers 0.0060 7
BSBM23C BSBM23C SVM I need to do well to get into the university I want 0.0036 5 
bridge_df_overall %>%
  mutate(Item = factor(Item, levels = rev(bridge_df_overall$Item))) %>%
  ggplot(aes(Item, Betweenness, fill = Scale)) +
  geom_col(alpha = 0.88, width = 0.75) +
  geom_text(aes(label = round(Betweenness, 3)),
            hjust = -0.15, size = 2.8) +
  scale_fill_manual(values = cols_map) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.22))) +
  coord_flip() +
  labs(title    = "Bridge Centrality — TIMSS 2023 Affective Items",
       subtitle = "Items ordered by betweenness centrality (GGM-EBIC)",
       x = NULL, y = "Betweenness Centrality (normalized)") +
  theme_classic(base_size = 10) +
  theme(legend.position  = "bottom",
        axis.text.y      = element_text(size = 8),
        plot.title       = element_text(face = "bold"))
Figure 4. Betweenness centrality for all 26 items.

Figure 4. Betweenness centrality for all 26 items.

7. Evaluation Functions

# ── Overall evaluation (full sample) ─────────────────────────
eval_multinomial <- function(subset_items, y, n_splits = 5, seed = 42) {
  set.seed(seed)
  X        <- as.matrix(data[, subset_items, drop = FALSE])
  n        <- nrow(X)
  idx      <- sample(seq_len(n))
  test_idx <- idx[seq_len(floor(n * 0.2))]
  dev_idx  <- idx[-seq_len(floor(n * 0.2))]

  folds  <- caret::createFolds(y[dev_idx], k = n_splits, list = TRUE)
  cv_au  <- numeric(n_splits)

  for (f in seq_along(folds)) {
    val_i <- dev_idx[folds[[f]]]
    tr_i  <- dev_idx[-folds[[f]]]
    fit   <- nnet::multinom(
      y ~ .,
      data  = data.frame(y = y[tr_i], X[tr_i, , drop = FALSE]),
      trace = FALSE)
    prob  <- predict(fit,
                     newdata = data.frame(X[val_i, , drop = FALSE]),
                     type    = "probs")
    cv_au[f] <- tryCatch({
      mean(sapply(levels(y), function(lv) {
        bin_y <- as.integer(y[val_i] == lv)
        if (length(unique(bin_y)) < 2) return(NA_real_)
        as.numeric(pROC::auc(bin_y, prob[, lv], quiet = TRUE))
      }), na.rm = TRUE)
    }, error = function(e) NA_real_)
  }

  fit_f  <- nnet::multinom(
    y ~ .,
    data  = data.frame(y = y[dev_idx], X[dev_idx, , drop = FALSE]),
    trace = FALSE)
  prob_t <- predict(fit_f,
                    newdata = data.frame(X[test_idx, , drop = FALSE]),
                    type    = "probs")
  pred_t <- predict(fit_f,
                    newdata = data.frame(X[test_idx, , drop = FALSE]))

  test_au <- tryCatch({
    mean(sapply(levels(y), function(lv) {
      bin_y <- as.integer(y[test_idx] == lv)
      if (length(unique(bin_y)) < 2) return(NA_real_)
      as.numeric(pROC::auc(bin_y, prob_t[, lv], quiet = TRUE))
    }), na.rm = TRUE)
  }, error = function(e) NA_real_)

  cm <- caret::confusionMatrix(pred_t, y[test_idx])
  list(cv_auroc   = mean(cv_au, na.rm = TRUE),
       cv_sd      = sd(cv_au,   na.rm = TRUE),
       test_auroc = test_au,
       kappa      = as.numeric(cm$overall["Kappa"]),
       acc        = as.numeric(cm$overall["Accuracy"]))
}

# ── Country-level evaluation ──────────────────────────────────
eval_on_df <- function(df_sub, subset_items, label_var,
                       n_splits = 3, seed = 42) {
  y <- df_sub[[label_var]]
  if (!is.factor(y))
    y <- factor(y, levels = c("High","Intermediate","Low"))
  if (length(levels(droplevels(y))) < 3) return(NULL)

  X <- as.matrix(df_sub[, subset_items, drop = FALSE])
  n <- nrow(X)
  if (n < 100) return(NULL)

  set.seed(seed)
  idx      <- sample(seq_len(n))
  test_idx <- idx[seq_len(floor(n * 0.2))]
  dev_idx  <- idx[-seq_len(floor(n * 0.2))]
  if (length(unique(y[dev_idx])) < 3) return(NULL)

  folds <- caret::createFolds(y[dev_idx], k = n_splits, list = TRUE)
  cv_au <- numeric(n_splits)

  for (f in seq_along(folds)) {
    val_i <- dev_idx[folds[[f]]]; tr_i <- dev_idx[-folds[[f]]]
    if (length(unique(y[tr_i])) < 3) next
    fit <- tryCatch(
      nnet::multinom(
        y ~ .,
        data  = data.frame(y = y[tr_i], X[tr_i, , drop=FALSE]),
        trace = FALSE),
      error = function(e) NULL)
    if (is.null(fit)) next
    prob <- predict(fit,
                    newdata = data.frame(X[val_i, , drop=FALSE]),
                    type    = "probs")
    cv_au[f] <- tryCatch(
      mean(sapply(levels(y), function(lv) {
        b <- as.integer(y[val_i] == lv)
        if (length(unique(b)) < 2) NA_real_
        else as.numeric(pROC::auc(b, prob[, lv], quiet = TRUE))
      }), na.rm = TRUE),
      error = function(e) NA_real_)
  }

  fit_f  <- tryCatch(
    nnet::multinom(
      y ~ .,
      data  = data.frame(y = y[dev_idx], X[dev_idx, , drop=FALSE]),
      trace = FALSE),
    error = function(e) NULL)
  if (is.null(fit_f)) return(NULL)

  prob_t <- predict(fit_f,
                    newdata = data.frame(X[test_idx, , drop=FALSE]),
                    type    = "probs")
  pred_t <- predict(fit_f,
                    newdata = data.frame(X[test_idx, , drop=FALSE]))

  test_au <- tryCatch(
    mean(sapply(levels(y), function(lv) {
      b <- as.integer(y[test_idx] == lv)
      if (length(unique(b)) < 2) NA_real_
      else as.numeric(pROC::auc(b, prob_t[, lv], quiet = TRUE))
    }), na.rm = TRUE),
    error = function(e) NA_real_)

  cm <- tryCatch(
    caret::confusionMatrix(pred_t, y[test_idx]),
    error = function(e) NULL)

  list(cv_auroc   = mean(cv_au, na.rm = TRUE),
       test_auroc = test_au,
       kappa      = if (!is.null(cm)) as.numeric(cm$overall["Kappa"]) else NA_real_,
       acc        = if (!is.null(cm)) as.numeric(cm$overall["Accuracy"]) else NA_real_)
}

8. Single-item Benchmarks

single_results <- lapply(names(item_structure), function(sname) {
  cd <- item_structure[[sname]]
  y  <- data[[cd$label_var]]
  do.call(rbind, lapply(cd$items, function(item) {
    r <- eval_multinomial(item, y)
    data.frame(Scale = sname, Item = item,
               AUROC = r$test_auroc, Kappa = r$kappa,
               CV = r$cv_auroc, CV_SD = r$cv_sd,
               stringsAsFactors = FALSE)
  }))
}) %>% bind_rows() %>% arrange(Scale, desc(AUROC))

best_single <- single_results %>%
  group_by(Scale) %>% slice_max(AUROC, n = 1) %>% ungroup()

kable(single_results %>%
        mutate(across(where(is.numeric), ~ round(., 3))),
      caption = "Table 5. Single Own-scale Item Performance") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE) %>%
  pack_rows("SCM", 1,  8) %>%
  pack_rows("SLM", 9, 17) %>%
  pack_rows("SVM",18, 26)
Table 5. Single Own-scale Item Performance
Scale Item AUROC Kappa CV CV_SD
SCM
SCM BSBM22D 0.839 0.536 0.840 0.002
SCM BSBM22C 0.830 0.530 0.830 0.001
SCM BSBM22E 0.810 0.475 0.808 0.002
SCM BSBM22G 0.803 0.313 0.802 0.002
SCM BSBM22B 0.798 0.316 0.797 0.001
SCM BSBM22A 0.786 0.294 0.785 0.002
SCM BSBM22H 0.782 0.414 0.783 0.002
SCM BSBM22F 0.750 0.222 0.750 0.004
SLM
SLM BSBM19E 0.908 0.714 0.908 0.000
SLM BSBM19I 0.904 0.679 0.905 0.001
SLM BSBM19G 0.890 0.686 0.890 0.002
SLM BSBM19A 0.878 0.630 0.879 0.001
SLM BSBM19H 0.877 0.610 0.879 0.001
SLM BSBM19F 0.866 0.612 0.867 0.002
SLM BSBM19D 0.828 0.522 0.828 0.002
SLM BSBM19C 0.723 0.384 0.722 0.001
SLM BSBM19B 0.671 0.316 0.673 0.002
SVM
SVM BSBM23F 0.853 0.610 0.856 0.001
SVM BSBM23D 0.843 0.597 0.843 0.001
SVM BSBM23G 0.825 0.570 0.825 0.001
SVM BSBM23A 0.817 0.537 0.819 0.002
SVM BSBM23B 0.812 0.521 0.810 0.003
SVM BSBM23C 0.811 0.543 0.810 0.001
SVM BSBM23E 0.810 0.420 0.812 0.001
SVM BSBM23I 0.793 0.499 0.792 0.001
SVM BSBM23H 0.709 0.361 0.707 0.001

9. Cross-construct Sequential Addition

plateau_all <- list()

for (sname in names(item_structure)) {
  cd      <- item_structure[[sname]]
  y       <- data[[cd$label_var]]
  base_au <- best_single$AUROC[best_single$Scale == sname]
  seq_res <- list()

  cat(sprintf("\n[%s] Baseline AUROC = %.3f\n", sname, base_au))

  for (k in seq_along(bridge_order_overall)) {
    subset  <- bridge_order_overall[seq_len(k)]
    res     <- eval_multinomial(subset, y)
    n_own   <- sum(subset %in% cd$items)
    n_cross <- sum(!subset %in% cd$items)
    seq_res[[k]] <- c(
      list(k = k, items = list(subset),
           n_own = n_own, n_cross = n_cross), res)

    cat(sprintf("  k=%2d [own=%d cross=%d]: AUROC=%.3f Kappa=%.3f\n",
                k, n_own, n_cross, res$test_auroc, res$kappa))

    if (!is.na(res$test_auroc) &&
        res$test_auroc >= base_au * 0.99) {
      cat(sprintf("  --> Plateau at k=%d\n", k))
      break
    }
  }
  plateau_all[[sname]] <- seq_res
}
## 
## [SLM] Baseline AUROC = 0.908
##   k= 1 [own=0 cross=1]: AUROC=0.742 Kappa=0.372
##   k= 2 [own=0 cross=2]: AUROC=0.778 Kappa=0.398
##   k= 3 [own=1 cross=2]: AUROC=0.873 Kappa=0.552
##   k= 4 [own=1 cross=3]: AUROC=0.874 Kappa=0.553
##   k= 5 [own=1 cross=4]: AUROC=0.893 Kappa=0.601
##   k= 6 [own=1 cross=5]: AUROC=0.894 Kappa=0.602
##   k= 7 [own=2 cross=5]: AUROC=0.939 Kappa=0.720
##   --> Plateau at k=7
## 
## [SCM] Baseline AUROC = 0.839
##   k= 1 [own=0 cross=1]: AUROC=0.650 Kappa=0.097
##   k= 2 [own=0 cross=2]: AUROC=0.658 Kappa=0.077
##   k= 3 [own=0 cross=3]: AUROC=0.685 Kappa=0.138
##   k= 4 [own=0 cross=4]: AUROC=0.688 Kappa=0.123
##   k= 5 [own=1 cross=4]: AUROC=0.808 Kappa=0.353
##   k= 6 [own=2 cross=4]: AUROC=0.914 Kappa=0.641
##   --> Plateau at k=6
## 
## [SVM] Baseline AUROC = 0.853
##   k= 1 [own=1 cross=0]: AUROC=0.810 Kappa=0.420
##   k= 2 [own=2 cross=0]: AUROC=0.889 Kappa=0.583
##   --> Plateau at k=2
seq_plots <- lapply(names(item_structure), function(sname) {
  seq_res <- plateau_all[[sname]]
  base_au <- best_single$AUROC[best_single$Scale == sname]
  col     <- item_structure[[sname]]$color
  k_plat  <- seq_res[[length(seq_res)]]$k

  df_seq <- data.frame(
    k         = sapply(seq_res, `[[`, "k"),
    auroc     = sapply(seq_res, `[[`, "test_auroc"),
    kappa     = sapply(seq_res, `[[`, "kappa"),
    n_own     = sapply(seq_res, `[[`, "n_own"),
    n_cross   = sapply(seq_res, `[[`, "n_cross"),
    last_item = sapply(seq_res, function(r) {
      items <- unlist(r$items); items[length(items)]})
  ) %>%
    mutate(item_scale = construct_vec[last_item],
           is_cross   = item_scale != sname)

  ggplot(df_seq) +
    geom_hline(yintercept = base_au,
               linetype = "dotted", color = "gray50", linewidth = 1.2) +
    geom_hline(yintercept = base_au * 0.99,
               linetype = "dotted", color = "#E74C3C", linewidth = 0.9) +
    geom_vline(xintercept = k_plat,
               linetype = "dashed", color = "#E74C3C", linewidth = 1) +
    geom_line(aes(k, auroc), color = col, linewidth = 2.2) +
    geom_point(aes(k, auroc, shape = is_cross, color = item_scale),
               size = 4) +
    geom_line(aes(k, kappa), color = col,
              linewidth = 1.5, linetype = "dashed", alpha = 0.7) +
    annotate("text", x = k_plat + 0.3,
             y = df_seq$auroc[nrow(df_seq)] - 0.06,
             label = sprintf("Plateau\nk=%d\n(own=%d, cross=%d)",
                             k_plat,
                             df_seq$n_own[nrow(df_seq)],
                             df_seq$n_cross[nrow(df_seq)]),
             color = "#E74C3C", fontface = "bold",
             size = 3.2, hjust = 0) +
    scale_color_manual(values = cols_map, name = "Item scale") +
    scale_shape_manual(values = c(`FALSE` = 16, `TRUE` = 17),
                       labels = c("Own-scale","Cross-scale"),
                       name   = "Item type") +
    scale_x_continuous(breaks = seq(1, n_items, 2)) +
    labs(title    = sprintf("%s — Cross-construct Sequential Addition", sname),
         subtitle = sprintf("Baseline = %.3f  |  Plateau k = %d",
                            base_au, k_plat),
         x = "Items added (bridge-first, all 26 pooled)",
         y = "AUROC (solid) / Kappa (dashed)") +
    ylim(0.4, 1.05) +
    theme_classic(base_size = 11) +
    theme(plot.title    = element_text(face = "bold", size = 11),
          plot.subtitle = element_text(size = 9),
          legend.position = "right")
})

gridExtra::grid.arrange(
  grobs = seq_plots, ncol = 3,
  top   = grid::textGrob(
    "Design B: Cross-construct Bridge-first Sequential Addition\n(circles = own-scale; triangles = cross-scale)",
    gp = grid::gpar(fontsize = 13, fontface = "bold")))
Figure 5. Cross-construct sequential addition curves. Circles = own-scale items; triangles = cross-scale items.

Figure 5. Cross-construct sequential addition curves. Circles = own-scale items; triangles = cross-scale items. 

summary_overall <- do.call(rbind, lapply(names(item_structure), function(sname) {
  plat  <- plateau_all[[sname]]
  last  <- plat[[length(plat)]]
  items <- unlist(last$items)
  base  <- best_single$AUROC[best_single$Scale == sname]
  data.frame(
    Scale          = sname,
    Full_items     = length(item_structure[[sname]]$items),
    Plateau_k      = last$k,
    Own            = last$n_own,
    Cross          = last$n_cross,
    Pct_cross      = paste0(round(last$n_cross / last$k * 100), "%"),
    Baseline_AUROC = round(base,             3),
    Plateau_AUROC  = round(last$test_auroc,  3),
    Delta          = round(last$test_auroc - base, 3),
    Plateau_Kappa  = round(last$kappa,       3),
    Selected       = paste(items, collapse = ", "),
    row.names      = NULL, stringsAsFactors = FALSE)
}))

kable(summary_overall,
      caption = "Table 6. Cross-construct Sequential Selection Results (Overall, N = 232,970)") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = TRUE) %>%
  column_spec(8, bold = TRUE) %>%
  column_spec(9, bold = TRUE, color = "white",
              background = ifelse(summary_overall$Delta >= 0,
                                  "#27AE60","#E74C3C"))
Table 6. Cross-construct Sequential Selection Results (Overall, N = 232,970)
Scale Full_items Plateau_k Own Cross Pct_cross Baseline_AUROC Plateau_AUROC Delta Plateau_Kappa Selected
SLM 9 7 2 5 71% 0.908 0.939 0.031 0.720 BSBM23E, BSBM23A, BSBM19D, BSBM23H, BSBM22A, BSBM22B, BSBM19A
SCM 8 6 2 4 67% 0.839 0.914 0.075 0.641 BSBM23E, BSBM23A, BSBM19D, BSBM23H, BSBM22A, BSBM22B
SVM 9 2 2 0 0% 0.853 0.889 0.036 0.583 BSBM23E, BSBM23A

10. Country-level Pre-screening

country_diagnostics <- data %>%
  group_by(CTY) %>%
  summarise(
    N    = n(),
    Region = unname(ifelse(
      as.character(first(CTY)) %in% names(region_map),
      region_map[as.character(first(CTY))], "Other")),
    SLM_H = round(mean(BSDGSLM == "High", na.rm=TRUE), 2),
    SLM_L = round(mean(BSDGSLM == "Low",  na.rm=TRUE), 2),
    SCM_H = round(mean(BSDGSCM == "High", na.rm=TRUE), 2),
    SCM_L = round(mean(BSDGSCM == "Low",  na.rm=TRUE), 2),
    SVM_H = round(mean(BSDGSVM == "High", na.rm=TRUE), 2),
    SVM_L = round(mean(BSDGSVM == "Low",  na.rm=TRUE), 2),
    .groups = "drop"
  ) %>%
  mutate(
    flag_n   = N < 500,
    flag_slm = pmin(SLM_H, SLM_L) < 0.05,
    flag_scm = pmin(SCM_H, SCM_L) < 0.05,
    flag_svm = pmin(SVM_H, SVM_L) < 0.05,
    any_flag = flag_n | flag_slm | flag_scm | flag_svm,
    Risk = case_when(
      flag_n                             ~ "High (small n)",
      flag_slm | flag_scm | flag_svm    ~ "Medium (imbalance)",
      TRUE                               ~ "Low")
  ) %>%
  arrange(Risk, N)

kable(country_diagnostics %>%
        select(CTY, Region, N, Risk, SLM_H, SLM_L,
               SCM_H, SCM_L, SVM_H, SVM_L),
      caption = "Table 7. Country Pre-screening Diagnostics") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE) %>%
  row_spec(which(country_diagnostics$Risk == "High (small n)"),
           background = "#FADBD8") %>%
  row_spec(which(country_diagnostics$Risk == "Medium (imbalance)"),
           background = "#FEF9E7")
Table 7. Country Pre-screening Diagnostics
CTY Region N Risk SLM_H SLM_L SCM_H SCM_L SVM_H SVM_L
ROM Europe 2706 Low 0.19 0.50 0.12 0.63 0.28 0.35
NZL Europe 3115 Low 0.12 0.55 0.12 0.54 0.29 0.25
CYP Europe 3144 Low 0.17 0.56 0.20 0.51 0.29 0.33
UZB Other 3164 Low 0.56 0.12 0.19 0.44 0.71 0.07
KWT Middle East/Africa 3174 Low 0.23 0.48 0.13 0.54 0.40 0.23
MLT Europe 3207 Low 0.16 0.54 0.16 0.54 0.31 0.27
JPN East Asia 3210 Low 0.11 0.58 0.05 0.74 0.14 0.30
ENG Europe 3340 Low 0.11 0.59 0.15 0.48 0.29 0.26
QAT Middle East/Africa 3368 Low 0.26 0.43 0.18 0.50 0.42 0.21
GEO Europe 3382 Low 0.25 0.38 0.18 0.51 0.43 0.19
FRA Europe 3453 Low 0.10 0.57 0.16 0.52 0.24 0.28
CHL Americas 3532 Low 0.17 0.49 0.11 0.61 0.37 0.15
BHR Middle East/Africa 3721 Low 0.28 0.41 0.18 0.50 0.43 0.21
ITA Europe 3821 Low 0.12 0.60 0.17 0.51 0.17 0.41
AUT Europe 3877 Low 0.10 0.62 0.14 0.49 0.12 0.47
HKG East Asia 3885 Low 0.18 0.45 0.09 0.64 0.20 0.31
SAU Middle East/Africa 3895 Low 0.35 0.33 0.17 0.50 0.52 0.14
AZE Other 3909 Low 0.41 0.28 0.16 0.57 0.53 0.15
KOR East Asia 3916 Low 0.12 0.60 0.11 0.61 0.23 0.28
SWE Europe 3948 Low 0.09 0.65 0.15 0.53 0.19 0.32
PSE Middle East/Africa 3995 Low 0.41 0.31 0.18 0.52 0.53 0.17
TUR Other 4150 Low 0.22 0.46 0.14 0.62 0.37 0.25
SGP East Asia 4293 Low 0.23 0.39 0.16 0.54 0.33 0.18
IRL Europe 4309 Low 0.11 0.59 0.13 0.52 0.25 0.29
ISR Other 4327 Low 0.17 0.52 0.21 0.46 0.43 0.20
LTU Europe 4352 Low 0.10 0.58 0.11 0.61 0.22 0.32
JOR Middle East/Africa 4490 Low 0.39 0.30 0.14 0.53 0.56 0.14
PRT Europe 4526 Low 0.13 0.58 0.12 0.66 0.27 0.30
HUN Europe 4604 Low 0.11 0.62 0.19 0.50 0.24 0.29
KAZ Other 4826 Low 0.18 0.39 0.11 0.56 0.25 0.29
TWN East Asia 4910 Low 0.12 0.60 0.10 0.67 0.18 0.41
ASH Europe 4929 Low 0.35 0.33 0.20 0.48 0.52 0.16
NOR Europe 4951 Low 0.09 0.65 0.14 0.54 0.22 0.31
FIN Europe 4964 Low 0.07 0.69 0.14 0.57 0.17 0.37
IRN Middle East/Africa 5149 Low 0.34 0.31 0.14 0.53 0.43 0.22
OMN Middle East/Africa 5152 Low 0.33 0.33 0.15 0.47 0.49 0.15
MAR Middle East/Africa 5600 Low 0.43 0.24 0.15 0.49 0.60 0.12
USA Americas 6572 Low 0.15 0.52 0.16 0.51 0.28 0.24
AUS Europe 6815 Low 0.13 0.54 0.15 0.49 0.31 0.24
CZE Europe 6915 Low 0.08 0.66 0.13 0.53 0.14 0.35
CIV Middle East/Africa 6959 Low 0.52 0.10 0.08 0.50 0.62 0.07
ZAF Middle East/Africa 9095 Low 0.41 0.19 0.11 0.54 0.69 0.06
BRA Americas 15077 Low 0.18 0.50 0.06 0.73 0.43 0.16
ARE Middle East/Africa 24557 Low 0.31 0.36 0.19 0.50 0.47 0.18
MYS Other 7686 Medium (imbalance) 0.20 0.32 0.04 0.70 0.36 0.15 
exclude_ctry <- country_diagnostics %>% filter(any_flag)  %>% pull(CTY)
include_ctry <- country_diagnostics %>% filter(!any_flag) %>% pull(CTY)
cat(sprintf("Included: %d  |  Flagged: %d (%s)\n",
            length(include_ctry), length(exclude_ctry),
            paste(exclude_ctry, collapse = ", ")))
## Included: 44  |  Flagged: 1 (MYS)

11. Country-level Analysis

country_list    <- sort(include_ctry)
country_results <- list()

cat(sprintf("Running loop: %d countries\n", length(country_list)))
## Running loop: 44 countries
for (ctry in country_list) {
  df_c <- data %>% filter(CTY == ctry)
  cat(sprintf("  [%s] n=%d\n", ctry, nrow(df_c)))

  net_c <- tryCatch(
    bootnet::estimateNetwork(
      df_c[, all_items],
      default  = "EBICglasso",
      tuning   = 0.5,
      missing  = "pairwise"),
    error = function(e) { message("GGM failed: ", ctry); NULL })
  if (is.null(net_c)) next

  pc_c <- net_c$graph
  ig_c <- igraph::graph_from_adjacency_matrix(
    ifelse(abs(pc_c) > 0.03, 1, 0),
    mode = "undirected", diag = FALSE)
  bc_c   <- igraph::betweenness(ig_c, normalized = TRUE)
  bord_c <- names(sort(bc_c, decreasing = TRUE))

  ctry_res <- list(pc = pc_c, bridge_order = bord_c,
                   n_edges = sum(abs(pc_c[upper.tri(pc_c)]) > 0.03),
                   n = nrow(df_c))

  for (sname in names(item_structure)) {
    cd     <- item_structure[[sname]]
    y_c    <- df_c[[cd$label_var]]
    single_c <- sapply(cd$items, function(item)
      tryCatch(eval_on_df(df_c, item, cd$label_var,
                          n_splits=3)$test_auroc,
               error = function(e) NA_real_))
    base_c <- max(single_c, na.rm = TRUE)

    seq_c <- list()
    for (k in seq_along(bord_c)) {
      res_c <- tryCatch(
        eval_on_df(df_c, bord_c[seq_len(k)], cd$label_var, n_splits=3),
        error = function(e) NULL)
      if (is.null(res_c)) next
      seq_c[[k]] <- c(list(k=k,
        n_own   = sum(bord_c[seq_len(k)] %in% cd$items),
        n_cross = sum(!bord_c[seq_len(k)] %in% cd$items)), res_c)
      if (!is.na(res_c$test_auroc) &&
          res_c$test_auroc >= base_c * 0.99) break
    }
    ctry_res[[sname]] <- list(baseline = base_c, plateau = seq_c)
  }
  country_results[[ctry]] <- ctry_res
}
##   [ARE] n=24557
##   [ASH] n=4929
##   [AUS] n=6815
##   [AUT] n=3877
##   [AZE] n=3909
##   [BHR] n=3721
##   [BRA] n=15077
##   [CHL] n=3532
##   [CIV] n=6959
##   [CYP] n=3144
##   [CZE] n=6915
##   [ENG] n=3340
##   [FIN] n=4964
##   [FRA] n=3453
##   [GEO] n=3382
##   [HKG] n=3885
##   [HUN] n=4604
##   [IRL] n=4309
##   [IRN] n=5149
##   [ISR] n=4327
##   [ITA] n=3821
##   [JOR] n=4490
##   [JPN] n=3210
##   [KAZ] n=4826
##   [KOR] n=3916
##   [KWT] n=3174
##   [LTU] n=4352
##   [MAR] n=5600
##   [MLT] n=3207
##   [NOR] n=4951
##   [NZL] n=3115
##   [OMN] n=5152
##   [PRT] n=4526
##   [PSE] n=3995
##   [QAT] n=3368
##   [ROM] n=2706
##   [SAU] n=3895
##   [SGP] n=4293
##   [SWE] n=3948
##   [TUR] n=4150
##   [TWN] n=4910
##   [USA] n=6572
##   [UZB] n=3164
##   [ZAF] n=9095
cat("Country loop complete.\n")
## Country loop complete.

12. Cross-national Consistency

12.1 Bridge Item Consistency Index

top_k <- 5

consistency_list <- lapply(names(country_results), function(ctry) {
  cr <- country_results[[ctry]]
  bo <- cr$bridge_order
  overlap <- if (is.null(bo) || length(bo) < top_k) NA_real_
             else mean(bridge_order_overall[seq_len(top_k)] %in%
                         bo[seq_len(top_k)])
  region <- tryCatch(unname(region_map[[ctry]]),
                     error = function(e) "Other")
  if (is.null(region) || length(region)==0 || is.na(region))
    region <- "Other"
  data.frame(Country     = as.character(ctry),
             Region      = as.character(region),
             N           = as.integer(cr$n),
             Consistency = round(overlap, 3),
             stringsAsFactors = FALSE)
})
consistency_df   <- do.call(rbind, consistency_list)
rownames(consistency_df) <- NULL

# Per-scale AUROC and k at plateau
for (sname in names(item_structure)) {
  consistency_df[[paste0(sname,"_AUROC")]] <- sapply(
    names(country_results), function(ctry) {
      cr <- country_results[[ctry]]
      if (is.null(cr[[sname]])) return(NA_real_)
      plat <- cr[[sname]]$plateau
      if (length(plat)==0) return(NA_real_)
      round(plat[[length(plat)]]$test_auroc, 3)
    })
  consistency_df[[paste0(sname,"_k")]] <- sapply(
    names(country_results), function(ctry) {
      cr <- country_results[[ctry]]
      if (is.null(cr[[sname]])) return(NA_integer_)
      plat <- cr[[sname]]$plateau
      if (length(plat)==0) return(NA_integer_)
      plat[[length(plat)]]$k
    })
}

cat(sprintf("Mean consistency index: %.3f (SD = %.3f)\n",
            mean(consistency_df$Consistency, na.rm=TRUE),
            sd(consistency_df$Consistency,   na.rm=TRUE)))
## Mean consistency index: 0.509 (SD = 0.158)
ggplot(consistency_df,
       aes(reorder(Country, Consistency), Consistency, fill = Region)) +
  geom_col(alpha = 0.85) +
  geom_hline(yintercept = mean(consistency_df$Consistency, na.rm=TRUE),
             linetype = "dashed", color = "#2C3E50", linewidth = 1) +
  geom_text(aes(label = round(Consistency, 2)),
            hjust = -0.15, size = 2.8) +
  scale_fill_manual(values = region_colors) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.22)),
                     limits = c(0, 1)) +
  coord_flip() +
  labs(title    = sprintf("Bridge Item Consistency Index (top-%d items)", top_k),
       subtitle = "Proportion of overall top-k bridge items also in country top-k",
       x = NULL, y = "Consistency Index") +
  theme_classic(base_size = 10) +
  theme(legend.position = "bottom",
        axis.text.y     = element_text(size = 8),
        plot.title      = element_text(face = "bold"))
Figure 6. Bridge item consistency index by country and region.

Figure 6. Bridge item consistency index by country and region.

12.2 Country Clustering

ctry_names <- names(country_results)
n_ctry     <- length(ctry_names)
sim_mat    <- matrix(NA, n_ctry, n_ctry,
                     dimnames = list(ctry_names, ctry_names))

for (i in seq_len(n_ctry)) {
  for (j in seq_len(n_ctry)) {
    if (i == j) { sim_mat[i,j] <- 1; next }
    pc_i <- country_results[[ctry_names[i]]]$pc
    pc_j <- country_results[[ctry_names[j]]]$pc
    if (is.null(pc_i) || is.null(pc_j)) next
    common <- intersect(rownames(pc_i), rownames(pc_j))
    if (length(common) < 3) next
    ut <- upper.tri(pc_i[common, common])
    sim_mat[i,j] <- tryCatch(
      cor(pc_i[common,common][ut], pc_j[common,common][ut]),
      error = function(e) NA_real_)
  }
}
sim_mat[is.na(sim_mat)] <- 0
dist_mat <- as.dist(1 - sim_mat)
hc       <- hclust(dist_mat, method = "ward.D2")

# Cluster assignment (k=4)
k_clusters     <- 4
cluster_labels <- cutree(hc, k = k_clusters)
consistency_df$Cluster <- paste0("C",
  cluster_labels[consistency_df$Country])

region_cols <- sapply(ctry_names, function(x)
  ifelse(x %in% names(region_map),
         region_colors[region_map[x]], "gray50"))
dend <- as.dendrogram(hc)
dend <- dendextend::color_labels(
  dend, col = region_cols[order.dendrogram(dend)])

cat(sprintf("Clustering done: k=%d, %d countries.\n",
            k_clusters, n_ctry))
## Clustering done: k=4, 44 countries.
plot(dend,
     main = "Country Clustering by GGM Network Structure\n(Ward linkage, 1 - Pearson r)",
     cex = 0.9, horiz = TRUE)
legend("topleft", legend = names(region_colors),
       fill = region_colors, cex = 0.8, bty = "n")
Figure 7. Hierarchical clustering of countries by GGM network similarity (Ward linkage, 1 - Pearson r).

Figure 7. Hierarchical clustering of countries by GGM network similarity (Ward linkage, 1 - Pearson r).

12.3 AUROC Distribution

consistency_df %>%
  select(Country, Region, SLM_AUROC, SCM_AUROC, SVM_AUROC) %>%
  pivot_longer(ends_with("AUROC"),
               names_to  = "Scale",
               values_to = "AUROC") %>%
  mutate(Scale = sub("_AUROC", "", Scale)) %>%
  ggplot(aes(Scale, AUROC, fill = Scale)) +
  geom_boxplot(alpha = 0.7, outlier.shape = NA) +
  geom_jitter(aes(color = Region), width = 0.15, size = 2.5, alpha = 0.8) +
  geom_hline(yintercept = 0.80, linetype = "dashed",
             color = "gray50", linewidth = 1) +
  geom_hline(yintercept = 0.90, linetype = "dotted",
             color = "#2C3E50", linewidth = 1) +
  scale_fill_manual(values  = cols_map) +
  scale_color_manual(values = region_colors) +
  labs(title    = "Cross-national Distribution of Plateau AUROC",
       subtitle = "Each point = one country | dashed = .80 | dotted = .90",
       x = NULL, y = "AUROC (plateau item set)") +
  theme_classic(base_size = 12) +
  theme(legend.position = "right",
        plot.title      = element_text(face = "bold"))
Figure 8. Distribution of plateau AUROC across countries by scale.

Figure 8. Distribution of plateau AUROC across countries by scale.

13. Representative Countries

rep_countries <- consistency_df %>%
  filter(Country %in% names(country_results)) %>%
  slice_max(Consistency, n = 4, with_ties = FALSE) %>%
  arrange(desc(Consistency), desc(N)) %>%
  pull(Country)

cat("Representative countries (top-4 by consistency index):\n")
## Representative countries (top-4 by consistency index):
print(consistency_df %>%
        filter(Country %in% rep_countries) %>%
        select(Country, Region, N, Consistency,
               SLM_AUROC, SCM_AUROC, SVM_AUROC) %>%
        arrange(desc(Consistency)))
##   Country    Region    N Consistency SLM_AUROC SCM_AUROC SVM_AUROC
## 1     CZE    Europe 6915         0.8     0.931     0.886     0.880
## 2     KOR East Asia 3916         0.8     0.927     0.920     0.902
## 3     NZL    Europe 3115         0.8     0.952     0.899     0.864
## 4     USA  Americas 6572         0.8     0.952     0.900     0.886
rep_plot_list <- list()

for (ctry in rep_countries) {
  cr  <- country_results[[ctry]]
  if (is.null(cr)) next
  reg <- ifelse(ctry %in% names(region_map), region_map[ctry], "Other")

  for (sname in names(item_structure)) {
    sres <- cr[[sname]]
    plat <- if (!is.null(sres)) sres$plateau else NULL
    base <- if (!is.null(sres)) sres$baseline else NULL
    if (is.null(plat) || length(plat) == 0) next
    if (is.null(base)  || is.na(base))       next

    last <- plat[[length(plat)]]
    k_pl <- last$k
    col  <- item_structure[[sname]]$color

    df_p <- data.frame(
      k     = sapply(plat, `[[`, "k"),
      auroc = sapply(plat, function(r)
        ifelse(is.null(r$test_auroc), NA_real_, r$test_auroc)),
      kappa = sapply(plat, function(r)
        ifelse(is.null(r$kappa), NA_real_, r$kappa))
    ) %>% filter(!is.na(auroc))
    if (nrow(df_p) == 0) next

    g <- ggplot(df_p) +
      geom_hline(yintercept = base,
                 linetype = "dotted", color = "gray50",   linewidth = 0.9) +
      geom_hline(yintercept = base * 0.99,
                 linetype = "dotted", color = "#E74C3C",  linewidth = 0.7) +
      geom_vline(xintercept = k_pl,
                 linetype = "dashed", color = "#E74C3C",  linewidth = 0.9) +
      geom_line(aes(k, auroc),  color = col, linewidth = 1.8) +
      geom_point(aes(k, auroc), color = col, size = 2.5) +
      geom_line(aes(k, kappa),  color = col,
                linewidth = 1.2, linetype = "dashed", alpha = 0.65) +
      annotate("text", x = k_pl + 0.4,
               y = min(df_p$auroc, na.rm=TRUE) + 0.04,
               label = sprintf("k=%d", k_pl),
               color = "#E74C3C", fontface = "bold", size = 3.2) +
      labs(title    = sprintf("%s | %s (%s)", ctry, sname, reg),
           subtitle = sprintf("Baseline=%.3f | Plateau k=%d", base, k_pl),
           x = "Number of items added (bridge-first)",
           y = "Performance") +
      ylim(0.4, 1.05) +
      theme_classic(base_size = 11) +
      theme(plot.title    = element_text(face = "bold", size = 11),
            plot.subtitle = element_text(size = 9),
            axis.title    = element_text(size = 10))

    rep_plot_list[[paste(ctry, sname, sep="_")]] <- g
  }
}

cat(sprintf("Plots generated: %d\n", length(rep_plot_list)))
## Plots generated: 12
if (length(rep_plot_list) > 0) {
  gridExtra::grid.arrange(
    grobs = rep_plot_list,
    nrow  = length(rep_countries),
    ncol  = length(names(item_structure)),
    top   = grid::textGrob(
      "Representative Countries — Cross-construct Sequential Addition\n(rows = countries, columns = scales)",
      gp = grid::gpar(fontsize = 13, fontface = "bold")))
}
Figure 9. Sequential addition curves for representative countries (top-4 by bridge item consistency index). Rows = countries; columns = scales. Solid = AUROC; dashed = Kappa.

Figure 9. Sequential addition curves for representative countries (top-4 by bridge item consistency index). Rows = countries; columns = scales. Solid = AUROC; dashed = Kappa.

14. Boundary Conditions

boundary_df <- consistency_df %>%
  left_join(
    data %>% group_by(CTY) %>%
      summarise(SLM_High_pct = mean(BSDGSLM=="High", na.rm=TRUE),
                SCM_High_pct = mean(BSDGSCM=="High", na.rm=TRUE),
                SVM_High_pct = mean(BSDGSVM=="High", na.rm=TRUE),
                .groups = "drop"),
    by = c("Country" = "CTY"))

ggplot(boundary_df, aes(SLM_High_pct, Consistency, color = Region, label = Country, size = N)) +
  geom_point(alpha = 0.85) +
  ggrepel::geom_text_repel(size = 2.8, max.overlaps = 20) +
  geom_smooth(method = "lm", se = TRUE, color = "gray40",
              linewidth = 1, linetype = "dashed",
              inherit.aes = FALSE,
              aes(x = SLM_High_pct, y = Consistency)) +
  scale_color_manual(values = region_colors) +
  scale_size_continuous(range = c(2, 7)) +
  labs(title    = "Boundary Conditions: SLM High% vs. Bridge Consistency",
       subtitle = "Countries with lower High% tend to show lower consistency",
       x = "Proportion SLM = High",
       y = "Bridge Consistency Index") +
  theme_classic(base_size = 11) +
  theme(legend.position = "bottom",
        plot.title      = element_text(face = "bold"))
Figure 10. Relationship between SLM High proportion and bridge item consistency index by country.

Figure 10. Relationship between SLM High proportion and bridge item consistency index by country. 

boundary_df %>%
  group_by(Cluster) %>%
  summarise(
    N_countries  = n(),
    Mean_consist = round(mean(Consistency, na.rm=TRUE), 3),
    Mean_SLM_AU  = round(mean(SLM_AUROC,  na.rm=TRUE), 3),
    Mean_SCM_AU  = round(mean(SCM_AUROC,  na.rm=TRUE), 3),
    Mean_SVM_AU  = round(mean(SVM_AUROC,  na.rm=TRUE), 3),
    .groups = "drop"
  ) %>%
  kable(caption = "Table 8. Boundary Conditions by Cluster") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE)
Table 8. Boundary Conditions by Cluster
Cluster N_countries Mean_consist Mean_SLM_AU Mean_SCM_AU Mean_SVM_AU
C1 30 0.533 0.936 0.892 0.884
C2 4 0.550 0.945 0.899 0.863
C3 8 0.400 0.921 0.825 0.879
C4 2 0.500 0.921 0.839 0.825

Supplementary: Full Country Results

supp_table <- consistency_df %>%
  select(Country, Region, N, Consistency, Cluster,
         SLM_AUROC, SLM_k, SCM_AUROC, SCM_k,
         SVM_AUROC, SVM_k) %>%
  arrange(Region, Country)

kable(supp_table,
      caption = "Table S1. Full Results — All Countries") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed","responsive"),
                full_width = FALSE) %>%
  scroll_box(height = "500px")
Table S1. Full Results — All Countries
Country Region N Consistency Cluster SLM_AUROC SLM_k SCM_AUROC SCM_k SVM_AUROC SVM_k
BRA Americas 15077 0.4 C1 0.946 4 0.869 5 0.876 2
CHL Americas 3532 0.6 C1 0.922 5 0.915 6 0.892 2
USA Americas 6572 0.8 C1 0.952 5 0.900 7 0.886 4
HKG East Asia 3885 0.4 C1 0.945 5 0.875 5 0.915 2
JPN East Asia 3210 0.4 C1 0.966 8 0.938 4 0.867 2
KOR East Asia 3916 0.8 C1 0.927 3 0.920 2 0.902 4
SGP East Asia 4293 0.4 C1 0.944 7 0.920 2 0.885 3
TWN East Asia 4910 0.6 C1 0.948 5 0.972 4 0.895 3
ASH Europe 4929 0.6 C1 0.898 5 0.858 3 0.866 1
AUS Europe 6815 0.6 C1 0.952 6 0.903 3 0.916 4
AUT Europe 3877 0.6 C2 0.959 6 0.903 7 0.829 3
CYP Europe 3144 0.4 C1 0.910 1 0.869 3 0.840 2
CZE Europe 6915 0.8 C2 0.931 5 0.886 4 0.880 4
ENG Europe 3340 0.6 C1 0.904 4 0.922 5 0.915 3
FIN Europe 4964 0.6 C1 0.951 3 0.924 4 0.864 5
FRA Europe 3453 0.6 C1 0.952 4 0.907 3 0.873 5
GEO Europe 3382 0.4 C1 0.945 5 0.833 3 0.861 6
HUN Europe 4604 0.4 C2 0.929 7 0.890 2 0.892 5
IRL Europe 4309 0.6 C1 0.931 3 0.918 8 0.894 2
ITA Europe 3821 0.4 C2 0.960 4 0.917 5 0.852 2
LTU Europe 4352 0.6 C1 0.970 7 0.918 4 0.845 4
MLT Europe 3207 0.6 C1 0.947 5 0.895 6 0.880 2
NOR Europe 4951 0.6 C1 0.953 5 0.870 4 0.887 2
NZL Europe 3115 0.8 C1 0.952 7 0.899 10 0.864 2
PRT Europe 4526 0.6 C1 0.905 5 0.943 8 0.891 2
ROM Europe 2706 0.2 C1 0.930 4 0.907 5 0.902 2
SWE Europe 3948 0.6 C1 0.958 7 0.947 4 0.879 2
ARE Middle East/Africa 24557 0.6 C1 0.959 5 0.821 7 0.907 2
BHR Middle East/Africa 3721 0.4 C3 0.933 3 0.897 9 0.892 2
CIV Middle East/Africa 6959 0.4 C4 0.907 3 0.848 5 0.810 4
IRN Middle East/Africa 5149 0.4 C1 0.914 3 0.852 7 0.884 2
JOR Middle East/Africa 4490 0.6 C3 0.937 9 0.823 2 0.903 3
KWT Middle East/Africa 3174 0.4 C3 0.918 2 0.834 4 0.840 5
MAR Middle East/Africa 5600 0.6 C3 0.900 9 0.785 6 0.847 2
OMN Middle East/Africa 5152 0.4 C3 0.921 3 0.782 4 0.830 2
PSE Middle East/Africa 3995 0.4 C3 0.906 4 0.884 6 0.918 3
QAT Middle East/Africa 3368 0.4 C1 0.936 4 0.820 6 0.894 2
SAU Middle East/Africa 3895 0.2 C3 0.936 5 0.806 2 0.892 4
ZAF Middle East/Africa 9095 0.6 C4 0.935 6 0.831 4 0.840 2
AZE Other 3909 0.6 C1 0.914 1 0.870 3 0.864 2
ISR Other 4327 0.6 C1 0.894 2 0.827 5 0.882 3
KAZ Other 4826 0.4 C1 0.945 6 0.871 2 0.921 5
TUR Other 4150 0.2 C1 0.910 3 0.889 2 0.881 4
UZB Other 3164 0.2 C3 0.917 3 0.790 6 0.912 3 
write.csv(supp_table,
          "supplementary_table_S1_country_results.csv",
          row.names = FALSE)
cat("Saved: supplementary_table_S1_country_results.csv\n")
## Saved: supplementary_table_S1_country_results.csv

Session Info

sessionInfo()
## R version 4.6.0 (2026-04-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=Korean_Korea.utf8  LC_CTYPE=Korean_Korea.utf8   
## [3] LC_MONETARY=Korean_Korea.utf8 LC_NUMERIC=C                 
## [5] LC_TIME=Korean_Korea.utf8    
## 
## time zone: Asia/Seoul
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] kableExtra_1.4.0  knitr_1.51        psych_2.6.5       dendextend_1.19.1
##  [5] factoextra_2.0.0  ggrepel_0.9.8     gridExtra_2.3     caret_7.0-1      
##  [9] lattice_0.22-9    pROC_1.19.0.1     nnet_7.3-20       igraph_2.3.1     
## [13] bootnet_1.8       ggplot2_4.0.3     qgraph_1.9.8      purrr_1.2.2      
## [17] tidyr_1.3.2       dplyr_1.2.1       haven_2.5.5      
## 
## loaded via a namespace (and not attached):
##   [1] mathjaxr_2.0-0       RColorBrewer_1.1-3   rstudioapi_0.18.0   
##   [4] jsonlite_2.0.0       shape_1.4.6.1        magrittr_2.0.5      
##   [7] jomo_2.7-6           farver_2.1.2         nloptr_2.2.1        
##  [10] rmarkdown_2.31       vctrs_0.7.3          minqa_1.2.8         
##  [13] base64enc_0.1-6      htmltools_0.5.9      forcats_1.0.1       
##  [16] polynom_1.4-1        plotrix_3.8-14       weights_1.1.2       
##  [19] broom_1.0.13         Formula_1.2-5        mitml_0.4-5         
##  [22] parallelly_1.47.0    sass_0.4.10          bslib_0.10.0        
##  [25] htmlwidgets_1.6.4    plyr_1.8.9           lubridate_1.9.5     
##  [28] cachem_1.1.0         lifecycle_1.0.5      iterators_1.0.14    
##  [31] pkgconfig_2.0.3      Matrix_1.7-5         R6_2.6.1            
##  [34] fastmap_1.2.0        future_1.70.0        rbibutils_2.4.1     
##  [37] digest_0.6.39        fdrtool_1.2.18       colorspace_2.1-2    
##  [40] textshaping_1.0.5    Hmisc_5.2-5          ellipse_0.5.0       
##  [43] labeling_0.4.3       timechange_0.4.0     nnls_1.6            
##  [46] gdata_3.0.1          mgcv_1.9-4           abind_1.4-8         
##  [49] IsingSampler_0.2.4   compiler_4.6.0       proxy_0.4-29        
##  [52] withr_3.0.2          doParallel_1.0.17    glasso_1.11         
##  [55] htmlTable_2.5.0      S7_0.2.2             backports_1.5.1     
##  [58] viridis_0.6.5        mgm_1.2-15           R.utils_2.13.0      
##  [61] pan_1.9              lava_1.9.1           MASS_7.3-65         
##  [64] corpcor_1.6.10       gtools_3.9.5         ModelMetrics_1.2.2.2
##  [67] tools_4.6.0          pbivnorm_0.6.0       foreign_0.8-91      
##  [70] otel_0.2.0           future.apply_1.20.2  R.oo_1.27.1         
##  [73] glue_1.8.1           quadprog_1.5-8       NetworkToolbox_1.4.4
##  [76] nlme_3.1-169         grid_4.6.0           checkmate_2.3.4     
##  [79] cluster_2.1.8.2      reshape2_1.4.5       generics_0.1.4      
##  [82] snow_0.4-4           recipes_1.3.2        gtable_0.3.6        
##  [85] R.methodsS3_1.8.2    class_7.3-23         data.table_1.18.4   
##  [88] hms_1.1.4            xml2_1.5.2           foreach_1.5.2       
##  [91] pillar_1.11.1        stringr_1.6.0        splines_4.6.0       
##  [94] smacof_2.1-7         networktools_1.6.0   survival_3.8-6      
##  [97] tidyselect_1.2.1     pbapply_1.7-4        reformulas_0.4.4    
## [100] svglite_2.2.2        IsingFit_0.4         stats4_4.6.0        
## [103] xfun_0.57            hardhat_1.4.3        timeDate_4052.112   
## [106] stringi_1.8.7        yaml_2.3.12          boot_1.3-32         
## [109] evaluate_1.0.5       codetools_0.2-20     wordcloud_2.6       
## [112] tibble_3.3.1         cli_3.6.6            mantar_0.2.0        
## [115] rpart_4.1.27         systemfonts_1.3.2    Rdpack_2.6.6        
## [118] jquerylib_0.1.4      lavaan_0.6-21        Rcpp_1.1.1-1.1      
## [121] globals_0.19.1       png_0.1-9            parallel_4.6.0      
## [124] gower_1.0.2          jpeg_0.1-11          listenv_0.10.1      
## [127] lme4_2.0-1           glmnet_5.0           viridisLite_0.4.3   
## [130] mvtnorm_1.3-7        ipred_0.9-15         prodlim_2026.03.11  
## [133] scales_1.4.0         e1071_1.7-17         eigenmodel_1.12     
## [136] rlang_1.2.0          mnormt_2.1.2         mice_3.19.0