This week I havent had much time for R stuff either, but my honours students have been writing preregistration documents and doing power analyses. Usually I point students toward the gpower calculator but for some reason my student couldn’t make it work so we ended up trying out the pwr package and it is awesome.
the pwr package
The pwr package (vignette here) makes it easy to do a power analysis to work out how many participants you will need to test, given the design of your experiment and statistical tests you have planned.
library(pwr)
library(tidyverse)
library(ggannotate)cohens D
There are different functions depending on what statistical analysis you are planning. In case you can’t remember the cohen’s effect size levels, there is a cohen.ES() function that will remind you.
cohen.ES(test = "t", size = "medium")##
## Conventional effect size from Cohen (1982)
##
## test = t
## size = medium
## effect.size = 0.5calculate sample size
When talking about power, we generally need information about the test we want to do, the effect size we are trying to detect, the significance level we are going to use and the sample size. Most often we are looking to “solve” for sample size (n), that is work out how many participants we need to test to detect an effect of X size with Y power. Other times (less commonly) we have run the study and want to determine what power we have to detect an medium effect given the sample we have tested, for example.
There are different functions for common stats tests (see below), for each, you can just leave out whichever piece of information you don’t know and it will solve for that.
In R, when you type a function in a chunk (or the console), put your cursor within the bracket and press TAB, it will show you which arguments the function takes.
independent samples t-test
Lets imagine we are trying to work out how many participants we need to test in each group if we are going to use an independent samples t-test. We don’t know n, so we leave that out.
We know that we are looking for a medium effect; you don’t even need to remember that for t-tests a medium effect = 0.5; the function take a character string “medium” to the test arguments directly. We always go with a significance level of 0.05 and 80% power. The type of t.test refers to one, two, or paired; here we are using independent samples so it is a type = “two”. Alternative refers to whether we want one or two sided tests; the default is two sided but we will spell it out here.
pwr.t.test(d = 0.5, sig.level = 0.05, power = 0.8, type = "two", alternative = "two.sided")##
## Two-sample t test power calculation
##
## n = 63.76561
## d = 0.5
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: n is number in *each* groupOK to detect a medium effect size (d = 0.5) with 80% power if we are testing different groups of participants, we need 64 participants PER GROUP.
paired samples t-test
Lets imagine we are going to run the same experiment with a within subjects design. Each participant will do condition 1 and 2 in a counterbalanced order and we will use a paired-sample t-test to test whether there is a difference in their performance by condition.
pwr.t.test(d = 0.5, sig.level = 0.05, power = 0.8, type = "paired", alternative = "two.sided")##
## Paired t test power calculation
##
## n = 33.36713
## d = 0.5
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: n is number of *pairs*ooooo NICE! with a within subjects design to detect the same magnitude of effect we only need 34 pairs, that is N = 34 total.
plot power by sample size
plot()
There is a plot() function that you can also visualise the relation between power and sample size.
# define a range of possible n values, from 10 to 100 in 10s
n <- seq(10,100,10)
# get the power to detect a medium effect size for each of those ns
p.out <- pwr.t.test(d = 0.5,
n = n,
sig.level = 0.05)
# plot
plot(p.out)## Warning in if (x$power < 0.5) {: the condition has length > 1 and only the first
## element will be used
Hmmm not sure how to move those captions so that the points at the top are visible. Maybe it would be better to turn the n and power values into a dataframe and use ggplot?
ggplot()
# make a dataframe of n by power
n_p <- data.frame(n, power = p.out$power)
n_p## n power
## 1 10 0.1850957
## 2 20 0.3379390
## 3 30 0.4778965
## 4 40 0.5981469
## 5 50 0.6968934
## 6 60 0.7752659
## 7 70 0.8358223
## 8 80 0.8816025
## 9 90 0.9155872
## 10 100 0.9404272Plot n by power… adding annotations with the ggannotate package.
n_p %>%
ggplot(aes(x = n, y = power)) +
geom_point() +
geom_line() +
geom_vline(xintercept = 64, linetype="dashed") +
geom_hline(yintercept = 0.8, linetype="dashed") +
geom_text(data = data.frame(x = 75.1352579178505, y = 0.68943731144711, label = "N = 64 per group \n power = 0.8"),
mapping = aes(x = x, y = y, label = label),
inherit.aes = FALSE) +
labs(title = "Two sample t.test, power to detect medium effect")
As always, ggplot is always easier to customise!!
The pwr package seems super easy to use and contains functions for all the kinds of tests we typically do in honours. Check it out!
