# Follow-up to simulation.R: does the negative binomial Wald test's false-positive rate settle as groups grow?
# Same data model (areas uniform 0.7 to 1.6 mm^2; negative binomial counts with theta = 2, i.e. overdispersed),
# no true difference, n = 6, 10, 20 and 40 per group. Seed fixed in advance.
# R 4.3.3, MASS 7.3-60. Run: Rscript simulation_n.R > simulation_n_R_output.txt
suppressPackageStartupMessages(library(MASS))
set.seed(20261111)
runs <- 5000

one_run <- function(n, rate, theta = 2) {
  grp <- factor(rep(c("vehicle", "drug"), each = n), levels = c("vehicle", "drug"))
  area <- runif(2 * n, 0.7, 1.6)
  y <- rnbinom(2 * n, mu = rate * area, size = theta)
  off <- log(area)
  qp <- glm(y ~ grp + offset(off), family = quasipoisson)
  nb <- tryCatch(suppressWarnings(glm.nb(y ~ grp + offset(off))), error = function(e) NULL)
  nb0 <- tryCatch(suppressWarnings(glm.nb(y ~ 1 + offset(off))), error = function(e) NULL)
  c(quasipoisson = summary(qp)$coefficients[2, 4],
    nb_wald = if (is.null(nb)) NA else summary(nb)$coefficients[2, 4],
    nb_lrt = if (is.null(nb) || is.null(nb0)) NA else
      pchisq(max(0, 2 * (as.numeric(logLik(nb)) - as.numeric(logLik(nb0)))), 1, lower.tail = FALSE))
}

cat(sprintf("runs per row: %d; theta 2; no true difference; seed 20261111\n", runs))
cat(sprintf("%-6s %-6s %13s %9s %8s\n", "n", "rate", "quasipoisson", "nb_wald", "nb_lrt"))
for (n in c(6, 10, 20, 40)) for (rate in c(3, 15)) {
  res <- t(replicate(runs, one_run(n, rate)))
  r <- 100 * colMeans(res < 0.05, na.rm = TRUE)
  cat(sprintf("%-6d %-6g %12.1f%% %8.1f%% %7.1f%%\n", n, rate, r[1], r[2], r[3]))
}