# False-positive rate of count models for the page's design: 6 vs 6 animals (and 20 vs 20 for contrast),
# one count per animal, negative binomial counts with the same mean in both groups (no true difference).
# Seed fixed in advance (not searched). Run: Rscript count_simulation.R
suppressPackageStartupMessages(library(MASS))
set.seed(20261110)
runs <- 5000
grid <- expand.grid(n = c(6, 20), mu = c(5, 10), size = c(2, 5))
res <- do.call(rbind, lapply(seq_len(nrow(grid)), function(i) {
  g <- grid[i, ]
  grp <- factor(rep(c("vehicle", "drug"), each = g$n), levels = c("vehicle", "drug"))
  rej <- replicate(runs, {
    y <- rnbinom(2 * g$n, mu = g$mu, size = g$size)
    qp <- summary(glm(y ~ grp, family = quasipoisson))$coefficients[2, 4]
    nb <- tryCatch(suppressWarnings(glm.nb(y ~ grp)), error = function(e) NULL)
    if (is.null(nb)) return(c(qp < 0.05, NA, NA))
    wald <- summary(nb)$coefficients[2, 4]
    nb0 <- tryCatch(suppressWarnings(glm.nb(y ~ 1)), error = function(e) NULL)
    lr <- if (is.null(nb0)) NA else pchisq(as.numeric(2 * (logLik(nb) - logLik(nb0))), 1, lower.tail = FALSE)
    c(qp < 0.05, wald < 0.05, lr < 0.05)
  })
  data.frame(n_per_group = g$n, mean = g$mu, nb_size = g$size,
             quasi_poisson_t = mean(rej[1, ]), nb_wald = mean(rej[2, ], na.rm = TRUE),
             nb_lr = mean(rej[3, ], na.rm = TRUE))
}))
print(res, row.names = FALSE, digits = 3)