# DASS Analysis Clinic, Case 003: "Can I power my grant on my pilot's effect size?"
# https://statisticalsolutions.org/analysis-clinic/pilot-effect-size-power
#
# Base R only. Simulates the common plan: run a small pilot, estimate d,
# power the main study on that estimate, and go ahead only if the required
# sample is affordable. Then compares the power the main study actually has.

set.seed(2026)
true_d   <- 0.40   # the real effect, unknown to the researcher
n_pilot  <- 15     # per group
max_n    <- 150    # largest main-study sample per group the budget allows
n_sims   <- 5000

planned_n <- function(d, power = 0.80) {
  if (d <= 0.05) return(Inf)
  ceiling(power.t.test(delta = d, sd = 1, power = power, sig.level = 0.05)$n)
}
true_power <- function(n) power.t.test(n = n, delta = true_d, sd = 1, sig.level = 0.05)$power

one_pilot <- function() {
  x <- rnorm(n_pilot, true_d); y <- rnorm(n_pilot, 0)
  sp <- sqrt((var(x) + var(y)) / 2)
  d_hat <- (mean(x) - mean(y)) / sp
  # Safeguard effect size: lower limit of a 60% CI for d (Perugini et al., 2014)
  se_d <- sqrt(2 / n_pilot + d_hat^2 / (4 * n_pilot))
  d_safe <- d_hat - qnorm(0.80) * se_d
  c(d_hat = d_hat, d_safe = d_safe)
}
pilots <- t(replicate(n_sims, one_pilot()))

n_naive <- sapply(pilots[, "d_hat"], planned_n)
go      <- n_naive <= max_n                    # follow-up only if affordable
power_naive <- sapply(n_naive[go], true_power)

cat(sprintf("True d = %.2f; n needed for 80%% power = %d per group\n", true_d, planned_n(true_d)))
cat(sprintf("Pilot d estimates: middle 95%% from %.2f to %.2f\n",
            quantile(pilots[, "d_hat"], .025), quantile(pilots[, "d_hat"], .975)))
cat(sprintf("Pilots that lead to a main study (planned n <= %d): %.0f%%\n", max_n, 100 * mean(go)))
cat(sprintf("Of those main studies: median true power %.2f; %.0f%% have power below .80\n",
            median(power_naive), 100 * mean(power_naive < 0.80)))
cat(sprintf("Promising pilots overestimate: mean pilot d among go-ahead studies = %.2f\n",
            mean(pilots[go, "d_hat"])))

# Alternative 1: power on the smallest effect size of interest (SESOI), chosen
# from what would matter in practice, e.g. d = 0.35.
cat(sprintf("Powering on a SESOI of d = 0.35: n = %d per group\n", planned_n(0.35)))

# Alternative 2: safeguard power, using the lower 60%% CI limit from the pilot.
n_safe <- sapply(pilots[, "d_safe"], planned_n)
go_safe <- n_safe <= max_n
cat(sprintf("Safeguard plan: %.0f%% of pilots fit the budget; median true power %.2f among them\n",
            100 * mean(go_safe), median(sapply(n_safe[go_safe], true_power))))