Activity: Analysis of Variance (ANOVA)
One-way ANOVA
Worksheet: Analysis of Variance (ANOVA)
How to use this worksheet
Work through each part in order, at your own pace. Type every line of code yourself into a plain R script — do not copy-paste. Blocks marked ▶ Run this are code you should type and execute. Blocks marked ✏️ Your turn ask you to write, modify, or answer something. Boxes marked 🚀 If you finish early are optional bonus material that goes a bit further than what we covered in lecture.
Part 1 · Setup and Data
▶ Run this in your Script:
library(performance)
library(FSA)
library(car) # functions for regression diagnostics, ANOVA, and VIF
library(emmeans) # calculates and compares adjusted means from statistical models
library(tidyverse) # includes ggplot2, dplyr, tidyr, etc.
c_df <- tibble(
treatment = rep(c("Control", "Knees", "Eyes"), times = c(8, 7, 7)),
phase_shift = c(0.53, 0.36, 0.20, -0.37, -0.60, -0.64, -0.68, -1.27, # Control
0.73, 0.31, 0.03, -0.29, -0.56, -0.96, -1.61, # Knees
-0.78, -0.86, -1.35, -1.48, -1.52, -2.04, -2.83) # Eyes
)
c_df▶ Run this — visualize the data:
c_plot <- ggplot(c_df, aes(x = treatment, y = phase_shift, color = treatment)) +
geom_point(position = position_jitter(width = 0.1)) +
stat_summary(fun = mean, geom = "point", size = 5, shape = 18) +
stat_summary(fun.data = "mean_se", geom = "errorbar", width = 0.2)
c_plot▶ Run this — fit the ANOVA (as a linear model, exactly like a regression) and get the ANOVA table:
model_aov <- lm(phase_shift ~ treatment, data = c_df)
summary(model_aov)
Anova(model_aov)✏️ Your turn: What is the F-statistic and p-value? At α = 0.05, do we reject H₀? ________________________
🚀 If you finish early: Compute group means, SDs, and n by hand with group_by() + summarize(), and confirm they match what you’d expect from c_plot.
# Write your code here:Part 2 · Assumptions and Diagnostics
ANOVA has the same assumptions as the two-sample t-test, but applied to all k groups: random samples, normality of Y in each population, homogeneity of variance, and independence.
▶ Run this — the base R diagnostic plots:
par(mfrow = c(2, 2))
plot(model_aov)
par(mfrow = c(1, 1))▶ Run this — the performance package’s all-in-one diagnostics:
check_model(model_aov)▶ Run this — Levene’s test (H₀: variances are homogeneous — you want a non-significant result):
levene_test <- leveneTest(phase_shift ~ treatment, data = c_df)
levene_test▶ Run this — Shapiro-Wilk on the residuals (H₀: normally distributed — you want a non-significant result):
shapiro_test <- shapiro.test(residuals(model_aov))
shapiro_test✏️ Your turn: Do these three groups pass both the normality and equal-variance checks? ________________________
🚀 If you finish early: In the base R diagnostic plots, which panel would you check first for equal variance, and which for normality? Name the two panels.
________________________________________________________Part 3 · Estimated Marginal Means and Post-Hoc Testing
When ANOVA rejects H₀, we need to determine which groups differ. First, we calculate estimated marginal means (EMMs).
▶ Run this:
emmeans_df <- emmeans(model_aov, "treatment")
emmeans_dfEstimated Marginal Means (EMMs), also called least-squares means, are model-based predictions of group means from your fitted model rather than directly from raw data. In a one-way ANOVA with balanced data and no covariates, EMMs = group means exactly — but the emmeans package’s approach generalizes cleanly to two-way ANOVA, ANCOVA, and more complex designs, and gives you built-in standard errors, confidence intervals, and comparisons for free.
▶ Run this — verify EMMs equal the simple group means:
group_means <- c_df %>%
group_by(treatment) %>%
summarize(
sample_mean = mean(phase_shift),
sample_sd = sd(phase_shift),
n = n()
)
emm_df <- as.data.frame(emmeans_df)
comparison <- group_means %>%
left_join(emm_df %>% select(treatment, emmean, SE),
by = "treatment") %>%
mutate(
difference = round(sample_mean - emmean, 10) # Round to show they're identical
)
comparison✏️ Your turn: Is difference zero (or extremely close to it, within rounding error) for every row? ________________________
▶ Run this — pairwise comparisons with a Sidak adjustment:
pairwise_comparisons <- pairs(emmeans_df, adjust = "sidak")
pairwise_comparisonsThe estimate is the difference between two group means. The SE is the standard error of that difference — it depends on the within-group variability and the sample sizes being compared. The t-ratio = estimate / SE tells you how many standard errors away from zero the difference is.
▶ Run this — letters showing which groups are (not) significantly different:
letter_groups <- multcomp::cld(emmeans_df, Letters = letters, adjust = "sidak")
letter_groups▶ Run this — the same thing in one pipeline:
cld_result <- emmeans(model_aov, "treatment") %>%
multcomp::cld(Letters = letters, adjust = "sidak")
cld_result▶ Run this — visualize the EMMs:
plot(emmeans_df)
emmip(model_aov, ~ treatment, CIs = TRUE) +
theme_minimal() +
labs(title = "Estimated Marginal Means with 95% CIs",
x = "Treatment", y = "Estimated Phase Shift")✏️ Your turn: Which two treatments share a compact-letter-display letter (meaning they’re not significantly different)? ________________________
🚀 If you finish early: Rerun pairs(emmeans_df, adjust = ...) with "bonferroni" and "tukey" instead of "sidak". Do the p-values change much for this dataset?
# Write your code here:Part 4 · Planned Comparisons
Instead of comparing every possible pair, sometimes you have specific, planned comparisons in mind.
▶ Run this:
levels(c_df$treatment) # See the order
emm <- emmeans(model_aov, "treatment")
planned_contrasts <- contrast(emm,
method = list(
"control vs eyes" = c(1, -1, 0),
"control vs knees" = c(1, 0, -1)))
planned_contrasts▶ Run this — the significance-groups plot:
ggplot(c_df, aes(x = treatment, y = phase_shift, color = treatment)) +
geom_jitter(width = 0.2, alpha = 0.7) +
stat_summary(fun = mean, geom = "point", size = 4, shape = 18) +
stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", width = 0.2) +
geom_text(data = as.data.frame(cld_result),
aes(x = treatment, y = -2.5, label = .group),
vjust = 0.5, size = 5)✏️ Your turn: The contrast vector c(1, -1, 0) picks out “Control minus Eyes.” Write the contrast vector for “Knees minus Eyes” instead (hint: check the order from levels(c_df$treatment)).
________________________________________________________Part 5 · Reporting the Results
Formal scientific writing example:
“The effect of light treatment on circadian rhythm phase shift was analyzed using a one-way ANOVA. There was a significant effect of treatment on phase shift (F(2, 19) = 7.29, p = 0.004, η² = 0.43). Post-hoc comparisons using Tukey’s HSD test indicated that the mean phase shift for the Eyes treatment (M = -1.55 h, SD = 0.71) was significantly different from both the Control treatment (M = -0.31 h, SD = 0.62) and the Knees treatment (M = -0.34 h, SD = 0.79). However, the Control and Knees treatments did not significantly differ from each other.”
✏️ Your turn: Using your own output from Parts 1–3, write a results sentence for the Sidak-adjusted pairwise comparisons instead of Tukey’s.
________________________________________________________Part 6 · Non-Parametric ANOVA
▶ Run this — build a version of the data that clearly violates ANOVA’s assumptions:
set.seed(42)
v_circ_df <- c_df %>%
mutate(
phase_shift_violated = case_when(
# Control: very tight, almost no variance
treatment == "Control" ~ phase_shift * 0.15,
# Knees: also tight variance
treatment == "Knees" ~ phase_shift * 1.3,
# Eyes: extreme spread - compress middle values, huge outliers
treatment == "Eyes" ~ {
n <- length(phase_shift)
c(phase_shift[1:(n-3)] * 1.2, # Very compressed normal values
phase_shift[(n-2)] * 2, # Moderate outlier
phase_shift[(n-1)] * 2.4, # Extreme outlier 1
phase_shift[n] * 4) # Extreme outlier 2
}
)
)
violated_model <- lm(phase_shift_violated ~ treatment,
data = v_circ_df)▶ Run this — visualize and check assumptions on the violated data:
ggplot(v_circ_df,
aes(x = treatment, y = phase_shift_violated,
color = treatment)) +
geom_jitter(width = 0.2, alpha = 0.7, size = 3) +
geom_boxplot(alpha = 0.3, outlier.shape = NA) +
geom_jitter(width = 0.2, alpha = 0.7, size = 3) +
theme_minimal() +
labs(title = "Modified Data with Violations",
subtitle = "Median with IQR",
x = "Light Treatment", y = "Phase Shift (hours)")
shapiro.test(resid(violated_model))
qqnorm(resid(violated_model), main = "Q-Q Plot: Violated Data")
qqline(resid(violated_model), col = "red", lwd = 2)
leveneTest(phase_shift_violated ~ treatment,
data = v_circ_df)✏️ Your turn: Do the Shapiro-Wilk and Levene’s tests confirm this data violates ANOVA’s assumptions? ________________________
▶ Run this — the Kruskal-Wallis test, on both the original and violated data:
kruskal_original_result <- kruskal.test(
phase_shift ~ treatment,
data = c_df
)
kruskal_original_result
anova(model_aov) # for comparison
kruskal_violated_result <- kruskal.test(
phase_shift_violated ~ treatment,
data = v_circ_df
)
kruskal_violated_result
anova(violated_model) # for comparison✏️ Your turn: Does the parametric ANOVA’s conclusion on the violated data agree with the Kruskal-Wallis conclusion? If they disagree, which would you trust more, and why? ________________________
▶ Run this — post-hoc tests: pairwise Wilcoxon, and the more appropriate Dunn’s test:
pairwise_result <- pairwise.wilcox.test(
x = v_circ_df$phase_shift_violated,
g = v_circ_df$treatment,
p.adjust.method = "bonferroni"
)
pairwise_result
dunn_result <- dunnTest(
phase_shift_violated ~ treatment,
data = v_circ_df,
method = "bonferroni"
)
dunn_resultWhich post-hoc test should you use after Kruskal-Wallis?
Use Dunn’s test — it’s the proper follow-up to Kruskal-Wallis, uses the same overall ranking from the omnibus test, and is standard in published research. Pairwise Wilcoxon re-ranks separately for each pair, throwing away that shared-ranking information.
🚀 If you finish early: Run Dunn’s test with method = "holm" instead of "bonferroni". Do any conclusions change? Which method is more conservative?
# Write your code here:Review and checkpoint
At this point you can:
📤 What to turn in before next class
Upload both of these to the course management system:
- Your code — the
scripts/folder (or just12_analysis_of_variance.R) - This worksheet, with your written answers
Part 7 · Take-Home Extension — Crayfish Growth Across Lakes
Due Monday, October 26 — before the Two-Way ANOVA class.
Same workflow as Parts 1–3 — fit, check assumptions, compare groups — new dataset, new question. This time you decide whether the standard ANOVA route or the Kruskal-Wallis route is the right call, and you pick the post-hoc test yourself; nobody walks you through it step by step.
Background
Sargent & Lodge (2014) reared young-of-year rusty crayfish (Orconectes rusticus) in enclosures across three northern Wisconsin lakes. Below are their daily growth rate measurements.
▶ Run this:
cray_df <- read_csv("data/sargent_lodge_crayfish.csv")
head(cray_df)Your Task
Question: Does crayfish growth rate (growth_per_day) differ among the three lakes (lake)?
✏️ Your turn: State your hypotheses. Check the assumptions yourself and decide whether a standard one-way ANOVA is appropriate here, or whether you should reach for Kruskal-Wallis instead. If groups differ, pick and justify a post-hoc test. Justify your choices, run the analysis, interpret it, and report the result properly.
H0 =
________________________________________________________
Ha =
________________________________________________________
Test I will use and why:
________________________________________________________# Write your code here:Interpretation:
________________________________________________________
________________________________________________________Final Figure
Produce one publication-quality figure comparing growth rate across lakes, with post-hoc groupings shown (e.g. a compact letter display) — proper axis labels, no default ggplot grey background — and export it with ggsave().
# Write your code here:📤 What to turn in — due Oct 26
Your write-up should include:
Reference
Sargent, L.W. & Lodge, D.M. (2014). Evolution of invasive traits in nonindigenous species: increased survival and faster growth in invasive populations of rusty crayfish (Orconectes rusticus). Evolutionary Applications, 7, 949–961.
Getting unstuck
When code breaks — and it will, that is normal:
- Read the error message out loud. R usually names the line and the problem.
- Check the usual suspects: did you run
library(performance),library(FSA),library(car),library(emmeans), andlibrary(tidyverse)? Spelling? A missing)or%>%at the start of a line? ?function_nameopens the built-in help page.- Bring the exact error (copy-paste it) to class or office hours.
💡 Key idea: Every working scientist googles error messages daily. Getting stuck is not failing — it is the job.