# Simulation for /guides/statistical-tests-for-count-data/ (constructed scenarios; seed fixed in advance).
# Two groups of 6 animals; each animal has its own observation area (uniform 0.7 to 1.6 mm^2).
# Counts are negative binomial with mean = rate x area. theta = 2 is clearly overdispersed;
# theta = 1000 is practically Poisson. For each run, every method is applied to the same data.
# Output: share of runs with p < 0.05 (false positives when the true rate ratio is 1; power when it is 2),
# and the median estimated rate ratio when the truth is 2.
# R 4.3.3, MASS 7.3-60. Run: Rscript simulation.R > simulation_R_output.txt
suppressPackageStartupMessages(library(MASS))
set.seed(20261110)
runs <- 10000
n <- 6
grp <- factor(rep(c("vehicle", "drug"), each = n), levels = c("vehicle", "drug"))

one_run <- function(rate, ratio, theta) {
  area <- runif(2 * n, 0.7, 1.6)
  mu <- rate * area * ifelse(grp == "drug", ratio, 1)
  y <- rnbinom(2 * n, mu = mu, size = theta)
  off <- log(area)
  out <- c(poisson = NA, quasipoisson = NA, nb_wald = NA, nb_lrt = NA, log1p_welch = NA, rate_welch = NA,
           rr_nb = NA, rr_log1p = NA, nb_failed = 0)
  po <- glm(y ~ grp + offset(off), family = poisson)
  out["poisson"] <- summary(po)$coefficients[2, 4]
  qp <- glm(y ~ grp + offset(off), family = quasipoisson)
  out["quasipoisson"] <- summary(qp)$coefficients[2, 4]
  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)) {
    out["nb_failed"] <- 1
  } else {
    out["nb_wald"] <- summary(nb)$coefficients[2, 4]
    lr <- 2 * (as.numeric(logLik(nb)) - as.numeric(logLik(nb0)))
    out["nb_lrt"] <- pchisq(max(lr, 0), df = 1, lower.tail = FALSE)
    out["rr_nb"] <- exp(coef(nb)[2])
  }
  l <- log(y + 1)
  out["log1p_welch"] <- if (sd(l[grp == "drug"]) + sd(l[grp == "vehicle"]) > 0)
    t.test(l[grp == "drug"], l[grp == "vehicle"])$p.value else 1
  r <- y / area
  out["rate_welch"] <- if (sd(r[grp == "drug"]) + sd(r[grp == "vehicle"]) > 0)
    t.test(r[grp == "drug"], r[grp == "vehicle"])$p.value else 1
  out["rr_log1p"] <- exp(mean(l[grp == "drug"]) - mean(l[grp == "vehicle"]))
  out
}

scenarios <- expand.grid(rate = c(3, 15), theta = c(2, 1000), ratio = c(1, 2))
methods <- c("poisson", "quasipoisson", "nb_wald", "nb_lrt", "log1p_welch", "rate_welch")
cat(sprintf("runs per scenario: %d; n = %d per group; seed 20261110\n", runs, n))
for (i in seq_len(nrow(scenarios))) {
  s <- scenarios[i, ]
  res <- t(replicate(runs, one_run(s$rate, s$ratio, s$theta)))
  rej <- colMeans(res[, methods] < 0.05, na.rm = TRUE)
  cat(sprintf("\nrate %g per mm2, theta %g, true ratio %g\n", s$rate, s$theta, s$ratio))
  cat(sprintf("  %-13s %5.1f%%\n", methods, 100 * rej), sep = "")
  cat(sprintf("  NB fits failed: %d\n", sum(res[, "nb_failed"])))
  if (s$ratio == 2) cat(sprintf("  median estimated ratio: NB %.2f; log(count+1) back-transformed %.2f\n",
                                median(res[, "rr_nb"], na.rm = TRUE), median(res[, "rr_log1p"])))
}