I am interested in the legitimacy and classification of the following approach to sample size determination. This question is inspired by “There is only one test!” by Allen Downey and “Statistics without the agonizing pain” by John Rauser.
Suppose there is a customer base, and each customer is described by some characteristic $x$, such as the total amount of spendings. Let $f$ be an arbitrary metric that we care about, such as the average value. Let also $T$ be a treatment that is hypothesized to improve the business as measured by $f$. We are then to run the following NHST:
$H_0$: $T$ does not affect the population as measured by $f$.
$H_1$: $T$ positively affects the population as measured by $f$.
Given an $\alpha$ and a minimum detectable effect size $\delta$, the goal is to determine the minimum sample size $n$ required to attain a prescribed $\beta$. At our disposal, we have a dataset $X = \{x_i\}$ describing the customer base during a time interval in the recent past. We resort to simulation as follows.
First, we assume some $n$ and draw $n$ samples with replacement from $X$ and apply $f$. We repeat it multiple times. By doing so, we obtain an empirical estimate of the sampling distribution under the null hypothesis. We then assume that, under the alternative hypothesis, the distribution will have the same shape but will be shifted by $\delta$. We repeat the sampling with $\delta$ as an offset and obtain an estimate the sampling distribution under the alternative hypothesis. We can then compute $\beta$.
Placing the above procedure in an optimization loop, we obtain the minimum $n$ required to attain the desired $\beta$. The following R snippet demonstrates the approach:
library(tidyverse)
set.seed(42)
data <- tibble(value = rlnorm(20000))
metric <- mean
alpha <- 0.05
beta <- 0.2
delta <- 0.05 * metric(data$value)
simulate <- function(n, N) {
control <- replicate(N, { metric(sample(data$value, n, replace = TRUE)) })
treatment <- replicate(N, { metric(sample(data$value, n, replace = TRUE)) }) + delta
critical <- quantile(control, 1 - alpha)
beta <- sum(treatment < critical) / N
list(control = control,
treatment = treatment,
critical = critical,
beta = beta)
}
N <- 1000
target <- function(n) beta - simulate(as.integer(n), N)$beta
n <- as.integer(uniroot(target, interval = c(1, 10000))$root)
result <- simulate(n, N)
tibble(control = result$control, treatment = result$treatment) %>%
gather(group, value) %>%
ggplot(aes(value, fill = group)) +
geom_density(alpha = 0.5, color = NA) +
geom_vline(xintercept = result$critical, linetype = 'dashed')
Is it a statistically sound procedure? How would you classify it? Does it have a name? Is it a variant of bootstrap? How does it compare with the classical t-test?
The last is particularly concerning, as it gives a totally different estimate:
power.t.test(delta = delta,
sd = sd(data$value),
sig.level = alpha,
alternative = 'one.sided',
power = 1 - beta)
The former gives 4822 per group, while the latter 8620 per group. I have also tried with rnorm(20000, mean = 100, sd = 50)
instead of rlnorm(20000)
and got 634 and 1261 per group, respectively. What accounts for this difference?
Aftermath
For those interested, I have written a blog post describing the technique in more detail.
power.t.test
interface. (These appear to be the two most common reasons people post this kind of question here). Search our site for power simulation for many examples. $\endgroup$