Activity: Nested ANOVA

Nested designs

Hands-on activity: nested ANOVA the traditional way (aov with an Error term, afex) and post-hoc comparisons on the sea urchin grazing / algae cover dataset.
Author

Bill Perry

Worksheet: Nested 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.

This worksheet uses the sea urchin grazing / algae cover dataset from lecture, where patch is nested within treat (each patch received only one urchin-density treatment). Next class, the Linear Mixed Models activity revisits this exact same dataset with a different (and more flexible) modeling approach — this worksheet stays entirely with the traditional method.

Part 1 · Setup and Data

▶ Run this in your Script:

library(janitor)
library(afex)
library(emmeans)       # For estimated marginal means
library(tidyverse)

options(scipen = 999)
options(contrasts = c("contr.treatment", "contr.poly"))

u_df <- read_csv("data/andrew.csv") %>% clean_names()

u_df <- u_df %>%
  mutate(treat = case_when(
         treat == "removal" ~ "none",
         treat == "control" ~ "high",
         treat == "dens_33" ~ "low",
         treat == "dens_66" ~ "medium",
         TRUE ~ "other")) %>%
  mutate(treat = factor(treat,
                        levels = c("none", "low", "medium", "high"),
                        labels = c("none", "low", "medium", "high")),
         patch = as_factor(patch))

head(u_df)
levels(u_df$treat)

u_df contains: patch (random patches 1–16 where treatments were applied), treat (urchin density treatment — none, low, medium, high), quad (replicate quadrats within each patch:treatment combination), algae (percentage cover of filamentous algae — the response variable).

▶ Run this — summary statistics and a first plot:

summary_stats <- u_df %>%
  group_by(treat) %>%
  summarise(
    n = sum(!is.na(algae)),
    mean = mean(algae, na.rm = TRUE),
    sd = sd(algae, na.rm = TRUE),
    se = sd / sqrt(n),
    min = min(algae, na.rm = TRUE),
    max = max(algae, na.rm = TRUE),
    .groups = 'drop'
  )
summary_stats

dodge_position <- position_dodge(width = 0.3)

u_df %>%
  ggplot(aes(treat, algae, color = patch)) +
  stat_summary(fun = "mean", geom = "point", position = dodge_position) +
  stat_summary(fun.data = "mean_se", geom = "errorbar", width = 0.2, position = dodge_position)

✏️ Your turn: How many patches are nested within each treatment? ________________________

Part 2 · Fitting the Nested Model — Base R and afex

Since patch is nested within treat, we need the correct error term — the traditional aov() approach lets us specify this directly.

▶ Run this:

# This gives you the correct F-test using PATCH within TREAT as error term
nested_model <- aov(algae ~ treat + Error(treat:patch), data = u_df)
summary(nested_model)

▶ Run this — the afex package, recommended for unbalanced designs:

model_afx <- aov_car(algae ~ treat + Error(patch),
                     data = u_df)
summary(model_afx)

✏️ Your turn: Do nested_model and model_afx agree on whether the treatment effect is significant? ________________________

Part 3 · Post-Hoc Comparisons

▶ Run this:

emm <- emmeans(nested_model, ~ treat)
summary(emm)

pairs_result <- pairs(emm, adjust = "sidak")
pairs_result

cld_result <- multcomp::cld(emm, alpha = 0.05, Letters = letters)
cld_result

✏️ Your turn: Mean algae cover for Control (1.30%) looks much lower than the reduced-density treatments (66%: 21.55%, 33%: 19.00%, Removed: 39.20%). Based on your compact letter display, are any of those differences actually statistically significant at α = 0.05? ________________________

Part 4 · Visualization and Reporting

▶ Run this — a publication-style boxplot:

ggplot(u_df, aes(x = treat, y = algae, fill = treat)) +
  geom_boxplot(alpha = 0.7, outlier.shape = NA) +
  geom_jitter(width = 0.2, alpha = 0.4, size = 1) +
  scale_fill_viridis_d(option = "D", end = 0.85) +
  labs(
    title = "Effect of Urchin Density on Filamentous Algae Cover",
    x = "Urchin Density Treatment",
    y = "Filamentous Algae Cover (%)",
    caption = "Boxplots showing the distribution of algal cover across urchin density treatments."
  ) +
  theme(
    legend.position = "none",
    plot.title = element_text(face = "bold", size = 14),
    axis.title = element_text(face = "bold", size = 12),
    plot.caption = element_text(hjust = 0, face = "italic", size = 10)
  )

▶ Run this — a mean ± SE plot:

ggplot(summary_stats, aes(x = treat, y = mean, group = 1)) +
  geom_point(size = 3, shape = 21, fill = "white") +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.2) +
  labs(
    title = "Mean Algae Cover by Urchin Density Treatment",
    x = "Urchin Density Treatment",
    y = "Mean Filamentous Algae Cover (%)",
    caption = "Mean (+/- SE) percentage cover of filamentous algae across urchin density treatments."
  ) +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    axis.title = element_text(face = "bold", size = 12),
    plot.caption = element_text(hjust = 0, face = "italic", size = 10)
  )
Note

Scientific interpretation

This nested ANOVA reveals substantial spatial heterogeneity in algae cover, with significant variation among patches within each treatment. The substantial patch-to-patch variability shown in Part 2’s ANOVA table is exactly why patch needed to be modeled explicitly, rather than being pooled into the residual — a good reminder that pseudoreplication doesn’t just bias your F-test, it hides real biological variation.

✏️ Your turn: Write a 2–3 sentence results paragraph for this analysis, using your own F-statistic, df, and p-value from Part 2.

________________________________________________________
________________________________________________________

Review and checkpoint

At this point you can:

Note

📤 What to turn in before next class

Upload both of these to the course management system:

  1. Your code — the scripts/ folder (or just 14_nested_anova.R)
  2. This worksheet, with your written answers

Getting unstuck

When code breaks — and it will, that is normal:

  1. Read the error message out loud. R usually names the line and the problem.
  2. Check the usual suspects: did you run library(janitor), library(afex), library(emmeans), and library(tidyverse)? Spelling? A missing ) or %>% at the start of a line?
  3. ?function_name opens the built-in help page.
  4. 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.