Think back to Type I errors from DAPR1. When we conduct an hypothesis test, and we set \(\alpha = .05\), we will reject the null hypothesis \(H_0\) when we find a \(p < .05\).
Now remember what a \(p\)-value represents: the chance of observing a statistic at least as extreme as the one we do have, assuming the null hypothesis to be true. This means that if\(H_0\)is true, then we will still observe a \(p < .05\) 5% of the time. So, our chance of making a Type I error = the threshold (\(\alpha\)) at which below a \(p\)-value results in us rejecting \(H_0\).
We accept a 5% error rate overall. But this error rate applies to each statistical test we run. So if we conduct an analysis in which we plan on conducting lots of tests of different comparisons, the chance of an error being made increases substantially. Across all the tests performed, that chance will be much higher than our accepted risk level of 5%.
When to correct for multiple comparisons
Whenever you are running more than one test using the same data and using them both to draw conclusions. This could be:
two linear models
one linear model and one t-test
one linear model and one chi-squared test
…
In this scenario, the simplest correction to use is Bonferroni (see below).
(We see this a lot in Honours dissertations: many tests are run to address several related RQs, but the p-values are not corrected, so the level of risk is higher than our accepted 5%.)
library(palmerpenguins)# predict body mass based on species, sex, and their interactionmdl_cc <-lm(body_mass_g ~ species * sex, data = penguins)# test the simple effects of species for each sexpairs( emmeans::emmeans(mdl_cc, ~ species * sex), simple ="species")
sex = female:
contrast estimate SE df t.ratio p.value
Adelie - Chinstrap -158 64.2 327 -2.470 0.0377
Adelie - Gentoo -1311 54.4 327 -24.090 <0.0001
Chinstrap - Gentoo -1153 66.8 327 -17.250 <0.0001
sex = male:
contrast estimate SE df t.ratio p.value
Adelie - Chinstrap 105 64.2 327 1.630 0.2357
Adelie - Gentoo -1441 53.7 327 -26.850 <0.0001
Chinstrap - Gentoo -1546 66.2 327 -23.350 <0.0001
P value adjustment: tukey method for comparing a family of 3 estimates
Within each sex, we are running three statistical tests, and so the \(p\)-values of those tests need to be adjusted.
By default, pairs() uses the Tukey correction.
CautionHow does the Tukey correction work?
Tukey’s HSD compares the means of the members of one pair in a pairwise comparison.
The difference between group means is represented as a q-statistic.
The q-statistic is compared to a critical value from the studentised range distribution to decide whether the difference is significant.
You can also tell pairs() to use the Bonferroni correction, which divides the \(\alpha\) level by the number of tests / multiplies the \(p\)-value by the number of tests.
pairs( emmeans::emmeans(mdl_cc, ~ species * sex), simple ="species",adjust ='bonferroni')