Power analysis

Planning sample sizes by simulation

Dr Massimiliano Canzi

Power analysis

Did you know: the peacock mantis shrimp smashes snail shells with a club that reaches 14 to 23 metres per second. The energy for the strike is stored in a saddle-shaped spring in its limb.

Source: Patek, Korff & Caldwell (2004), Nature

Where we are

Every data set in this course so far was simulated. We chose the true effects, generated the data, and checked whether the model found them again.

Today we turn that around. Before a study is run, the same kind of simulation can tell us how often a planned design would detect an effect of a given size: its power. From there we can work out how many participants and items we need.

Today

  1. What power is
  2. Power by simulation
  3. Power with {simr}
  4. Why post hoc power tells us nothing

What power is

Did you know: a mantis shrimp strike hits twice. The club strikes first, and about 0.4 ms later the collapse of cavitation bubbles between club and shell delivers a second force peak.

Source: Patek & Caldwell (2005), Journal of Experimental Biology

Four outcomes

No true effect True effect
p < .05 Type I error (α, usually 5%) correct: the probability is the power
p ≥ .05 correct Type II error (β)

Power = 1 − β: the probability that a test rejects the null hypothesis when an effect of a given size exists. The usual target of 80% is an arbitrary convention (Green & MacLeod, 2016).

What power depends on

  • The size of the effect we want to detect
  • The noise around it: residual variance, and in mixed models the variances of participants and items
  • The number of participants, the number of items, and observations per participant and item
  • The α level and the test

Only the sample sizes are under our control once the design and the effect of interest are fixed.

Low power distorts estimates

2,000 simulated studies with 20 participants per group. Only 15% are significant, and the significant positive estimates average 0.79 (dashed), 2.6 times the true effect (Gelman & Carlin, 2014).

The significance filter

If only significant results get published, published effects from underpowered studies are too large. Vasishth et al. (2018) ran seven direct replications of a psycholinguistic finding, with 268 participants in total. The original estimates turned out too noisy to be informative, and a larger-sample study gave a smaller effect.

A power analysis based on a published effect size then promises more power than the next study will have.

Power by simulation

Did you know: peak impact forces of peacock mantis shrimp strikes reached 1,501 newtons in one study, thousands of times the animal’s body weight.

Source: Patek & Caldwell (2005), Journal of Experimental Biology

The recipe

  1. Write down the model, with values for every parameter: fixed effects, random-effect variances, residual variance
  2. Simulate a data set from it, with the planned sample size
  3. Fit the model to the simulated data and run the test
  4. Repeat many times. The proportion of significant tests estimates the power

Because the effect exists in every simulated data set, every non-significant result is a Type II error (Green & MacLeod, 2016).

A first example

A two-group comparison with 30 participants per group and a true difference of 0.5 SD:

sim_p <- function(n, d) t.test(rnorm(n, d), rnorm(n, 0), var.equal = TRUE)$p.value
set.seed(1)
p <- replicate(2000, sim_p(n = 30, d = 0.5))
mean(p < 0.05)
[1] 0.483
power.t.test(n = 30, delta = 0.5)$power
[1] 0.477841

The simulation agrees with the formula, up to simulation error.

A power curve

For d = 0.5, about 64 participants per group reach 80% power.

Why simulate?

  • For a t-test, a formula is quicker. For participants and items as crossed random effects, random slopes, binary outcomes or GAMMs, there is no simple formula
  • Simulation works for any model that can be fitted, with the analysis we will actually run (Green & MacLeod, 2016; Brysbaert & Stevens, 2018)

Writing the simulation also forces us to state every assumption about effect sizes, variances and the design.

Power with {simr}

Did you know: some mantis shrimp have 12 types of colour photoreceptors, against three in humans. In tests, they still told colours apart worse than humans do.

Source: Thoen, How, Chiou & Marshall (2014), Science

The plan

We want to replicate the lexical decision study from lecture 1, and the question is the interaction: does word frequency matter more for L2 than for L1 speakers?

  • {simr} simulates from an lmer() or glmer() model, refits it and tests a fixed effect, many times (Green & MacLeod, 2016). The model can come from a pilot study, or be built from scratch with makeLmer()
  • fixef<- sets the effect size; extend() and powerCurve() vary the sample size

A pilot study

set.seed(2035)
participants <- tibble(participant = factor(sprintf("P%02d", 1:30)),
                       group = rep(c(-0.5, 0.5), 15),
                       p0 = rnorm(30, 0, 60), p1 = rnorm(30, 0, 12))
words <- tibble(word = factor(sprintf("w%02d", 1:40)),
                freq_c = runif(40, 1, 6) - 3.5, w0 = rnorm(40, 0, 30))
pilot <- crossing(participants, words) |>
  mutate(rt = 750 + p0 + w0 + 120 * group + (-30 - 20 * group + p1) * freq_c + rnorm(n(), 0, 70)) |>
  select(participant, group, word, freq_c, rt)
m_pilot <- lmer(rt ~ freq_c * group + (1 + freq_c | participant) + (1 | word), data = pilot)
round(fixef(m_pilot), 1)
 (Intercept)       freq_c        group freq_c:group 
       756.4        -29.2        127.5        -29.1

group is coded −0.5 (L1) and 0.5 (L2), so freq_c:group is the difference between the two frequency slopes. The groups alternate, so any first k participants are balanced.

Choosing the effect size

The pilot estimate of the interaction is noisy, and it may be large by chance. Planning with it repeats the significance filter.

Instead, decide on the smallest effect of interest: the smallest difference that would matter for the theory. Here we say the slopes should differ by at least 10 ms per log unit.

fixef(m_pilot)["freq_c:group"] <- -10

Retrospective “observed power”, with the effect size taken from the data, gives misleading results (Hoenig & Heisey, 2001; Green & MacLeod, 2016). More on this at the end.

Power for the pilot design

ps30 <- powerSim(m_pilot, test = fixed("freq_c:group", "z"), nsim = 200, seed = 1)
ps30
Power for predictor 'freq_c:group', (95% confidence interval):
      51.00% (43.85, 58.12)

Test: z-test
      Effect size for freq_c:group is -10.

Based on 200 simulations, (0 errors, 3 warnings, 0 messages)

Time elapsed: 0 h 0 m 9 s

51% power for 30 participants and 40 words. The z-test is fast, which suits exploring designs, but it is too liberal with few participants. For the final calculation use "sa" (Satterthwaite) or "kr" (Kenward–Roger), and 1,000 simulations or more.

More participants or more words?

m_part <- extend(m_pilot, along = "participant", n = 80)
pc_part <- powerCurve(m_part, test = fixed("freq_c:group", "z"), along = "participant",
                      breaks = c(10, 20, 30, 40, 60, 80), nsim = 100, seed = 1)
m_word <- extend(m_pilot, along = "word", n = 120)
pc_word <- powerCurve(m_word, test = fixed("freq_c:group", "z"), along = "word",
                      breaks = c(20, 40, 60, 80, 120), nsim = 100, seed = 1)

extend() adds new participants (or words) by copying the design of existing ones, and powerCurve() fits the model to the first 10, 20, … of them.

Two power curves

The interaction compares groups of participants, and participants differ in their frequency slopes. Extra participants raise power much faster than extra words, which only make each participant’s slope a little more precise.

From scratch

Without a pilot study, makeLmer() builds the model from values we choose, for example from earlier studies:

design <- select(pilot, participant, group, word, freq_c)
m_scratch <- makeLmer(rt ~ freq_c * group + (1 + freq_c | participant) + (1 | word),
                      fixef = c(750, -30, 120, -10),
                      VarCorr = list(participant = matrix(c(60^2, 0, 0, 12^2), 2), word = 30^2),
                      sigma = 70, data = design)

The random-effect variances matter as much as the fixed effect. Because they are estimates too, try a plausible range and see how much the power changes.

Post hoc power

Did you know: mantis shrimp were the first animals shown to see circularly polarised light. Some species also reflect it from body parts they display in courtship.

Source: Chiou et al. (2008), Current Biology

What post hoc power is

Post hoc power, also called observed or retrospective power, is computed after a study, with the effect size estimated from its own data. Reviewers often ask for it (Heinsberg & Weeks, 2022).

The calculation is the same as in a power analysis. The next slides show why the result is still uninformative.

Two universities

Heinsberg and Weeks (2022) simulated the same mouse experiment, run in parallel at two universities (n = 100 each), with different true effects.

Estimate (95% CI) p Post hoc power True effect True power
A 0.290 (−0.002, 0.582) .052 0.584 0.402 0.841
B 0.275 (−0.002, 0.553) .052 0.537 0.201 0.322

The two studies look almost identical, and so does their post hoc power. Their true power is very different, and nothing in the data reveals it.

It only restates the p-value

1,000 studies of the same design, all with 28% true power. Post hoc power falls on one curve of the p-value: about 50% at p = .05, and below 50% whenever p > .05 (Hoenig & Heisey, 2001).

Why it is uninformative

  • It treats the observed effect as the true effect, which is exactly what we do not know
  • It adds nothing to the p-value, and a non-significant result always comes with low post hoc power
  • It cannot tell a study that missed a real effect from a study of an effect that does not exist
  • Running powerSim() on your final model without changing the effect size computes the same thing

A reviewer asks

You report a non-significant interaction. A reviewer asks for a post hoc power analysis “to show whether the study was underpowered”.

How would you answer? What could you report instead that would help the reader judge the null result?

What to report instead

  • The estimate with its confidence interval: it shows which effect sizes the data are compatible with
  • The power analysis you ran before the study, with the effect size you chose and why
  • A sensitivity analysis: the smallest effect the design had 80% power to detect. It uses the design and assumed variances, not the observed effect

Reporting a power analysis

We determined the sample size by simulation with {simr} (Green & MacLeod, 2016), using a linear mixed model with the structure of the planned analysis. The smallest interaction of interest was a difference of 10 ms per log unit between the frequency slopes of the two groups; the random-effect variances came from a pilot study. With 40 words, 60 participants gave an estimated power of 79%.

Also report the test, the number of simulations (1,000 or more in a real study) and how you chose every value.

Summary

  • Power is the probability of detecting an effect of a given size with a given design
  • Underpowered studies miss real effects, and their significant results exaggerate them
  • Simulation estimates power for any model: simulate, refit, test, repeat
  • With {simr}, set the smallest effect of interest with fixef<-, then vary participants and items with extend() and powerCurve()
  • Plan power before the study. Post hoc power is a restatement of the p-value; report confidence intervals instead

References

  • Brysbaert, M., & Stevens, M. (2018). Power analysis and effect size in mixed effects models: A tutorial. Journal of Cognition, 1(1), 9. https://doi.org/10.5334/joc.10
  • Gelman, A., & Carlin, J. (2014). Beyond power calculations: Assessing Type S (sign) and Type M (magnitude) errors. Perspectives on Psychological Science, 9(6), 641–651.
  • Green, P., & MacLeod, C. J. (2016). SIMR: An R package for power analysis of generalized linear mixed models by simulation. Methods in Ecology and Evolution, 7(4), 493–498. https://doi.org/10.1111/2041-210X.12504
  • Heinsberg, L. W., & Weeks, D. E. (2022). Post hoc power is not informative. Genetic Epidemiology, 46(7), 390–394. https://doi.org/10.1002/gepi.22464
  • Hoenig, J. M., & Heisey, D. M. (2001). The abuse of power: The pervasive fallacy of power calculations for data analysis. The American Statistician, 55(1), 19–24.
  • Vasishth, S., Mertzen, D., Jäger, L. A., & Gelman, A. (2018). The statistical significance filter leads to overoptimistic expectations of replicability. Journal of Memory and Language, 103, 151–175. https://doi.org/10.1016/j.jml.2018.07.004

Sources of the animal facts

  • Chiou, T.-H., Kleinlogel, S., Cronin, T., Caldwell, R., Loeffler, B., Siddiqi, A., Goldizen, A., & Marshall, J. (2008). Circular polarization vision in a stomatopod crustacean. Current Biology, 18(6), 429–434. https://doi.org/10.1016/j.cub.2008.02.066
  • Patek, S. N., & Caldwell, R. L. (2005). Extreme impact and cavitation forces of a biological hammer: Strike forces of the peacock mantis shrimp Odontodactylus scyllarus. Journal of Experimental Biology, 208(19), 3655–3664. https://doi.org/10.1242/jeb.01831
  • Patek, S. N., Korff, W. L., & Caldwell, R. L. (2004). Deadly strike mechanism of a mantis shrimp. Nature, 428(6985), 819–820. https://doi.org/10.1038/428819a
  • Thoen, H. H., How, M. J., Chiou, T.-H., & Marshall, J. (2014). A different form of color vision in mantis shrimp. Science, 343(6169), 411–413. https://doi.org/10.1126/science.1245824

Slides and code

QR code linking to these slides

These slides

mcanzi.de/ss27_advanced_statistics/08_power_analysis.html

Code and all lectures

github.com/mcanzi/ss27_advanced_statistics