# What is ‘item sampling’ and why is it useful? Or how can we can measure person characteristics with large item pools?
# R code from https://philippmasur.de/blog/2024-06-19-item-sampling/

# Set parameters
alpha <- .9
k_s <- 5
k <- 10

# Formula from above
new_alpha <- ((k_s/k)*alpha)/(1+((k_s/k)-1)*alpha)
round(new_alpha, 2)

library(tidyverse)

# Custom function to compute reliability and correlation
reli_new <- function(alpha, k_s, k, N = 100) {

  # Ratio between sample items and total items
  n <- k_s/k

  # Spearman-Brown formula
  alpha_new <- ((n)*alpha)/(1+((n)-1)*alpha)

  # Standard error of alpha_new (van Zyl, Neudecker & Nel, 2000)
  se_a <- (1 - alpha_new) * sqrt(2 * k_s / ((k_s - 1) * N))

  # Correlation between scales (part-whole correlation for parallel items)
  r <- sqrt(alpha_new/alpha)
  se_r <- (1-r^2)/sqrt(N)

  # Create output table
  tibble(N = N,
         k = k,
         k_s = k_s,
         n = n,
         alpha = alpha,
         alpha_new = alpha_new,
         se_a = se_a,
         r = r,
         se_r = se_r)
}

# Check function
(check <- reli_new(alpha = .9, k_s = 5, k = 10, N = 100))

# Transparent background and muted text, so the figures work in light and dark mode
theme_blog <- function() {
  theme_minimal(base_size = 12) +
    theme(plot.background  = element_rect(fill = "transparent", colour = NA),
          panel.background = element_rect(fill = "transparent", colour = NA),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_line(colour = "#8a8f9840", linewidth = .4),
          strip.text = element_text(colour = "#7a7f88", face = "bold"),
          axis.text  = element_text(colour = "#7a7f88"),
          axis.title = element_text(colour = "#7a7f88"),
          legend.text  = element_text(colour = "#7a7f88"),
          legend.title = element_text(colour = "#7a7f88"),
          legend.position = "bottom")
}
blues <- c("#c4e4f6", "#8cc8ec", "#4aa8e0", "#1680bd", "#0b5a87")   # sequential: alpha .75 -> .95

# Expected reliability for a 20-item pool across sampling ratios and original alphas
analytic <- expand.grid(alpha = seq(.75, .95, by = .05), k_s = 2:20, N = c(100, 1000)) |>
  pmap_dfr(reli_new, k = 20) |>
  mutate(N = factor(paste("N =", format(N, big.mark = ",")), levels = c("N = 100", "N = 1,000")))

ggplot(analytic, aes(n, alpha_new, colour = factor(alpha), fill = factor(alpha))) +
  geom_hline(yintercept = .8, linetype = "dashed", colour = "#8a8f98") +
  geom_ribbon(aes(ymin = alpha_new - 1.96 * se_a, ymax = alpha_new + 1.96 * se_a),
              alpha = .25, colour = NA) +
  geom_line(linewidth = 1) +
  facet_wrap(~ N) +
  scale_colour_manual(values = blues) +
  scale_fill_manual(values = blues) +
  scale_x_continuous(labels = scales::percent) +
  labs(x = "Share of items each person answers (k_s / k)",
       y = "Expected reliability (alpha)",
       colour = "Alpha of the full scale", fill = "Alpha of the full scale") +
  theme_blog()

ggplot(analytic, aes(n, r, colour = factor(alpha), fill = factor(alpha))) +
  geom_hline(yintercept = c(.8, .9), linetype = "dashed", colour = "#8a8f98") +
  geom_ribbon(aes(ymin = r - 1.96 * se_r, ymax = pmin(1, r + 1.96 * se_r)), alpha = .25, colour = NA) +
  geom_line(linewidth = 1) +
  facet_wrap(~ N) +
  scale_colour_manual(values = blues) +
  scale_fill_manual(values = blues) +
  scale_x_continuous(labels = scales::percent) +
  labs(x = "Share of items each person answers (k_s / k)",
       y = "Expected correlation with the full scale",
       colour = "Alpha of the full scale", fill = "Alpha of the full scale") +
  theme_blog()

library(lavaan)
library(mice)
library(furrr)

set.seed(123)  # for reproducibility

sim_items <- function(N, k, alpha, answer_options = 5) {

  # Average inter-item correlation implied by alpha, and the resulting loading
  rbar   <- alpha / (k - alpha * (k - 1))
  lambda <- sqrt(rbar)

  # True scores, item difficulties and continuous responses
  theta      <- rnorm(N)
  difficulty <- runif(k, -.75, .75)
  latent <- theta %o% rep(lambda, k) +
    matrix(rnorm(N * k, sd = sqrt(1 - lambda^2)), N, k) -
    matrix(difficulty, N, k, byrow = TRUE)

  # Cut into Likert-type answer options (1 to 5)
  cuts  <- qnorm(seq(0, 1, length.out = answer_options + 1))[2:answer_options]
  items <- matrix(findInterval(latent, cuts) + 1, N, k,
                  dimnames = list(NULL, paste0("i", 1:k)))

  list(items = as.data.frame(items), theta = theta)
}

# Check function
test <- sim_items(N = 100, k = 10, alpha = .9)
head(test$items)

psych::alpha(test$items)$total

sample_items <- function(items, n) {
  k   <- ncol(items)
  k_s <- round(n * k)
  for (i in seq_len(nrow(items))) {
    items[i, sample(k, k - k_s)] <- NA
  }
  items
}

# Check function
test_na <- sample_items(test$items, n = .8)
head(test_na)

mean(is.na(test_na))

score_mean <- function(items) rowMeans(items, na.rm = TRUE)

score_adjusted <- function(items) {
  centered <- sweep(as.matrix(items), 2, colMeans(items, na.rm = TRUE))
  rowMeans(centered, na.rm = TRUE)
}

score_cfa <- function(items) {
  model <- paste("G =~", paste(names(items), collapse = " + "))
  scores <- as.numeric(lavPredict(cfa(model, data = items, std.lv = TRUE)))
  # The sign of a factor is arbitrary: align it with the sum score
  # (all items are positively keyed)
  if (cor(scores, rowSums(items)) < 0) scores <- -scores
  scores
}

score_imputed <- function(items) {
  score_cfa(complete(mice(items, m = 1, printFlag = FALSE)))
}

# Reliability of each score = squared correlation with the true score
sapply(list(mean     = score_mean(test_na),
            adjusted = score_adjusted(test_na),
            imputed  = score_imputed(test_na)),
       \(score) cor(score, test$theta)^2) |> round(2)

alpha_full <- function(items) {
  k <- ncol(items)
  k / (k - 1) * (1 - sum(apply(items, 2, var)) / var(rowSums(items)))
}

sim_func <- function(N, k, alpha, n) {
  d <- sim_items(N, k, alpha)
  sampled <- sample_items(d$items, n)
  a <- alpha_full(d$items)

  full <- list(mean = score_mean(d$items), adjusted = score_mean(d$items),
               imputed = score_cfa(d$items))
  part <- list(mean = score_mean(sampled), adjusted = score_adjusted(sampled),
               imputed = score_imputed(sampled))

  tibble(N = N, k = k, k_s = round(n * k), n = n, alpha = alpha, alpha_emp = a,
         method = names(part),
         r_full = mapply(cor, part, full),
         reliability = sapply(part, cor, y = d$theta)^2,
         reliability_full = sapply(full, cor, y = d$theta)^2,
         spearman_brown = n * a / (1 + (n - 1) * a))
}

# Check sim function
sim_func(N = 500, k = 20, alpha = .85, n = .5)

# Setup simulation grid
setup <- expand.grid(
  N = c(250, 500, 1000),
  k = c(10, 20, 40),
  alpha = seq(.75, .95, by = .05),
  n = c(.3, .4, .5, .6, .7, .8)
)

# The simulation takes a while: its results are saved next to this post and
# re-used (delete sims/results.rds to run it again)
n_sim <- 200
if (file.exists("sims/results.rds")) {
  results <- readRDS("sims/results.rds")
} else {
  plan(multisession, workers = parallelly::availableCores(omit = 1))
  results <- setup[rep(seq_len(nrow(setup)), n_sim), ] |>
    future_pmap_dfr(sim_func, .options = furrr_options(seed = 2024))
  dir.create("sims", showWarnings = FALSE)
  saveRDS(results, "sims/results.rds")
}

results

method_labels <- c(adjusted = "Difficulty-adjusted mean", mean = "Mean of answered items",
                   imputed = "Imputation (PMM) + factor scores")
method_cols <- c(adjusted = "#1497d6", mean = "#7fc3ec", imputed = "#d9772b")

summ <- results |>
  group_by(N, k, alpha, n, method) |>
  summarise(rel = mean(reliability), rel_lo = quantile(reliability, .025), rel_hi = quantile(reliability, .975),
            r = mean(r_full), r_lo = quantile(r_full, .025), r_hi = quantile(r_full, .975),
            sb = mean(spearman_brown), alpha_emp = mean(alpha_emp), .groups = "drop") |>
  mutate(method = factor(method, levels = names(method_labels)))

# Panel labels: input alpha and the (slightly lower) alpha of the simulated Likert data
alpha_labs <- summ |> group_by(alpha) |> summarise(a = mean(alpha_emp)) |>
  mutate(lab = sprintf("alpha = %.2f\n(Likert data: %.2f)", alpha, a))
summ <- left_join(summ, select(alpha_labs, alpha, lab), by = "alpha")

# Numbers used in the text
ex <- filter(summ, N == 500, k == 20, alpha == .85, n == .5)
adj <- filter(summ, method == "adjusted")
cell <- function(a, nn, what = "rel") f2(filter(adj, N == 500, k == 20, alpha == a, n == nn)[[what]])
loss_by_n <- filter(summ, N == 500, k == 20) |>
  group_by(n) |> summarise(loss = mean(rel[method == "adjusted"] - rel[method == "imputed"]))
gap_sb <- summ |> group_by(method) |> summarise(diff = mean(rel - sb))
loss_imp <- with(filter(summ, N == 500, k == 20),
                 mean(rel[method == "adjusted"] - rel[method == "imputed"]))
f2 <- function(x) sprintf("%.2f", x)

d <- filter(summ, N == 500, k == 20)

ggplot(d, aes(n, rel, colour = method)) +
  geom_line(aes(y = sb), colour = "#8a8f98", linetype = "dashed", linewidth = .8) +
  geom_pointrange(aes(ymin = rel_lo, ymax = rel_hi), size = .2,
                  position = position_dodge(width = .04)) +
  geom_line(position = position_dodge(width = .04)) +
  facet_wrap(~ lab, nrow = 1) +
  scale_colour_manual(values = method_cols, labels = method_labels) +
  scale_x_continuous(labels = scales::percent, breaks = c(.3, .5, .7)) +
  labs(x = "Share of items each person answers",
       y = "Reliability of the person scores",
       colour = NULL,
       caption = "Points: mean of 200 replications, bars: 95% of replications. Dashed: Spearman–Brown prediction.") +
  theme_blog() +
  theme(plot.caption = element_text(colour = "#7a7f88"))

ggplot(d, aes(n, r, colour = method)) +
  geom_line(aes(y = sqrt(sb / alpha_emp)), colour = "#8a8f98", linetype = "dashed", linewidth = .8) +
  geom_hline(yintercept = c(.8, .9), colour = "#8a8f9860") +
  geom_pointrange(aes(ymin = r_lo, ymax = r_hi), size = .2,
                  position = position_dodge(width = .04)) +
  geom_line(position = position_dodge(width = .04)) +
  facet_wrap(~ lab, nrow = 1) +
  scale_colour_manual(values = method_cols, labels = method_labels) +
  scale_x_continuous(labels = scales::percent, breaks = c(.3, .5, .7)) +
  labs(x = "Share of items each person answers",
       y = "Correlation with the full-scale score",
       colour = NULL) +
  theme_blog()

summ |>
  filter(method == "adjusted", alpha == .85, n == .5) |>
  ggplot(aes(factor(N), rel, colour = factor(k))) +
  geom_hline(aes(yintercept = mean(sb)), colour = "#8a8f98", linetype = "dashed") +
  geom_pointrange(aes(ymin = rel_lo, ymax = rel_hi), position = position_dodge(width = .5)) +
  scale_colour_manual(values = c("10" = "#8cc8ec", "20" = "#1497d6", "40" = "#0b5a87")) +
  labs(x = "Sample size", y = "Reliability", colour = "Items in the pool") +
  theme_blog()

summ |>
  filter(N == 500, k == 20, n == .6, alpha == .8) |>
  mutate(method = method_labels[as.character(method)]) |>
  select(method, alpha_emp, `Spearman-Brown` = sb, reliability = rel, r_full = r)

# Mathematical results (with the alpha of the simulated Likert-type data)
a <- unique(round(filter(summ, N == 500, k == 20, n == .6, alpha == .8)$alpha_emp, 3))
reli_new(alpha = a, k_s = 12, k = 20, N = 500) |>
  mutate(ll = r - 1.96*se_r,
         ul = r + 1.96*se_r) |>
  dplyr::select(N, k, k_s, alpha, alpha_new, r, ll, ul)

sessionInfo()