# Simulation for /guides/comparing-percentages-between-groups/.
# Share of experiments with p < 0.05 for four analyses of "positive cells out of cells counted" per mouse.
# Data model: cells counted per mouse uniform on 30..200; vehicle baseline 15% positive; mouse-to-mouse
# variation normal on the log-odds scale (SD given); drug odds ratio 1 (no effect) or 2.
# R 4.3.3, lme4 1.1-35.1. Seed fixed in advance: 20261110. 2,000 runs per row.
suppressPackageStartupMessages(library(lme4))
set.seed(20261110)
lo <- function(p) log(p / (1 - p)); inv <- function(x) 1 / (1 + exp(-x))
one <- function(m, sd, or) {
  g <- rep(c(0, 1), each = m)
  counted <- sample(30:200, 2 * m, replace = TRUE)
  p <- inv(lo(0.15) + log(or) * g + rnorm(2 * m, 0, sd))
  pos <- rbinom(2 * m, counted, p)
  d <- data.frame(g = factor(g), pos, counted, mouse = factor(seq_len(2 * m)))
  pooled <- suppressWarnings(chisq.test(rbind(c(sum(pos[g == 0]), sum(counted[g == 0] - pos[g == 0])),
                                              c(sum(pos[g == 1]), sum(counted[g == 1] - pos[g == 1]))),
                                        correct = FALSE)$p.value)
  pct <- pos / counted
  welch <- t.test(pct[g == 1], pct[g == 0])$p.value
  qb <- summary(glm(cbind(pos, counted - pos) ~ g, family = quasibinomial, data = d))$coefficients[2, 4]
  gm <- tryCatch(suppressWarnings(suppressMessages(
    summary(glmer(cbind(pos, counted - pos) ~ g + (1 | mouse), family = binomial, data = d))$coefficients[2, 4])),
    error = function(e) NA)
  c(pooled = pooled, welch = welch, quasibinomial = qb, glmm = gm)
}
rows <- list(c(5, 0, 1), c(5, 0.4, 1), c(5, 0.8, 1), c(3, 0.4, 1), c(8, 0.4, 1), c(5, 0.4, 2), c(8, 0.4, 2))
cat("mice_per_group mouse_sd odds_ratio | pooled_chisq welch_pct quasibinomial glmm_z (share p<0.05; runs=2000)\n")
for (r in rows) {
  res <- replicate(2000, one(r[1], r[2], r[3]))
  rate <- rowMeans(res < 0.05, na.rm = TRUE)
  cat(sprintf("%d %.1f %.0f | %.3f %.3f %.3f %.3f  (glmm NA: %d)\n", r[1], r[2], r[3],
              rate[1], rate[2], rate[3], rate[4], sum(is.na(res[4, ]))))
}