# Sample size determination using historical data and simulation

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. • This seems like an account of the concept of power as applied to a location-shift test. If you are getting "totally different" sample size estimates then likely there's a bug in your code or you are misunderstanding the 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. – whuber Sep 13 '19 at 20:13 • Then I am also asking for help with understanding my potential mistake. – Ivan Sep 15 '19 at 13:32 ## 2 Answers You're doing what's called a "bootstrap" power calculation. It's a perfectly reasonable method. Though not as common as methods that make assumptions about the distribution of the test statistic, it is valid, and likely to be more appropriate in cases where the necessary assumptions are not true enough. That is, the assumptions are almost never exactly true, but usually they're plenty good enough. I didn't find any references about this straightforward situation in a quick search of the literature, but here's a blog post from statistics you can probably trust that refers to Common Errors in Statistics (and How to Avoid Them) by Phillip Good and James Hardin, which I've not personally read but have heard good things about, and both authors are respected statisticians. Anyway, the problem with your code is that it's doing a one-sample t-test, not a two-sample t-test. power.t.test for the one-sided situation gives about 631, close to your 634. Try comparing with this. power.t.test(delta=delta, sig.level=alpha, power=1-beta, sd=sd(data$value),
alternative='one.sided', type="one.sample")


That is, the test statistic you're using is the mean of the test group, and comparing it with a fixed known control. This is a different situation than if you're using the difference between two random groups as the test statistic.

Here's a stab at what equivalent code for a two-sample test would look like; I get 1306, close to the expected 1261.

simulate <- function(n, N) {
control_under_null <- replicate(N, { metric(sample(data$$value, n, replace = TRUE)) }) treatment_under_null <- replicate(N, { metric(sample(datavalue, n, replace = TRUE)) }) difference_under_null <- treatment_under_null - control_under_null control <- replicate(N, { metric(sample(datavalue, n, replace = TRUE)) }) treatment <- replicate(N, { metric(sample(data$$value, n, replace = TRUE)) }) + delta
difference <- treatment - control
critical <- quantile(difference_under_null, 1- alpha)
beta <- sum(difference < critical) / N
list(null = difference_under_null,
alternative = difference,
critical = critical,
beta = beta)
}


PS. The Downey blog is misleading about the current state of statistical education; teaching via simulation and bootstrap is now very common, and many consider it best practices. For an example see one of the OpenIntro statistics textbooks: https://www.openintro.org/stat/textbook.php?stat_book=isrs

• Indeed! Thank you for investigating this! I think the problem with the code can now be considered solved. Zooming out a bit, how would you classify this approach to sample size determination, and is it a sound strategy? Could you please state this at the beginning of your answer? I will then accept it. – Ivan Sep 19 '19 at 11:56

What you have coded is essentially a bootstrap test, as it is sampling from existing data with replacement.

The reason you get different answers here is that power.t.test is assuming your input data come from a linear model (i.e. a Normal distribution, after accounting for any covariates). Your actual data come from a lognormal distribution, which has much fatter tails and so needs more sample to detect a given effect.

Try repeating your code but with your initial sample from rnorm instead of rlnorm. You should find the results are approximately the same.

In real life as well, it would be fairly unusual to have an additive treatment effect on lognormally distributed data (multiplicative would be much more common)

• I have tried with rnorm(20000, mean = 100, sd = 50) and got 634 and 1261 per group, respectively. It is interesting that they remain different by a factor of 2. – Ivan Sep 17 '19 at 5:13
• Regarding bootstrap, can it be classified as such if the number of samples drawn is different from the size of the original sample? According to what I have seen, the sizes usually match. Another concern is that here it is about sample size determination, not the actual testing. Plus there is this location-shift assumption mentioned by whuber. – Ivan Sep 17 '19 at 5:20
• What whuber calls a location-shift assumption is what I am calling an additive treatment effect, same thing. A bootstrap does assume the samples are the same size, but that's the bootstrap sample size and the size of the putative future experiment, not the reference data you are using here. – JDL Sep 17 '19 at 7:57
• I still cannot find the reason the two approaches give different results, even for normal data. Regarding “needs more sample to detect a given effect,” should it not affect both methods? When bootstrapping, it should also lead to a larger sample size. Or am I missing something? – Ivan Sep 18 '19 at 5:41
• "should it not affect both methods?" — no, because power.t.test is assuming the tails are normally-shaped tails, whereas the bootstrap method is not making that assumption. – JDL Sep 18 '19 at 7:53