Power Calculations Overview
Statistical power is the probability that,
- If there is an effect
- That your model will detect that effect
- At a given confidence level
Typically, power will depend on:
- The size of the effect (the larger the effect, the easier it is to gain statistical power)
- The sample size (the larger the sample size, the easier it is to gain statistical power)
- The amount of “noise” (the more noise, the harder it is to gain statistical power)
Let’s take an easy example. Let’s say that you’re estimating how much height impacts basketball ability. You take the statistics from your local rec league, and ask the administrator to give you statistics for a sample of players along with their heights. The question is how many players you will need to detect an effect.
set.seed(1)
N <- 100
effect <- 7 # for every foot of height, we expect 7 additional points per game
height <- rnorm(N, 5.8, 0.3) ## player heights
ppg <- pmax(3 + effect * (height-mean(height)) + rnorm(N, 0, 3), 0) # we truncate ppg at zero because it is bounded below
summary(lm(ppg ~ 1 + height))
The result of this call is:
Call:
lm(formula = ppg ~ 1 + height)
Residuals:
Min 1Q Median 3Q Max
-4.4927 -1.7958 -0.6446 1.5550 6.8470
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -27.6960 5.5725 -4.970 2.84e-06 ***
height 5.3079 0.9544 5.562 2.32e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.559 on 98 degrees of freedom
Multiple R-squared: 0.2399, Adjusted R-squared: 0.2321
F-statistic: 30.93 on 1 and 98 DF, p-value: 2.317e-07
We estimated an effect of height at 5.3079 and our confidence level that this coefficient is above zero is one minus 2.32e-07 (the p-value of the coefficient), which is effectively 100%.
In this one simulation, we picked up an effect. But power is the probability that we will pick up this effect, so in order to understand that, we need to repeatedly run this experiment:
N <- 100
set.seed(1)
confidence_level = 0.95
effect <- 7 # for every foot of height, we expect 7 additional points per game
is_significant <- c()
N_replicates <- 100
for(i in 1:N_replicates) {
height <- rnorm(N, 5.8, 0.3) ## player heights
ppg <- pmax(3 + effect * (height-mean(height)) + rnorm(N, 0, 3), 0) # we truncate ppg at zero because it is bounded below
is_significant[i] = (1-summary(lm(ppg ~ 1 + height))$coef[2, 4]) > confidence_level
}
mean(is_significant)
The mean of is_significant is 1, meaning that at this sample size, confidence level, and effect size, we will detect an effect with confidence ~100% of the time.
Let’s change the parameters, specifically to reduce the sample size to just 20.
N <- 20
set.seed(1)
confidence_level = 0.95
effect <- 7 # for every foot of height, we expect 7 additional points per game
is_significant <- c()
N_replicates <- 100
for(i in 1:N_replicates) {
height <- rnorm(N, 5.8, 0.3) ## player heights
ppg <- pmax(3 + effect * (height-mean(height)) + rnorm(N, 0, 3), 0) # we truncate ppg at zero because it is bounded below
is_significant[i] = (1-summary(lm(ppg ~ 1 + height))$coef[2, 4]) > confidence_level
}
mean(is_significant)
Now, the mean of is_significant is just 0.83. Meaning that with a sample size of 20, we will only detect this effect with confidence ~83% of the time.
Additionally, is_significant only has 100 trials. It itself is a statistic, with associated uncertainty. Here, we can use the uncertainty of proportions:
phat +/- z * sqrt(phat * (1-phat) / n)
Where phat is the observed proportion, z is the critical statistic (1.96 for 95%) and n is the sample size:
phat = mean(is_significant)
phat + c(-1, 1) * 1.96 * (sqrt(phat * (1-phat)/length(is_significant)))
Which gives us a confidence interval of 0.756376 0.903624 for the power of this experiment.
Let’s say that we’re doing a difference-in-difference (DiD) experiment. In a DiD, we separate our population into two groups, the control and treatment, and we survey both groups before and after we do the treatment.
The model for this is:
lm(outcome ~ treated + wave + control:treatment)
Where
treatedrepresents the differences between the treatment and control groups before the experimentwaverepresents the differences between the groups before and after the treatmenttreated:waverepresents the effect of the treatment
With simple random sampling, this would be a straightforward simulation that is only marginally different from the one above. But with the Rhetoric project, we are not sampling randomly; we are sampling from ZIP codes, and it is the ZIP codes that have been allocated to the treatment and control groups. To properly assess the power of the experiment, need to take this into account.
Here, the parameters that will matter are:
- How many ZIP codes we sample from
- How many respondents we get in each zip code
- How different the ZIP codes are from each other
- How large of an effect the treatment has
- How much noise there is in individual responses
We are doing this experiment in California, so we need some base statistics:
N_ZIPS_IN_CALI = 1741
POP_CALI = 39e6
And since we are developing a budget based on a total number of expected responses, we set that here:
TOTAL_SAMPLE_PER_WAVE = 4000
Since we are sampling from within each ZIP code, and each ZIP code has a varying number of people in it, we have to model how many people come out of each ZIP code explicitly:
First, we develop model the heterogeneity in the population across ZIP Codes (we could look up the population statistics, but this will work for now)
pop_in_zips = rnbinom(N_ZIPS, mu = POP_CALI/N_ZIPS_IN_CALI, size = 1)
And here we develop an expected sample per zip code, per wave:
expected_sample = TOTAL_SAMPLE_PER_WAVE * pop_in_zips/sum(pop_in_zips)
Finally, we can model what our sample actually ends up being using a poisson (count) distribution with a mean of our expected sample in each zip code:
sample_from_zips_1 = unlist(map(expected_sample, ~ rpois(1, .x)))
sample_from_zips_2 = unlist(map(expected_sample, ~ rpois(1, .x)))
Now that we’ve modeled the sampling process, we now need to model our observations.
Here, we set up the statistics that will control our observations
variation_between_zips = .5 # how much variation is there between zip codes?
effect_size = 0.03 # how much effect does our treatemnt have
noise_to_signal = 1 # how much noise is there relative to our treatment when respondents answer?
common_trend = 0.03 # how much of a temporal effect was there between the first and second wave?
zip_means = rnorm(N_ZIPS, 0, variation_between_zips)
zip_index_1 = map(1:N_ZIPS, function(i) {
if(sample_from_zips_1[i] == 0) return(NULL)
data.frame(
zip = i,
observation = rep(1, sample_from_zips_1[i])
)
}) %>%
reduce(bind_rows) %>%
mutate(observation = cumsum(observation))
zip_index_2 = map(1:N_ZIPS, function(i) {
if(sample_from_zips_2[i] == 0) return(NULL)
data.frame(
zip = i,
observation = rep(1, sample_from_zips_2[i])
)
}) %>%
reduce(bind_rows) %>%
mutate(observation = cumsum(observation))
control_zips = which(cumsum(sample_from_zips_1) <= TOTAL_SAMPLE_PER_WAVE/2)
treated_zips = setdiff(1:N_ZIPS, control_zips)
treated_zips_idx = rep(0, N_ZIPS); treated_zips_idx[treated_zips] <- 1
observations_1 = zip_index_1 %>%
mutate(wave = 0) %>%
mutate(treated = treated_zips_idx[zip]) %>%
mutate(observation = zip_means[zip] + rnorm(nrow(.), 0, noise_to_signal*effect_size))
observations_2 = zip_index_1 %>%
mutate(wave = 1) %>%
mutate(treated = treated_zips_idx[zip]) %>%
mutate(observation = zip_means[zip] + effect_size + rnorm(nrow(.), 0, noise_to_signal*effect_size))
observations = bind_rows(observations_1, observations_2)
Now, the data looks like this:
zip observation wave treated
1 55 0.13463329 1 0
2 62 0.20774712 0 1
3 32 1.10128628 0 0
4 62 0.25152150 1 1
5 92 -0.05326673 1 1
6 17 0.40706148 1 0
7 98 0.91711999 1 1
8 98 0.86775583 0 1
9 26 -0.31065658 0 0
10 69 -0.57515213 0 1
Finally, our statistical model will look like this:
mod = lmer(observation ~ (1 | zip) + wave + treated, data = observations)
summary(mod)$coefficients[3, ]
On this particular run, I get the following:
Estimate Std. Error t value
0.1263883 0.1178889 1.0720968
And since the t value is less than 1.96 (the critical value for 95% confidence), on this run, we did not detect the effect (even though it was present). To estimate the power for this experiment, you would need to run the model several (hundreds of) times and take the proportion of times that your results were significant.
❗ While this model captures the sampling process, it likely does not capture the observational model, as we almost certainly would be using a likert-scale response to capture survey data. In order to model that, you would need to first generate the data according to an ordinal model, and then to use an ordinal model to evaluate the simulated data.