# DASS Analysis Clinic, Case 001: "Reviewer 2 asked for a post hoc power analysis"
# https://statisticalsolutions.org/analysis-clinic/post-hoc-power
#
# Base R only. Shows that "observed power" is just the p-value in another unit,
# then runs the analyses that actually answer the reviewer: a confidence
# interval, an equivalence test (TOST), and a design sensitivity analysis.

# --- 1. Observed power is a function of the p-value --------------------------
# For a two-sided z test at alpha = .05, plug the observed effect back in as if
# it were the true effect. The p-value alone determines the answer.
observed_power <- function(p, alpha = 0.05) {
  z_obs  <- qnorm(1 - p / 2)          # |z| that produced this p
  z_crit <- qnorm(1 - alpha / 2)
  pnorm(z_obs - z_crit) + pnorm(-z_obs - z_crit)
}
p_values <- c(0.01, 0.03, 0.05, 0.10, 0.20, 0.40, 0.80)
print(data.frame(p = p_values, observed_power = round(observed_power(p_values), 3)))
# p = .05 gives observed power of about .50; every nonsignificant p gives about .50 or less.

# --- 2. A worked example ------------------------------------------------------
# Two groups of 40, outcome on a 0-100 scale. Nonsignificant result.
set.seed(2026)
n <- 40
control   <- rnorm(n, mean = 50, sd = 12)
treatment <- rnorm(n, mean = 52, sd = 12)
res <- t.test(treatment, control, var.equal = TRUE)
print(res)
cat(sprintf("Observed power at this p (%.3f): %.2f\n", res$p.value, observed_power(res$p.value)))

# --- 3. What to report instead -----------------------------------------------
# 3a. The confidence interval: which differences are still compatible with the data?
ci <- res$conf.int
cat(sprintf("95%% CI for the difference: [%.2f, %.2f] points\n", ci[1], ci[2]))

# 3b. Equivalence test (two one-sided tests, TOST). Choose the smallest
# difference that would matter BEFORE looking at the result: here 6 points.
tost <- function(x, y, bound) {
  lower <- t.test(x, y, mu = -bound, alternative = "greater", var.equal = TRUE)$p.value
  upper <- t.test(x, y, mu =  bound, alternative = "less",    var.equal = TRUE)$p.value
  c(p_lower = lower, p_upper = upper, p_tost = max(lower, upper))
}
print(round(tost(treatment, control, bound = 6), 4))
# The 90% CI (alpha = .05 for each one-sided test) must sit inside [-6, 6].
ci90 <- t.test(treatment, control, var.equal = TRUE, conf.level = 0.90)$conf.int
cat(sprintf("90%% CI: [%.2f, %.2f]; equivalence at +/-6 %s\n", ci90[1], ci90[2],
            ifelse(ci90[1] > -6 && ci90[2] < 6, "supported", "not supported")))

# 3c. Design sensitivity: what standardized effect could n = 40 per group detect
# with 80% power? Describes the design; it is not evidence for the null.
sens <- power.t.test(n = n, power = 0.80, sig.level = 0.05, type = "two.sample")
cat(sprintf("With n = %d per group, 80%% power for d = %.2f or larger\n", n, sens$delta))