# Sample size for a four-group dose study (vehicle and three doses), from the planned analysis.
# Expected means 100, 95, 90, 85 (vehicle, low, mid, high); SD 12 within groups; two-sided alpha 0.05; power 80%.
# Power from the noncentral t and F distributions. R 4.3.3, base R.
mu <- c(100, 95, 90, 85); sd <- 12; k <- length(mu); alpha <- 0.05; target <- 0.80
power_t <- function(diff, n, a) {
  df <- k * n - k; ncp <- diff / (sd * sqrt(2 / n)); tc <- qt(1 - a / 2, df)
  pt(tc, df, ncp, lower.tail = FALSE) + pt(-tc, df, ncp)
}
fval <- sqrt(mean((mu - mean(mu))^2)) / sd
power_f <- function(n) { df1 <- k - 1; df2 <- k * n - k; pf(qf(1 - alpha, df1, df2), df1, df2, fval^2 * k * n, lower.tail = FALSE) }
power_trend <- function(n, cc = c(-3, -1, 1, 3)) {
  df <- k * n - k; ncp <- sum(cc * mu) / (sd * sqrt(sum(cc^2) / n)); tc <- qt(1 - alpha / 2, df)
  pt(tc, df, ncp, lower.tail = FALSE) + pt(-tc, df, ncp)
}
smallest_n <- function(fun) { for (n in 2:499) if (fun(n) >= target) return(n) }
cat(sprintf("Cohen's f from the expected means: %.3f\n", fval))
n <- smallest_n(power_f); cat(sprintf("Omnibus ANOVA F-test: n = %d per group (power %.3f)\n", n, power_f(n)))
n <- smallest_n(power_trend); cat(sprintf("Linear trend contrast: n = %d per group (power %.3f)\n", n, power_trend(n)))
for (diff in c(15, 10, 5)) for (i in 1:2) {
  a <- c(alpha, alpha / 3)[i]; lab <- c("alpha 0.05", "alpha 0.05/3 (Bonferroni for 3 comparisons)")[i]
  n <- smallest_n(function(n) power_t(diff, n, a))
  cat(sprintf("Dose vs vehicle, difference %d, %s: n = %d per group (power %.3f)\n", diff, lab, n, power_t(diff, n, a)))
}
ns <- Filter(function(n) k * n - k >= 10 && k * n - k <= 20, 2:49)
cat(sprintf("Resource equation (error df 10 to 20): n = %d to %d per group\n", min(ns), max(ns)))
for (n in ns) cat(sprintf("  n = %d: power for difference 10 at alpha 0.05/3 = %.3f; for difference 15 = %.3f\n", n, power_t(10, n, alpha / 3), power_t(15, n, alpha / 3)))