Power analysis through simulation in R
Niklas Johannes
Age predicts grumpiness with a large effect. But the sample is too small for significance.
set.seed(1)
age <- runif(10, 20, 80)
grumpiness <- 50 + 0.5 * age + rnorm(10, 0, 20)
cor.test(age, grumpiness)
Pearson's product-moment correlation
data: age and grumpiness
t = 1.5624, df = 8, p-value = 0.1568
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2100550 0.8533538
sample estimates:
cor
0.4835198 Age predicts grumpiness with a super tiny effect, but we have a sample of a million, so the effect is significant.
age <- runif(1e6, 20, 80)
grumpiness <- 50 + 0.01 * age + rnorm(1e6, 0, 20)
cor.test(age, grumpiness)
Pearson's product-moment correlation
data: age and grumpiness
t = 8.5288, df = 999998, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.006568625 0.010488268
sample estimates:
cor
0.008528479 Age doesn’t predict grumpiness. Can a nonsignificant p-value tell us that?
Pearson's product-moment correlation
data: age and grumpiness
t = 0.093827, df = 999998, p-value = 0.9252
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.001866138 0.002053791
sample estimates:
cor
9.382711e-05 Remember \(P(data|H)\), not \(P(H|data)\)?
“…If one in twenty does not seem high enough odds, we may, if we prefer it, draw the line at one in fifty or one in a hundred. Personally, the writer prefers to set a low standard of significance at the 5 per cent point, and ignore entirely all results which fails to reach this level. A scientific fact should be regarded as experimentally established only if a properly designed experiment rarely fails to give this level of significance.”
We have three independent groups: control, treatment A, and treatment B. The pesky ethics board asks us to do a power analysis. You head to GPower.
Thankfully, there’s a previous study! It had n = 20 per condition and the conditions are only somewhat similar to our planned experiment, but they do report an effect size: \(\eta_2 = .21\). Off to GPower!
set.seed(42)
d <- data.frame(
id = 1:60,
condition = rep(c("control", "Treatment A", "Treatment B"), times = 20),
score = rnorm(60, mean = c(0, 10, 20), sd = 15)
)
model <-
aov(
score ~ condition, data = d
)
effectsize::eta_squared(model)# Effect Size for ANOVA
Parameter | Eta2 | 95% CI
-------------------------------
condition | 0.21 | [0.06, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
Cohen (1988)
\(d = \frac{M_1-M_2}{pooled\ \sigma}\)
\(pooled\ \sigma = \sqrt{\frac{(sd_1^2 + sd_2^2)}{2}}\)
Control group has a mean of 100 and an SD of 20. The treatment group has a mean of 105 and an SD of 10. The difference in the means is \(105-100 = 5\) (simplified). The pooled SD is (simplified!) \(\frac{20+10}{2} = 15\). So our difference is \(5/15\) or simply \(d = 0.33\). In other words, our difference is a third of a standard deviation unit.
Cohen suggested (and later very much regretted) some rules of thumb if a researcher has no better idea:
In small samples, Cohen’s d will be biased. Use Hedge’s g instead. In fact, you should probably always use g. (Software does it for you anyway.)
\(d = \frac{M_1-M_2}{pooled\ \sigma^*}\)
\(r = B_{xy} \frac{\sigma_x}{\sigma_y}\)
We predict grumpiness with age. The regression slope is 2: With each year, people score 2 higher on grumpiness. The SD of grumpiness is 30. The SD of age is 10. The correlation coefficient is \(2*10/30 = .67\). We could’ve also just standardized both variables and run a regression.
Cohen also provides a formula how to get \(r\) from \(d\). Remember, use Hedge’s \(g\) instead of \(d\).
\(r = \frac{d}{\sqrt{d^2 + 4}}\)
Back to that medium effect size:
\(r = \frac{0.5}{\sqrt{0.5^2 + 4}} = 0.24\)
Strength of association is just another way of saying magnitude of shared variance between variables. Or: Does the blue line do better than the black line?
\(Variance \ explained = \frac{\sigma_{effect}}{\sigma_{total}}\)
In the ANOVA context, we often use \(\eta^2\), because it has been standard in SPSS output (Lakens 2013).
\(\eta^2 = \frac{SS_{effect}}{SS_{total}}\)
\(\eta^2_p = \frac{SS_{effect}}{SS_{total} + SS_{error}}\)
\(f\) mostly used for one-way ANOVAs
\(f^2\) mostly used for regressions, but also one-way, or multi-way ANOVAs
These effect sizes of shared variance are often biased. Instead, use \(\omega^2\) or \(\epsilon^2\). Don’t panic: Smart people have provided spreadsheets.
Effect size converter: https://osf.io/vbdah/
All you need to remember:
Squaring the r is not merely uninformative; for purposes of evaluating effect size, the practice is actively misleading. (Funder and Ozer 2019, 3)
The moment we move beyond two groups or bivariate relationships:
As a rule, reports of effect size should focus on 1 df effects. (Baguley 2009, 614)
So… is \(r\) = .21 big then? (Meyer et al. 2001)
When we correlate variables that are specifically selected not to be related, we still reach \(r\) ~ .10.
Several considerations (Funder and Ozer 2019):
Smallest effect size of interest (Anvari et al. 2021)
Minimally detectable difference
Often used in medicine:
Global ratings of change methods:
Procedure (Anvari and Lakens 2021):
My study has 80% power to detect a medium sized effect, as shown by the meta-analysis by XYZ.
Translation: If this doesn’t work, we have learned close to nothing.
I designed my study to be able to detect an effect of a certain size with 95% power. Anything smaller than that is uninteresting. Don’t waste resources if you’re hoping to find an effect this large.
Translation: I thought about what I want and I’m putting that part of the process up for debate.
Unstandardized measures have several advantages:
