Motivating Example: You’ve decided a paired design is the right experimental setup! You know that pairing can add statistical power, but you want to know how to do that best. You also want to know the sample size you need to come to a reasonable scientific conclusion.
Learning Goals: By the end of this section, you will be able to:
pwr package to find the power of a planned paired t-test.presize package to find the expected precision (confidence interval width) of a planned paired t-test.I will not go over design for a one-sample t-test, because these come up very rarely – and besides, the power analysis logic below still holds.
Experimental considerations for paired t-tests. Before we end this chapter, we will consider designing an experiment with a “paired” design. Remember that the benefit of a paired design is that we can increase power and precision by minimizing extraneous variation unrelated to the treatment. But in doing this, we should ensure that no bias arises, and that the pairs are adequately matched such that we decrease extraneous variability (otherwise we do not gain from a paired approach). Then, we can consider the sample size needed for the power and precision we desire.
Like all experiments, a paired design is susceptible to bias. And the ways to guard against this bias are just specific instances of the best practices in experimental design: random assignment, and avoiding experimental artifacts.
Paired designs can provide us with more power than unpaired designs. But this is only true if pairing actually reduces extraneous variation unrelated to treatment. If pairs are made arbitrarily, a paired design could have less power than an unpaired one.
Pairing units after the fact, based on a measured covariate, can be a legit technique. But be careful. Pairing based on the response variable itself, after you’ve already seen it, is not ok. This is wrong, and makes you a bad statistician and an even worse person.
Recall that power quantifies the probability that we will correctly reject a false null hypothesis. Precision, by contrast, describes how certain we can be about our estimate – quantified, for example, by the width of its confidence interval. First I will show you how to quickly find both with R packages. Then I will give you more intuition by simulating.
A weird thing about power and precision analyses is that you need to specify an effect size before you do your study. When I was a student I found this super confusing. Like, if I already know the effect, why am I doing the experiment? The effect size we specify in these analyses isn’t the one that actually exists in the world – we don’t know that yet. It’s the smallest effect we’d want to be able to detect, if it’s really there.
How big of a sample do you need? Well, of course, the answer depends on what you need the sample for.
Say you wanted to detect a small effect, because even a little bit of harm from something nearly everyone is exposed to (e.g. food dye, microplastics etc.) would add up to a real public health cost. A Cohen’s d of 0.2 sits just over the hump of what’s conventionally called a “small” effect. So, how big a sample would you need to reject the null, and estimate this effect relatively precisely? The tools below can help!
pwr package for power analysesThe pwr package can tell us the power of a paired t-test directly. You provide the sample size (number of pairs), the effect size (Cohen’s d = \(\frac{\text{mean difference between pairs}}{\text{sd of the difference between pairs}}\)) you want to be able to detect, and your significance threshold, \(\alpha\) (by tradition \(\alpha = 0.05\)), and it finds your power.
library(pwr)
pwr.t.test(d = 0.2,
n = 60,
sig.level = 0.05,
type = "paired",
alternative = "two.sided")
Paired t test power calculation
n = 60
d = 0.2
sig.level = 0.05
power = 0.3316786
alternative = two.sided
NOTE: n is number of *pairs*So, we will reject the null in about one third of paired experiments, if we’re planning around a Cohen’s d of 0.2. This is not enough power for me, and we would actually substantially overestimate the effect when we did reject the null. I would bump my sample size up to 270 to reach ninety percent power.
presize package to plan for precisionThe presize package specializes in finding the precision for a planned experiment. While it doesn’t explicitly have a “paired” option, we know that a paired t-test is just a one-sample t-test for the difference between pairs, with a null of zero. So we can use prec_mean() to find the expected width of the confidence interval around that mean difference. Recall that here (and above) n is the number of pairs, not the number of individuals.
Note that the effect size is playing a slightly different role here than it did in our power calculation. There, d = 0.2 was the smallest effect we wanted to be able to detect. Here, mean = 0.2 is the effect size we’re hypothetically planning around – if the true difference really is about this size, how tightly would this sample size pin it down?
So, for example, we can find the expected 95% CI for our estimate of the mean difference, planning around a mean of \(\mu = 0.2\), \(\sigma = 1\), and a sample of 270 pairs (our larger study that gave 90% power), as given by:
library(presize)
prec_mean(
mean = 0.2,
sd = 1,
n = 270,
conf.level = 0.95
)# A tibble: 1 × 7
mean sd n conf.width conf.level lwr upr
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0.2 1 270 0.240 0.95 0.0802 0.320The output tells us our 95% CI would run from about 0.08 to 0.32 – a width of 0.24. In other words, even with a well-powered sample of 270 pairs, we’d still expect our estimate of the true mean difference to be uncertain by roughly ±0.12. Planning for power tells you whether you’re likely to detect an effect at all; planning for precision tells you how well you’ll actually pin it down once you do – two different, complementary questions worth asking about the same design.
Now try this yourself. Use the web-R environment below (which has presize loaded) to compare the precision of our underpowered study of 60 pairs, using the same Cohen’s d as above.
With 270 pairs, the width of the 95% CI is 0.32 − 0.08 = 0.24. Now find the width of the 95% CI from our original design (n = 60).
Q1) What is the width of the 95% CI at n = 60?
Width = upper bound − lower bound, same calculation as the n = 270 example above, just using your own n = 60 output this time.
Q2) A sample size of 270 is 4.5 times a sample size of 60. The 95% CI from the n = 60 design is about times as wide as the CI from the n = 270 design.
Notice that a four-and-a-half-times increase in sample size only shrunk the CI width by a bit more than two times. This is what we expect – CI width doesn’t shrink in direct proportion to sample size, it shrinks with \(\sqrt{n}\) (note \(\sqrt{4.5} = 2.1\)). This reflects our rule of thumb that \(n \propto \sigma^2/\text{effect}^2\) (see the section on Power and Precision in Chapter 12). This is why it’s so much work and money to precisely estimate small effects.
For many questions, small effects don’t matter much in practice, so this cost isn’t worth paying. But for some topics, even a small effect can have huge implications – worth knowing which kind of question you’re actually asking before you commit to a sample size.
A webapp for power: I found the n = 270 (the actual answer is 263, but 270 is a bit easier to write) above from the webapp: Inference for a Mean: Comparing a Mean to a Known Value – this is a one-sample t-test, which is equivalent to a paired t-test with a \(\mu_0 = 0\). So, I set \(\mu_0\) (the null difference) to zero, and \(\mu_1\) (the difference we hope to be able to find) to 0.2 and \(\sigma = 1\). I kept \(\alpha\) at its traditional 0.05, and bumped desired power up to 0.90.
A webapp for precision: The same presize package we used above is also available as a webapp here.
For me, the math never quite hits right. I find it much more satisfying to “see” chance via computer simulation.
To simulate for power we:
Make a sample of n pairs.
rnorm() (note: This isn’t an assumption of the paired t-test, but it’s easy).Run a statistical test and summarize the data.
Do this many times.
See our power and precision from this simulation.
Steps 1 - 3 are implemented in the functions I wrote, below:
library(dplyr)
library(broom)
# ONE SIMULATION
simulatePairs <- function(n_pairs, diff, sd){
# simulate paired data
pairs <- tibble(control = rnorm(n_pairs)) |>
mutate(treatment = control + rnorm(n_pairs, mean = diff, sd = sd))
# paired t-test
test_result <- with(pairs,
t.test(treatment, control, paired = TRUE))
# return the quantities we're interested in
glance(test_result) |>
mutate(cohens_d_estimate = statistic / sqrt(parameter + 1)) |>
dplyr::select(diff_estimate = estimate, cohens_d_estimate, p.value)
}
# MANY REPS (1000 is default)
pairedPower <- function(n_pairs, diff, sd, n_reps = 1000){
many_sims <- replicate( n = n_reps,
expr = simulatePairs(n_pairs, diff, sd),
simplify = FALSE) |>
bind_rows() |>
mutate(diff_param = diff,
replicate = 1:n(),
n_pairs = n_pairs,
cohens_d_param = diff / sd)
return(many_sims)
}Now, we can run the simulation:
power_sim <- pairedPower(n_pairs = 270, diff = 0.2, sd = 1)# A tibble: 1,000 × 7
diff_estimate cohens_d_estimate p.value diff_param replicate n_pairs
<dbl> <dbl> <dbl> <dbl> <int> <dbl>
1 0.217 0.231 0.000186 0.2 1 270
2 0.173 0.170 0.00570 0.2 2 270
3 0.159 0.162 0.00827 0.2 3 270
4 0.113 0.111 0.0686 0.2 4 270
5 0.245 0.256 0.0000348 0.2 5 270
6 0.275 0.283 0.00000505 0.2 6 270
7 0.177 0.194 0.00161 0.2 7 270
8 0.212 0.218 0.000397 0.2 8 270
9 0.190 0.189 0.00209 0.2 9 270
10 0.197 0.192 0.00178 0.2 10 270
# ℹ 990 more rows
# ℹ 1 more variable: cohens_d_param <dbl>We can now, e.g. find power (at \(\alpha = 0.05\)) as the proportion of simulations with p < 0.05:
power_sim |>
mutate(significant = p.value < 0.05)|>
summarise(power = mean(significant))# A tibble: 1 × 1
power
<dbl>
1 0.905This simulation verifies our (ok the R packages’ and webapps’) math – we have 90% power to reject the null of no difference between pairs, when we have 270 pairs and a true mean difference of Cohen’s d = 0.2!