# Review check (2026-11-10): does the negative binomial likelihood-ratio statistic referred to F(1, residual df)
# instead of chi-square(1) hold its false-positive rate at small n? Same data model as simulation_n.R
# (theta = 2, areas uniform 0.7 to 1.6 mm^2, no true difference). Seed fixed in advance.
# R 4.3.3, MASS 7.3-60. Run: Rscript simulation_lrt_f.R > simulation_lrt_f_R_output.txt
suppressPackageStartupMessages(library(MASS))
set.seed(20261112)
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)
  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)
  if (is.null(nb) || is.null(nb0)) return(c(chisq = NA, f = NA))
  lr <- max(0, 2 * (as.numeric(logLik(nb)) - as.numeric(logLik(nb0))))
  c(chisq = pchisq(lr, 1, lower.tail = FALSE), f = pf(lr, 1, nb$df.residual, lower.tail = FALSE))
}

cat(sprintf("runs per row: %d; theta 2; no true difference; seed 20261112\n", runs))
cat(sprintf("%-6s %-6s %14s %14s\n", "n", "rate", "LR vs chi2(1)", "LR vs F(1,df)"))
for (n in c(6, 10)) for (rate in c(3, 15)) {
  r <- 100 * colMeans(t(replicate(runs, one_run(n, rate))) < 0.05, na.rm = TRUE)
  cat(sprintf("%-6d %-6g %13.1f%% %13.1f%%\n", n, rate, r[1], r[2]))
}