# R replication of simulation.py for /guides/should-you-test-for-normality/ (R 4.3.3, base stats only).
# Independent random numbers (R's own generator, seed 20261110), so results agree with Python within
# simulation error (about +/- 0.3 percentage points at 5% with 20,000 runs), not digit for digit.
# Run: Rscript simulation.R > simulation_R_output.txt
set.seed(20261110)
REPS <- 20000
ALPHA <- 0.05

draw <- function(dist, n) switch(dist,
  "normal" = rnorm(n),
  "lognormal (sigma 0.5)" = rlnorm(n, 0, 0.5),
  "t, 5 df (heavier tails)" = rt(n, 5),
  "lognormal (sigma 0.1, nearly normal)" = rlnorm(n, 0, 0.1))

cat("Shapiro-Wilk flag rate (alpha 0.05)\n")
for (dist in c("normal", "lognormal (sigma 0.5)", "t, 5 df (heavier tails)", "lognormal (sigma 0.1, nearly normal)")) {
  for (n in c(5, 10, 30, 100, 1000)) {
    hits <- sum(replicate(REPS, shapiro.test(draw(dist, n))$p.value < ALPHA))
    cat(sprintf("%-40s n = %4d  flagged %.4f\n", dist, n, hits / REPS))
  }
}

scen <- list(
  list("normal, equal SD, 5 vs 5", 5, 5, 1, 1, "normal"),
  list("normal, equal SD, 10 vs 10", 10, 10, 1, 1, "normal"),
  list("lognormal, same distribution, 10 vs 10", 10, 10, 1, 1, "lognormal"),
  list("normal, SD 2 vs 1, 10 vs 10", 10, 10, 2, 1, "normal"),
  list("normal, SD 2 vs 1, 5 vs 15 (smaller group more variable)", 5, 15, 2, 1, "normal"))
grp <- function(n, sd, shape) if (shape == "normal") rnorm(n, 0, sd) else rlnorm(n, 0, 0.5)

cat("\nTwo-stage (Shapiro -> Student t or Mann-Whitney) vs always Welch; no true difference\n")
for (s in scen) {
  res <- replicate(REPS, {
    a <- grp(s[[2]], s[[4]], s[[6]]); b <- grp(s[[3]], s[[5]], s[[6]])
    pass <- shapiro.test(a)$p.value >= ALPHA && shapiro.test(b)$p.value >= ALPHA
    p2 <- if (pass) t.test(a, b, var.equal = TRUE)$p.value else wilcox.test(a, b)$p.value
    c(pass, p2 < ALPHA, pass && p2 < ALPHA, t.test(a, b)$p.value < ALPHA,
      t.test(a, b, var.equal = TRUE)$p.value < ALPHA)
  })
  cat(sprintf("%-58s sent to t %.4f | two-stage FP %.4f | t FP when chosen %.4f | Welch FP %.4f | Student FP %.4f\n",
              s[[1]], mean(res[1, ]), mean(res[2, ]), sum(res[3, ]) / sum(res[1, ]), mean(res[4, ]), mean(res[5, ])))
}