Power and model building

p-values, power, uncertainty and building models

Dr Massimiliano Canzi

What p-values are not

Did you know: slow lorises are the only known venomous primates. Glands on the inside of their elbows make a secretion that turns venomous when mixed with saliva.

Source: San Diego Zoo Wildlife Alliance

Today

The key points of three chapters of Winter (2019):

Chapter	Topics
10	what p-values are not; power; Type I, II, M and S errors; multiple testing; stopping rules
11	confidence intervals, testing predictors with several levels, plotting uncertainty
16	strategies for model building, reproducible research

Most of today is simulation: we make up a world where we know the truth, and see how our methods behave in it.

So far

Lecture 9 introduced p-values, standard errors and confidence intervals, and showed that significance says nothing about size. Lectures 10 to 12 added mixed models with lmer() and glmer(), tested with likelihood ratio tests.

Today: how p-values can mislead, and how to plan and report models more carefully.

Common misreadings

A p-value of .03 does not mean:

that there is a 3% chance that the null hypothesis is true
that there is a 97% chance that the effect is real
that the effect is large or important

It means: if the null hypothesis were true, data at least this extreme would turn up 3% of the time.

Absence of evidence

p > .05 does not mean that there is no effect.

The data are compatible with no effect, and they may be just as compatible with a large one, so look at the confidence interval. A small study can also miss a real effect; that is a question of power, which comes next.

Significant vs not significant

Two studies of the same effect (Gelman & Stern, 2006):

	Estimate	SE	p
Study A	25	10	.01	significant
Study B	10	10	.32	not significant

Do the two studies disagree? Test the difference:

diff    <- 25 - 10
se_diff <- sqrt(10^2 + 10^2)
2 * pnorm(-diff / se_diff)

[1] 0.2888444

The difference between them is not significant either. Compare the effects directly, not their p-values.

Power and errors

Did you know: a threatened slow loris clasps its arms over its head. With its face markings, the pose looks like the hood of an angry cobra.

Source: San Diego Zoo Wildlife Alliance

Two ways to be wrong

	No effect in reality	Effect in reality
Significant	Type I error (false alarm)	correct
Not significant	correct	Type II error (miss)

We choose the Type I error rate, α, and usually set it to .05. Power is the probability of finding an effect that is really there, 1 minus the Type II error rate. It grows with the effect size and the sample size, and shrinks as the noise grows.

Simulating a study

A world we control: group b’s vowels are 10 ms longer than group a’s, SD = 40 ms.

sim_study <- function(n, diff = 10, sd = 40) {
  d <- tibble(group    = rep(c("a", "b"), each = n),
              duration = rnorm(2 * n, mean = 200 + diff * (group == "b"), sd = sd))
  fit <- summary(lm(duration ~ group, data = d))$coefficients
  tibble(estimate = fit[2, 1], p = fit[2, 4])
}

set.seed(42)
one <- sim_study(n = 20)
one

# A tibble: 1 × 2
  estimate     p
     <dbl> <dbl>
1    -8.52 0.583

The true difference is +10 ms. This study, with 20 speakers per group, estimates -8.5 ms (p = 0.58): one study can be far off.

Run it 1,000 times

set.seed(42)
sims20 <- map(1:1000, \(i) sim_study(n = 20)) |> list_rbind()
mean(sims20$p < .05)

[1] 0.113

The power is about 11%: the effect is real, yet almost 9 studies in 10 miss it. A formula gives the same answer:

power.t.test(n = 20, delta = 10, sd = 40)$power

[1] 0.1171781

Power and sample size

For 80% power, this 10 ms effect needs about 253 speakers per group. Plan the sample size before collecting data.

Type M and Type S errors

Gelman & Carlin (2014): with low power, the significant results are the misleading ones.

Only the large estimates reach significance.

Type M and Type S errors

sig <- filter(sims20, p < .05)
mean(abs(sig$estimate))

[1] 30.47138

mean(sig$estimate < 0)

[1] 0.02654867

The true difference is 10 ms. The significant estimates are on average 3 times too large, a Type M (magnitude) error, and 2.7% of them point the wrong way, a Type S (sign) error. So a significant result from a small study probably exaggerates the effect.

Power for mixed models

For mixed models, there is no simple formula. Do what we just did:

Simulate data with your design: participants, items, plausible effect sizes and random effects (lectures 10 to 12)
Fit the model to each simulated data set
Count how often the effect is significant

The {simr} package automates this for lmer() and glmer() models (Green & MacLeod, 2016).

Multiple tests and stopping rules

Did you know: slow lorises also use their venom on each other: in one study, 57% of males and nearly a third of females showed signs of at least one venomous bite.

Sources: Nekaris et al. (2020), Current Biology; Mongabay

Multiple testing

Each test on a null effect has a 5% chance of a false alarm. With k tests:

tibble(k = c(1, 5, 10, 20)) |>
  mutate(at_least_one_false_alarm = 1 - 0.95^k)

# A tibble: 4 × 2
      k at_least_one_false_alarm
  <dbl>                    <dbl>
1     1                   0.0500
2     5                   0.226 
3    10                   0.401 
4    20                   0.642

If you test 20 effects that do not exist, you will probably find at least one “significant” result.

A simulation

1,000 studies, each running 10 t-tests on pure noise:

set.seed(3)
any_sig <- map_lgl(1:1000, \(i) {
  p <- map_dbl(1:10, \(j) t.test(rnorm(20), rnorm(20))$p.value)
  any(p < .05)
})
mean(any_sig)

[1] 0.38

About 38% of the studies find “something”.

Correcting for multiple tests

p <- c(0.003, 0.012, 0.031, 0.048, 0.21)
p.adjust(p, method = "bonferroni")

[1] 0.015 0.060 0.155 0.240 1.000

p.adjust(p, method = "holm")

[1] 0.015 0.048 0.093 0.096 0.210

Bonferroni multiplies each p-value by the number of tests; Holm does something similar, step by step, and loses less power.

Better still: test fewer things. Decide which comparisons answer your question before you look at the data.

Stopping rules

A tempting plan: test after every 10 participants, and stop as soon as p < .05.

20 studies of no effect: p-values wander, and some dip below .05 by chance.

Peeking inflates false alarms

peek <- function(max_n = 100, step = 10) {
  a <- rnorm(max_n); b <- rnorm(max_n)
  p <- map_dbl(seq(step, max_n, by = step), \(k) t.test(a[1:k], b[1:k])$p.value)
  any(p < .05)
}
set.seed(5)
false_alarm <- mean(map_lgl(1:1000, \(i) peek()))
false_alarm

[1] 0.197

Both groups come from the same distribution, and peek() returns TRUE if any of the 10 tests was below .05.

With 10 peeks, the false alarm rate is about 20%, not 5%. Fix the sample size in advance, and report how you chose it.

Discussion

After 20 participants, your result is p = .07. Your supervisor suggests testing 10 more and checking again.

What is the problem with this plan?
What could you do instead, and how would you report what happened honestly?

Uncertainty in regression

Did you know: in cool months, pygmy slow lorises can go into torpor, a short form of hibernation. One study recorded inactivity lasting up to 63 hours.

Source: San Diego Zoo Wildlife Alliance

Confidence intervals for coefficients

pengs <- penguins |> drop_na(body_mass_g, flipper_length_mm, species)
m <- lm(body_mass_g ~ flipper_length_mm, data = pengs)
coef(summary(m))

                     Estimate Std. Error   t value      Pr(>|t|)
(Intercept)       -5780.83136 305.814504 -18.90306  5.587301e-55
flipper_length_mm    49.68557   1.518404  32.72223 4.370681e-107

confint(m)

                        2.5 %      97.5 %
(Intercept)       -6382.35801 -5179.30471
flipper_length_mm    46.69892    52.67221

Each mm of flipper goes with about 50 g of body mass, and the 95% CI runs from 46.7 to 52.7 g: roughly the estimate ± 2 SE.

Predictors with several levels

m_sp <- lm(body_mass_g ~ species, data = pengs)
coef(summary(m_sp))

                   Estimate Std. Error    t value      Pr(>|t|)
(Intercept)      3700.66225   37.61935 98.3712321 2.488024e-251
speciesChinstrap   32.42598   67.51168  0.4803018  6.313226e-01
speciesGentoo    1375.35401   56.14797 24.4951686  5.420612e-77

Two coefficients, two treatment contrasts against Adelie (lecture 7). But is species as a whole a useful predictor?

Testing the whole predictor

Compare models with and without it, as in lectures 11 and 12:

m_0 <- lm(body_mass_g ~ 1, data = pengs)
anova(m_0, m_sp)

Analysis of Variance Table

Model 1: body_mass_g ~ 1
Model 2: body_mass_g ~ species
  Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
1    341 219307697                                  
2    339  72443483  2 146864214 343.63 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

For lm(), the comparison is an F-test; for lmer() and glmer(), a likelihood ratio test. Same idea: does the predictor improve the model?

Showing uncertainty: categories

new_sp <- tibble(species = levels(pengs$species))
preds  <- bind_cols(new_sp, predict(m_sp, newdata = new_sp, interval = "confidence"))
preds

# A tibble: 3 × 4
  species     fit   lwr   upr
  <chr>     <dbl> <dbl> <dbl>
1 Adelie    3701. 3627. 3775.
2 Chinstrap 3733. 3623. 3843.
3 Gentoo    5076. 4994. 5158.

interval = "confidence" adds the lower (lwr) and upper (upr) limits of the 95% CI.

Showing uncertainty: categories

ggplot(preds, aes(species, fit)) +
  geom_jitter(data = pengs, aes(y = body_mass_g), width = 0.15,
              colour = kingfisher, alpha = 0.4) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.15,
                colour = hornbill, linewidth = 1.2) +
  geom_point(colour = hornbill, size = 4) +
  labs(x = NULL, y = "Body mass (g)")

Showing uncertainty: continuous

Predict over a grid of flipper lengths, within the range of the data:

grid <- tibble(flipper_length_mm = seq(172, 231, by = 1))
grid <- bind_cols(grid, predict(m, newdata = grid, interval = "confidence"))
grid

# A tibble: 60 × 4
  flipper_length_mm   fit   lwr   upr
              <dbl> <dbl> <dbl> <dbl>
1               172 2765. 2669. 2861.
2               173 2815. 2721. 2908.
3               174 2864. 2774. 2955.
4               175 2914. 2826. 3002.
5               176 2964. 2878. 3049.
# ℹ 55 more rows

Showing uncertainty: continuous

ggplot(grid, aes(flipper_length_mm, fit)) +
  geom_point(data = pengs, aes(y = body_mass_g), colour = kingfisher, alpha = 0.4) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), fill = hornbill, alpha = 0.35) +
  geom_line(colour = hornbill, linewidth = 1.3) +
  labs(x = "Flipper length (mm)", y = "Body mass (g)")

The band shows the uncertainty in the line, which is wider at the edges. geom_smooth(method = "lm") draws the same.

Building models

Did you know: the illegal pet trade is one of the worst threats to slow lorises. Traffickers pull out their teeth so that they cannot bite.

Sources: San Diego Zoo Wildlife Alliance; Mongabay

What you have learned

Outcome	Repeated measures?	Model
continuous	no	lm()
continuous	yes	lmer()
binary (0/1)	no	glm(family = binomial)
binary (0/1)	yes	glmer(family = binomial)

We fitted glmer() in lecture 12; glm() is the same model without random effects. One framework: an outcome, a linear predictor (continuous and categorical predictors, interactions), a distribution, and random effects where the data are grouped. Counts would use Poisson regression (Winter, Ch 13).

The cookbook approach

Many textbooks teach a cookbook: two groups → t-test; three groups → ANOVA; two numbers → correlation.

Each of these tests is a special case of the linear model; lecture 9 showed this for the t-test. The cookbook hides the assumptions behind each recipe, and it has no recipe once the data have repeated measures or several predictors. Winter recommends thinking in terms of models instead.

Stepwise regression

Put in every predictor, and let the computer drop the useless ones:

one_step <- function() {
  noise <- as_tibble(matrix(rnorm(50 * 21), ncol = 21,
                            dimnames = list(NULL, c("y", paste0("x", 1:20)))))
  best <- step(lm(y ~ ., data = noise), trace = 0)
  p <- summary(best)$coefficients[-1, 4]
  sum(p < .05)
}

The data have 50 rows: an outcome y and 20 predictors, all random noise. y ~ . uses every other column as a predictor, and step() drops predictors one at a time, as long as that improves the AIC.

Stepwise regression finds noise

set.seed(11)
n_sig <- map_dbl(1:300, \(i) one_step())
mean(n_sig >= 1)

[1] 0.8466667

In 85% of the data sets, the final model keeps at least one “significant” predictor, though nothing is real. Its p-values ignore all the tests that went into choosing it, a hidden case of multiple testing. Choose predictors from your research question and theory instead.

Theory-driven modelling

Winter (Ch 16): modelling is subjective, and that is fine, as long as it is principled and transparent.

Decide the model before seeing the data: outcome, predictors, interactions, random effects
Keep confirmatory (planned) and exploratory (after the fact) analyses apart, and label them
Consider preregistering the analysis
Report the choices you made, and why

Reproducible research

Someone else, or you in a year, should get the same results from your data.

Practice	Instead of
the whole analysis as code, e.g. in R Markdown	clicking and copying results by hand
`set.seed()` before simulations or `sample()`	different random numbers on every run
R projects and relative paths	`C:/Users/me/Desktop/...`
`sessionInfo()` to record package versions	guessing which version gave the result
data and code shared with the paper, e.g. on the OSF	“available on request”

Your turn

Exercises

Use sim_study() to estimate the power to detect a 20 ms difference (SD = 40) with 20 and with 50 speakers per group. Check your answers with power.t.test()
A study reports these p-values for 6 tests: 0.001, 0.008, 0.02, 0.04, 0.045, 0.3. Which remain below .05 after a Holm correction?
Fit flipper_length_mm ~ species on pengs. Test species with a model comparison, and plot the predicted means with their 95% CIs

Solutions

set.seed(1)
mean(map_dbl(1:1000, \(i) sim_study(n = 20, diff = 20)$p) < .05)

[1] 0.331

mean(map_dbl(1:1000, \(i) sim_study(n = 50, diff = 20)$p) < .05)

[1] 0.717

power.t.test(n = c(20, 50), delta = 20, sd = 40)$power

[1] 0.3377084 0.6968888

Doubling the effect makes a huge difference: with 50 speakers per group, power is about 70%.

Solutions

p_ex <- c(0.001, 0.008, 0.02, 0.04, 0.045, 0.3)
p.adjust(p_ex, method = "holm")

[1] 0.006 0.040 0.080 0.120 0.120 0.300

Only the first two (0.006 and 0.04) stay below .05.

Solutions

pengs_f <- penguins |> drop_na(flipper_length_mm, species)
m_fl <- lm(flipper_length_mm ~ species, data = pengs_f)
anova(lm(flipper_length_mm ~ 1, data = pengs_f), m_fl)$`Pr(>F)`[2]

[1] 1.35171e-111

preds_fl <- bind_cols(new_sp, predict(m_fl, newdata = new_sp, interval = "confidence"))
ggplot(preds_fl, aes(species, fit)) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.15, colour = hornbill, linewidth = 1.2) +
  geom_point(colour = hornbill, size = 4) + labs(x = NULL, y = "Flipper length (mm)")

Summary

A p-value is not the probability that the null is true; p > .05 is not “no effect”; compare effects directly
Low power misses real effects, and exaggerates the ones it finds (Type M and S errors)
Simulate to find out the power of a design, including mixed models
Many tests and peeking at the data inflate false alarms: plan in advance
Report confidence intervals, and plot uncertainty
Build models from theory instead of stepwise selection; keep your analysis reproducible

References

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. https://doi.org/10.1177/1745691614551642
Gelman, A., & Stern, H. (2006). The difference between “significant” and “not significant” is not itself statistically significant. The American Statistician, 60, 328–331. https://doi.org/10.1198/000313006X152649
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. https://doi.org/10.1111/2041-210X.12504
Winter, B. (2019). Statistics for linguists: An introduction using R. Routledge.

Slides and code

These slides

mcanzi.de/ws26_statistics/13_power_and_models.html

Code, data and all lectures

github.com/mcanzi/ws26_statistics