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

Part 1 · Setup and Data

▶ Run this in your Script:

# ── Part 1.1 · Setup and Data ───────────────────────────────────────────
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 %>%
  # Map to Quinn & Keough's Table 9.1 densities. NOTE the file's dens_33 /
  # dens_66 labels are the reverse of the book's percentages -- verified by
  # matching patch means: dens_33 is the book's 66% group, dens_66 its 33%.
  mutate(treat = case_when(
         treat == "removal" ~ "0%",     # book  0% density, mean algae 39.2
         treat == "dens_66" ~ "33%",    # book 33% density, mean algae 19.0
         treat == "dens_33" ~ "66%",    # book 66% density, mean algae 21.6
         treat == "control" ~ "100%",   # book 100% density, mean algae 1.3
         TRUE ~ "other")) %>%
  mutate(treat = factor(treat,
                        levels = c("0%", "33%", "66%", "100%"),
                        labels = c("0%", "33%", "66%", "100%")),
         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 — 0%, 33%, 66%, 100% of natural density), 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:

# ── Part 1.2 · Setup and Data ───────────────────────────────────────────
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:

# ── Part 2.1 · Fitting the Nested Model ─────────────────────────────────
# This gives you the correct F-test using PATCH as the error term.
# Patch IDs are unique (1-16) and each patch got exactly one treatment,
# so Error(patch) is enough -- Error(treat:patch) gives the same numbers
# but warns "Error() model is singular".
nested_model <- aov(algae ~ treat + Error(patch), data = u_df)
summary(nested_model)

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

# ── Part 2.2 · Fitting the Nested Model ─────────────────────────────────
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:

# ── Part 3.1 · Post-Hoc Comparisons ─────────────────────────────────────
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:

# ── Part 4.1 · Visualization and Reporting ──────────────────────────────
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:

# ── Part 4.2 · Visualization and Reporting ──────────────────────────────
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