Seven sites, one question: do these possums differ in size?
One factor, more than two groups, done properly from end to end on real data found on Kaggle: 104 Australian brushtail possums trapped at seven forest sites. Why one test instead of twenty-one t-tests, naming and ordering groups with factor(), fitting lm() and reading car::Anova(type = 2), checking the assumptions AFTER the fit with performance::check_model() and car::leveneTest(), measuring how big the effect is with eta-squared, and finding out WHICH sites differ with emmeans and compact letter displays.
lm() fit a model and gave us slopes and intercepts✅ Transition
Today the question is the same shape — do these groups differ? — but there are seven groups, not two. That one change breaks the t-test, and the thing that replaces it is ANOVA. 🐾
factor() to give groups names and an orderlm(), test with car::Anova(type = 2)performance::check_model()car::leveneTest()emmeans + lettersTools today:
tidyversecar — Anova(), leveneTest()performance — check_model()emmeans — follow-up testsmultcomp — the lettersData: 🐾 104 brushtail possums, 7 sites, off Kaggle
Five chunks: why → look → test → check → which.
Five short chunks. After each one you switch to the activity and run the code yourself.
Every line is already written for you. Your job is to run it, read the output, and argue with it.
⚠️ Install first
install.packages("performance")
install.packages("see")
install.packages("emmeans")
install.packages("multcomp")performance pulls in a dozen dependencies. Do it before class, not during.
The homework is a separate file
19_anova_one_factor_homework.R — the same possums, measured a different way.
We will cover: where the data came from, and the one arithmetic fact that kills the obvious approach.
Your final projects need data. So here is how today’s data was actually found.
Searched kaggle.com/datasets for possum. One hit:
So the trail runs: Kaggle → textbook → peer-reviewed field study.
📖 Always walk the trail back
Kaggle is a mirror, not a source. Anyone can upload anything.
Before you trust a dataset, find out who collected it. If you cannot, that is your answer.
Mountain brushtail possum, Trichosurus cunninghami 🐾
Our question today:
Do possums from different sites differ in total length?
🧪 One factor, seven levels
total_l, total length in cm — a numbersite, which forest they were trapped in — seven groupsOne factor. That is the whole design, and the name of the test.
Seven groups. We know how to compare two. So: compare every pair?
Seven groups make 21 pairs. Twenty-one t-tests.
Each test at \(\alpha = 0.05\) has a 5% chance of a false positive. Run 21 of them on data where nothing is going on:
\[1 - 0.95^{21} = 0.66\]
A 66% chance of at least one “significant” result from pure noise.
⚠️ This is the whole problem
You would find something. You would find it two times out of three. And it would mean nothing.
More tests is not more rigour. More tests is more chances to fool yourself.
ANOVA asks one question first:
Is there any difference among these seven means?
One question, one p-value, one 5% risk.
It answers by comparing two kinds of variation:
\[F = \frac{\text{variation between groups}}{\text{variation within groups}}\]
🧠 What F means
“Analysis of variance” is a literal name.
We will cover: reading the file, turning site numbers into real names with factor(), and plotting before testing.
# Load packages ---------------------------------------
library(tidyverse)
library(car) # Anova(), leveneTest()
library(performance) # check_model()
library(emmeans) # follow-up comparisons
library(multcomp) # the compact letter display
# Read the file ---------------------------------------
possum_raw <- read_csv("data/possum.csv")
dim(possum_raw)[1] 104 8names(possum_raw)[1] "site" "pop" "sex" "age" "head_l" "skull_w"
[7] "total_l" "tail_l" multcompmultcomp loads MASS, and MASS has its own select(). It masks the dplyr one.
So after today’s library() block, this quietly breaks:
possum_df %>% select(site, total_l)
#> Error: unused argumentsThe fix is to say which one you mean:
possum_df %>% dplyr::select(site, total_l)⚠️ Watch the startup messages
R told you this when you loaded the package. Nobody reads it.
select() is the usual casualty. filter() is the other one.
None of today’s code uses select() — but your code might.
possum_raw %>% count(site)# A tibble: 7 × 2
site n
<dbl> <int>
1 1 33
2 2 13
3 3 7
4 4 7
5 5 13
6 6 13
7 7 18⚠️ Those are numbers, and R will treat them as numbers
site arrives as the integers 1–7. To R that is a quantity — as if site 6 were “twice site 3”. It is not. It is a name that happens to be written with a digit.
Leave it like this and lm() will fit a straight line through site number, which is nonsense.
factor()possum_df <- possum_raw %>%
filter(!is.na(total_l)) %>%
mutate(site_name = factor(site,
levels = 1:7,
labels = c("Cambarville", "Bellbird",
"Whian Whian", "Byrangery",
"Conondale", "Allyn River",
"Bulburin")))
possum_df %>% count(site_name)# A tibble: 7 × 2
site_name n
<fct> <int>
1 Cambarville 33
2 Bellbird 13
3 Whian Whian 7
4 Byrangery 7
5 Conondale 13
6 Allyn River 13
7 Bulburin 18factor(site,
levels = 1:7,
labels = c("Cambarville", ...))levels = the values as they appear in the datalabels = what you want them calledThe two lists must be the same length and in the same order. Get the order wrong and every single result is quietly mislabelled — no error, no warning.
The site names come from the OpenIntro codebook. We did not invent them.
📋 The seven sites
⚠️ Unequal groups
33 possums at Cambarville, 7 at Byrangery. That is fine for ANOVA — but it is why we will not hand-calculate anything.
ggplot(possum_df, aes(x = site_name, y = total_l)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(width = 0.15, alpha = 0.5)
Before a single test:
So we expect: a significant overall test, driven by the ends, with the middle a mess.
🎯 Why this matters
Predicting the answer before you test is not cheating. It is the only way to notice when the output is wrong.
If the test comes back saying nothing differs, you now know to go looking for a mistake.
possum_df %>%
group_by(site_name) %>%
summarise(n = n(),
mean_total = round(mean(total_l), 1),
sd_total = round(sd(total_l), 2))# A tibble: 7 × 4
site_name n mean_total sd_total
<fct> <int> <dbl> <dbl>
1 Cambarville 33 89.7 2.99
2 Bellbird 13 81.8 3.9
3 Whian Whian 7 88.1 3.35
4 Byrangery 7 92.2 2.48
5 Conondale 13 86.9 3.01
6 Allyn River 13 84.5 3.65
7 Bulburin 18 85.7 3.12Seven sites on a north–south gradient. Tempting story:
Colder in the south → bigger animals. Bergmann’s rule.
Look again. The two Victorian sites are the southernmost — and they are at opposite extremes:
Two sites, same state, 8 cm apart. Latitude cannot explain that.
⚠️ Do not let the test write your story
A significant F tells you the sites differ. It says nothing about why.
“Sites differ” is the finding. “Because of latitude” would be a second claim, and this data does not support it.
We will cover: lm() for more than two groups, then car::Anova() and how to read its table.
lm() — Same Function as Alwaystotal_model <- lm(total_l ~ site_name, data = possum_df)🧠 Nothing new here
response ~ predictor, exactly as in regression. The only change is that the predictor is a factor with seven levels instead of a number.
lm() does not need a different function for ANOVA. ANOVA is a linear model. That is not a simplification — it is literally the same machinery.
summary() Firstsummary(total_model) # six rows of "site_nameBellbird" etc.summary() gives you six slopes, each comparing one site to Cambarville — because Cambarville is level 1, so it became the baseline.
That answers “is each site different from Cambarville?”, which is not the question. Nobody asked about Cambarville.
⚠️ The wrong table
Six p-values, all against an arbitrary reference group. Change which site is level 1 and all six change.
We want one p-value for site_name as a whole.
car::Anova()Anova(total_model, type = 2)Anova Table (Type II tests)
Response: total_l
Sum Sq Df F value Pr(>F)
site_name 907.59 6 14.582 8.452e-12 ***
Residuals 1006.24 97
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1✅ This is the result
F(6, 97) = 14.58, p = 8.45 × 10⁻¹²
One line. One p-value. One question answered: the seven site means are not all the same.
Sum Sq 907.6 — variation explained by siteSum Sq 1006.2 — variation left over within sitesDf 6 — seven groups minus oneDf 97 — 104 possums minus 7 groupsF value 14.58 — the ratio, after dividing each Sum Sq by its DfPr(>F) — how often an F this big happens when nothing is going on\(F = 14.58\) means between-site variation is 14 and a half times within-site variation.
🔢 Why type = 2?
With one factor, Type I, II and III all give the identical answer. It genuinely does not matter here.
We write type = 2 anyway, as a habit, because the moment you add a second factor it starts mattering a great deal.
effectsize::eta_squared(Anova(total_model, type = 2), partial = FALSE)# Effect Size for ANOVA (Type II)
Parameter | Eta2 | 95% CI
-------------------------------
site_name | 0.47 | [0.34, 1.00]
- One-sided CIs: upper bound fixed at [1.00].Eta-squared (\(\eta^2\)) = 0.47.
Which site a possum came from accounts for 47% of all the variation in total length.
The other 53% is possums being individuals.
🎯 p ≠ size
p answers “is it real?” \(\eta^2\) answers “is it big?”
With a big enough sample, a trivial difference gets a tiny p. Always report both.
We will cover: what ANOVA assumes, why you check it after fitting, and two packages that do the work.
Plus one you check by thinking, not testing:
⚠️ Order matters
You cannot check residuals before you have residuals. So:
fit first, then check.
If the check fails you go back and change the model — you do not delete the data.
check_model()check_model(total_model, check = c("homogeneity", "qq"))
Those are the two panels that match our two assumptions:
Both look healthy here. That is a boring plot, and boring is exactly what you want.
📺 On your screen you get five
Plain check_model(total_model) draws all five:
We asked for two here only so they fit on a slide.
check_model() Throws an ErrorFive panels need room. If the RStudio Plots pane is small, check_model() fails with a complaint about the plot size rather than drawing anything.
Two fixes:
check_model(total_model,
check = c("homogeneity", "qq"))⚠️ Not your code
This error is about your window, not your model. Nothing is wrong with the analysis.
It catches everybody once.
# Equal variances across sites
leveneTest(total_l ~ site_name, data = possum_df)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 6 0.5032 0.8046
97 # Normal residuals
shapiro.test(residuals(total_model))
Shapiro-Wilk normality test
data: residuals(total_model)
W = 0.98939, p-value = 0.588check_homogeneity(total_model)OK: There is not clear evidence for different variances across groups (Bartlett Test, p = 0.867).check_normality(total_model)OK: residuals appear as normally distributed (p = 0.602).✅ Both assumptions pass
High p-values are the good outcome here. You are hoping not to reject.
For the main test, a small p was the interesting result.
For assumption tests it is the opposite — a small p means your test is in trouble.
That inversion catches people constantly. Say it out loud each time:
“Levene p = 0.805. Big. Variances are fine. Carry on.”
⚠️ And do not over-trust them
With small groups these tests are weak — they pass almost anything.
With huge samples they are brutal — they fail on differences too small to matter.
Look at check_model() as well, always.
We will cover: why F is not the end, emmeans, Tukey adjustment, and reading a compact letter display.
We have F(6, 97) = 14.58, p < 0.001.
That is: at least one site differs from at least one other.
It does not tell you:
For that you need post-hoc comparisons — “after the fact”.
🧠 Why this is allowed now
Doing 21 comparisons was wrong as a first move.
Doing them after a significant F, with the p-values adjusted, is the standard, defensible approach.
The gate came first.
emmeanssite_means <- emmeans(total_model, ~ site_name)
site_means site_name emmean SE df lower.CL upper.CL
Cambarville 89.7 0.561 97 88.6 90.8
Bellbird 81.8 0.893 97 80.0 83.5
Whian Whian 88.1 1.220 97 85.7 90.5
Byrangery 92.2 1.220 97 89.8 94.7
Conondale 86.9 0.893 97 85.2 88.7
Allyn River 84.5 0.893 97 82.8 86.3
Bulburin 85.7 0.759 97 84.2 87.2
Confidence level used: 0.95 pairs(site_means, adjust = "tukey") contrast estimate t.ratio p.value
1 Cambarville - Bellbird 7.94 7.53 < 1e-04
2 Cambarville - Whian Whian 1.64 1.22 0.88310
3 Cambarville - Byrangery -2.53 -1.89 0.49301
4 Cambarville - Conondale 2.79 2.64 0.12429
5 Cambarville - Allyn River 5.17 4.91 < 1e-04
6 Cambarville - Bulburin 4.02 4.26 0.00092
7 Bellbird - Whian Whian -6.30 -4.17 0.00124
8 Bellbird - Byrangery -10.47 -6.94 < 1e-04First 8 of 21 rows. Tukey widens every p-value to pay for running 21 tests.
Of the 21 pairs, 10 are significant and 11 are not.
Which is exactly what the boxplot predicted:
Twenty-one rows is a lot to read. There is a better way to show it.
🎯 Report the pattern
Nobody wants 21 p-values in a results section.
They want to know which groups cluster together. That is what letters do.
cld(site_means, Letters = letters) site_name emmean SE df lower.CL upper.CL .group
Bellbird 81.8 0.893 97 80.0 83.5 a
Allyn River 84.5 0.893 97 82.8 86.3 ab
Bulburin 85.7 0.759 97 84.2 87.2 b
Conondale 86.9 0.893 97 85.2 88.7 bc
Whian Whian 88.1 1.220 97 85.7 90.5 bcd
Cambarville 89.7 0.561 97 88.6 90.8 cd
Byrangery 92.2 1.220 97 89.8 94.7 d
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 7 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Share a letter → not shown to differ. No letter in common → differ.
a — alone at the bottomab — not distinguishable from either neighbourd — alone at the topa vs d, nothing shared → differA clean ladder from 81.8 cm up to 92.2 cm.
⚠️ Read R’s own note
“If two or more means share the same grouping symbol, then we cannot show them to be different. But we also did not show them to be the same.”
Sharing a letter is failure to detect a difference — not evidence of equality. Especially at Byrangery, where n = 7.
letters_df <- as.data.frame(cld(site_means, Letters = letters)) %>%
mutate(group = trimws(.group))
ggplot(possum_df, aes(x = site_name, y = total_l)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(width = 0.15, alpha = 0.4) +
geom_text(data = letters_df,
aes(x = site_name, y = 97, label = group),
size = 5)
📝 The whole result, in three sentences
Total length differed significantly among the seven trapping sites (one-way ANOVA, F(6, 97) = 14.58, p < 0.001), and site accounted for 47% of the variation in length (η² = 0.47). Residuals were normally distributed (Shapiro–Wilk p = 0.59) with homogeneous variances across sites (Levene’s test p = 0.81). Tukey-adjusted comparisons separated Bellbird (mean 81.8 cm) from Byrangery (92.2 cm), with the remaining five sites forming an overlapping middle group.
Every number has a job:
F(6, 97); the df let a reader reconstruct your designp < 0.001, not p = 0⚠️ Never write “p = 0.000”
Software rounds. Probability does not reach zero.
Write p < 0.001 and move on.
summary() instead of Anova() → six wrong p-values🎯 The recipe, memorised
factor() the groupslm()Anova(type = 2)check_model() + leveneTest()eta_squared()emmeans() → pairs() → cld()The file has four measurements on every animal. We used total_l. For homework you run the whole recipe again on tail_l — tail length.
Same 104 animals, same seven sites, same eight steps. A different ruler.
The question is whether the pattern holds: do the sites that had the longest bodies also have the longest tails?
📂 What you will do
All eight steps on tail_l, then one paragraph comparing your letters to the ones we got today.
Plus one optional part where an assumption fails on purpose.
Next time: New Zealand tree wētā, a large flightless cricket relative whose adult males come in three distinct forms — small-jawed sneaks, intermediates, and huge-jawed fighters.
We will find that the three forms differ enormously in ear size.
Then we will divide ear size by body size, run the identical test, and the answer will change.
🎯 Why that matters
ANOVA answers the question you asked it.
Making sure that is the question you meant is the whole job of next class.
factor(levels =, labels =) turns group codes into named, ordered groupslm() fits it — ANOVA is a linear model; car::Anova(type = 2) gives the one right tablecheck_model(), leveneTest(), shapiro.test(); here a big p is goodemmeans() → pairs() → cld() says which — with p-values adjusted✅ Next
Open the activity and run all of it. Homework is tail_l. Next class, the wētā. 🦗