library(tidyverse)
library(broom)
library(emmeans)
set.seed(42)
theme_set(theme_minimal(base_size = 12))Week 2, Session 1 — One-way ANOVA and contrasts (emmeans)
Course 2 — #courses
Note
Inference labs use the five-step template: Hypothesis → Visualise → Assumptions → Conduct → Conclude.
Learning objectives
- Fit a one-way ANOVA with
aov()and interpret the omnibus F-test. - Use
emmeansto obtain estimated marginal means and pairwise contrasts with controlled family-wise error. - Read a Tukey HSD output and translate it into a sentence.
Prerequisites
Course 1 t-tests; Course 2 Week 1 regression basics.
Background
Analysis of variance is linear regression with a categorical predictor and a tradition. The omnibus F-test asks whether any group means differ; the follow-up contrasts ask which ones. The correct order is omnibus first, then targeted contrasts, preferably with a multiple testing correction that respects the family of tests you intended before seeing the data.
The emmeans package has made the post-hoc machinery substantially cleaner. emmeans(fit, ~ group) returns the estimated marginal means on the response scale; pairs() produces the pairwise differences with Tukey-adjusted p-values by default.
A common mistake is to run pairwise t-tests without correction and report the one with the smallest p-value. A less common but more insidious mistake is to run the omnibus F, see it is significant, and stop — reporting a single p-value without naming which groups differ and by how much.
Setup
1. Hypothesis
Using ToothGrowth, ask whether guinea-pig tooth length depends on dose of vitamin C.
Null: mean tooth length is equal across the three doses. Alternative: at least one dose differs.
2. Visualise
tg <- ToothGrowth |> mutate(dose = factor(dose))
ggplot(tg, aes(dose, len, fill = dose)) +
geom_boxplot(alpha = 0.6, colour = "grey30") +
labs(x = "Dose (mg/day)", y = "Tooth length") +
theme(legend.position = "none")
3. Assumptions
Approximate normality within groups and equal variances.
fit <- aov(len ~ dose, data = tg)
par(mfrow = c(2, 2))
plot(fit)
par(mfrow = c(1, 1))
bartlett.test(len ~ dose, data = tg)
Bartlett test of homogeneity of variances
data: len by dose
Bartlett's K-squared = 0.66547, df = 2, p-value = 0.7174. Conduct
summary(fit) Df Sum Sq Mean Sq F value Pr(>F)
dose 2 2426 1213 67.42 9.53e-16 ***
Residuals 57 1026 18
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm <- emmeans(fit, ~ dose)
emm dose emmean SE df lower.CL upper.CL
0.5 10.6 0.949 57 8.71 12.5
1 19.7 0.949 57 17.84 21.6
2 26.1 0.949 57 24.20 28.0
Confidence level used: 0.95 pairs(emm, adjust = "tukey") contrast estimate SE df t.ratio p.value
dose0.5 - dose1 -9.13 1.34 57 -6.806 <0.0001
dose0.5 - dose2 -15.49 1.34 57 -11.551 <0.0001
dose1 - dose2 -6.37 1.34 57 -4.745 <0.0001
P value adjustment: tukey method for comparing a family of 3 estimates 5. Concluding statement
Tooth length increased with dose of vitamin C (one-way ANOVA F = 67.42; p = 9.5^{-16}). Tukey-adjusted contrasts indicated that the 2 mg/day dose gave longer teeth than both lower doses (all adjusted p < 0.05).
Mention that aov() and lm() fit the same model; the choice is mostly about what functions you want on top.
Common pitfalls
- Stopping at a significant F without reporting which contrasts drive it.
- Reporting unadjusted pairwise p-values after the omnibus test.
- Treating a non-significant Bartlett test as proof of equal variances.
Further reading
- Lenth RV. emmeans package vignette.
- Oehlert GW. A First Course in Design and Analysis of Experiments.
- Maxwell SE, Delaney HD. Designing Experiments and Analyzing Data.
Session info
sessionInfo()R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.4 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
locale:
[1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
[4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
[7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
time zone: UTC
tzcode source: system (glibc)
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] emmeans_2.0.3 broom_1.0.12 lubridate_1.9.5 forcats_1.0.1
[5] stringr_1.6.0 dplyr_1.2.1 purrr_1.2.2 readr_2.2.0
[9] tidyr_1.3.2 tibble_3.3.1 ggplot2_4.0.3 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] gtable_0.3.6 jsonlite_2.0.0 compiler_4.4.1 tidyselect_1.2.1
[5] scales_1.4.0 yaml_2.3.12 fastmap_1.2.0 lattice_0.22-6
[9] coda_0.19-4.1 R6_2.6.1 labeling_0.4.3 generics_0.1.4
[13] knitr_1.51 backports_1.5.1 htmlwidgets_1.6.4 pillar_1.11.1
[17] RColorBrewer_1.1-3 tzdb_0.5.0 rlang_1.2.0 stringi_1.8.7
[21] xfun_0.57 S7_0.2.2 otel_0.2.0 estimability_1.5.1
[25] timechange_0.4.0 cli_3.6.6 withr_3.0.2 magrittr_2.0.5
[29] digest_0.6.39 grid_4.4.1 xtable_1.8-8 mvtnorm_1.3-7
[33] hms_1.1.4 lifecycle_1.0.5 vctrs_0.7.3 evaluate_1.0.5
[37] glue_1.8.1 farver_2.1.2 rmarkdown_2.31 tools_4.4.1
[41] pkgconfig_2.0.3 htmltools_0.5.9