# DASS Analysis Clinic, Case 002: "The reviewer wants a multiple-comparisons correction"
# https://statisticalsolutions.org/analysis-clinic/multiple-comparisons
#
# Base R only. How fast the chance of a false positive grows, and what
# Bonferroni, Holm, and Benjamini-Hochberg do to the same set of p-values.

# --- 1. Familywise error with independent tests under the null ---------------
k <- c(1, 3, 5, 10, 20)
print(data.frame(tests = k, P_at_least_one_false_positive = round(1 - 0.95^k, 3)))

# Check by simulation: 10 independent null tests, 20,000 replications.
set.seed(2026)
sims <- replicate(20000, any(runif(10) < 0.05))
cat(sprintf("Simulated P(at least one p < .05 among 10 null tests): %.3f\n", mean(sims)))

# --- 2. The same p-values under three corrections -----------------------------
# Example: one primary outcome and eight secondary outcomes from a study.
p <- c(primary = 0.004,
       sec_1 = 0.011, sec_2 = 0.019, sec_3 = 0.032, sec_4 = 0.041,
       sec_5 = 0.090, sec_6 = 0.210, sec_7 = 0.480, sec_8 = 0.730)

# If the reviewer's "family" is all nine tests:
all9 <- data.frame(raw = p,
                   bonferroni = p.adjust(p, "bonferroni"),
                   holm = p.adjust(p, "holm"),
                   BH_fdr = p.adjust(p, "BH"))
print(round(all9, 3))
cat("Significant at .05 -> raw:", sum(p < .05),
    " Bonferroni:", sum(all9$bonferroni < .05),
    " Holm:", sum(all9$holm < .05),
    " BH:", sum(all9$BH_fdr < .05), "\n")

# If the prespecified family is the primary outcome alone, it is tested at .05,
# and the secondary family gets its own correction:
sec <- p[-1]
print(round(data.frame(raw = sec, holm = p.adjust(sec, "holm"), BH_fdr = p.adjust(sec, "BH")), 3))