Activity: Nested ANOVA

Nested designs

Hands-on activity: nested ANOVA on the sea urchin grazing / algae cover dataset. Fit it the wrong way first, then with patch as the error term, and prove the nested F-test is really a test on the patch means.

Author

Bill Perry

In-class Activity 14: Nested ANOVA

Today’s Objectives

Tell a nested design from a factorial one by looking at the data
Count how many times a treatment was actually applied, not how many rows you have
Fit the wrong model first — two different ways of getting it wrong
Fit the right one, with patch as the error term
Cross-check it with afex::aov_car()
Run post-hoc comparisons on the correct model
Prove the nested F-test is really a one-way ANOVA on the patch means — and watch that break when the design goes unbalanced
Report an effect size with effectsize, and tell a null result from an underpowered one
Name what goes wrong when you forget the design, and report the result honestly

📂 Files for today

📥 Urchin data 📥 Script skeleton (.R)

Save the skeleton as scripts/14_nested_anova.R, and andrew.csv in your data/ folder.

What you hand in

All of today’s code is right here in this activity. Type it into your script as we go and run it chunk by chunk. The skeleton has the plumbing — libraries, the data read, the treatment recode, the summarise block — and nothing else. Every line of statistics you type yourself.
The Q boxes are things to think about, not inline homework. Read them, talk them over, answer them in your head. You do not have to type answers into the script, and they are not graded.
The write-up at the bottom — sections 8 and 9 — is graded. Those answers go into the script as # comments.
Hand in scripts/14_nested_anova.R. Run it top to bottom (Ctrl/Cmd + Shift + Enter) first — it must run with no errors.

🔮 Predict before you run

Today has a trap in it, on purpose. Section 3 gives you an answer that looks great and is wrong. Before you run each model, say out loud what you expect the p-value to do — then run it. The gap between your prediction and the output is the lesson.

The invisible(NULL) lines in the skeleton

Every section ends with invisible(NULL). It is a bumper. Press Ctrl/Cmd + Enter on a section’s last real line and your cursor jumps to the next line that has code on it — without the bumper, that is the next section, way down past the next heading. With it, the cursor stops right there and you keep your place. It does nothing and prints nothing.

Keep the plots plain today — no labels, colors, or themes. The goal is to read them. The one exception is the final figure in section 9, which is the graded deliverable.

🧩 Block 1 — What makes a design nested (after lecture Parts 1–3)

Sections 1–2: read the data, then see the nesting in it.

Section 1 · The data

Written for you in the skeleton — run it and read it.

# ---- Libraries -------------------------------------------
library(janitor)
library(afex)       # aov_car()
library(emmeans)
library(effectsize) # eta_squared(), omega_squared()
library(tidyverse)

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

Andrew’s sea urchin grazing experiment. Sixteen patches of reef each got one urchin-density treatment; five quadrats were then scored for algae cover inside every patch.

Column	What it is
`treat`	urchin density — 0%, 33%, 66%, 100% of natural
`patch`	reef patch 1–16, each with one treatment
`quad`	quadrat within a patch
`algae`	% cover of filamentous algae — the response

# ---- 1 The data ------------------------------------------
u_df <- read_csv("data/andrew.csv") %>%
  clean_names() %>%
  mutate(
    treat = case_when(
      treat == "removal" ~ "0%",     # book   0% density
      treat == "dens_66" ~ "33%",    # book  33% density
      treat == "dens_33" ~ "66%",    # book  66% density
      treat == "control" ~ "100%"    # book 100% density
    ),
    treat = factor(treat, levels = c("0%", "33%", "66%", "100%")),
    patch = as_factor(patch)
  )

head(u_df)

u_df %>%
  group_by(treat) %>%
  summarise(
    n_quads = n(),
    n_patches = n_distinct(patch),
    mean_algae = mean(algae),
    sd_algae = sd(algae),
    .groups = "drop"
  )

⚠️ Watch out! The file’s dens_33 / dens_66 labels are the reverse of the book’s percentages. That is verified by matching patch means against Quinn & Keough Table 9.1 — dens_33 is the book’s 66% group and dens_66 its 33%. The recode above fixes it, which is why it is written for you.

✏️ Q1: You have 80 rows. How many times was the treatment actually applied? Those two numbers are not the same — what is each one?

✏️ Q2: Read the mean_algae column. The treatment means run from about 1.3% to 39.2%. Eyeballing it, does urchin density look like it matters?

Section 2 · See the nesting

🔮 Predict first: if you tabulate patch against treat, how many of the 64 cells will be non-zero?

▶ Type this yourself:

# ---- 2 See the nesting -----------------------------------
u_df %>%
  count(treat, patch) %>%
  pivot_wider(names_from = treat, values_from = n, values_fill = 0)

▶ Type this yourself:

u_df %>%
  ggplot(aes(x = treat, y = algae, color = patch)) +
  geom_point(position = position_dodge(width = 0.5))

✏️ Q3: Does patch 3 ever appear under the 0% treatment? Does any patch appear under more than one treatment?

✏️ Q4: In a factorial design every level of one factor appears with every level of the other. By that test, is this factorial or nested? Write the one-sentence rule you would use to tell them apart.

✏️ Q5: Look at the spread of points within a single treatment — pick 0%. Are the four patches in that treatment similar to each other, or wildly different?

💡 Key idea: patch is not a second factor you are interested in. It is the level at which the treatment was applied — the experimental unit. The quadrats inside a patch are subsamples of that unit, not independent replicates of the treatment.

🧩 Block 2 — The wrong model, then the right one (after lecture Parts 4–6, 8)

Sections 3–4: fit it ignoring the nesting, then fix it.

Section 3 · The wrong way — two of them

There are two ways to get this wrong, and you are going to make both on purpose so you can recognize them later.

🔮 Predict first: with 80 observations and treatment means from 1.3% to 39.2%, will either of these come out significant?

▶ Type this yourself — the naive factorial. patch is in the model, but treatment gets tested against the leftover quadrat-to-quadrat variation:

# ---- 3 The wrong way - two of them -----------------------
naive_model <- lm(algae ~ treat + patch, data = u_df)
anova(naive_model)

▶ Type this yourself — the cruder mistake: leave patch out altogether, so all 80 quadrats are treated as independent replicates of the treatment:

ignore_model <- aov(algae ~ treat, data = u_df)
summary(ignore_model)

✏️ Q6: Write down F, both df, and p for treat from both models. Significant at α = 0.05? By how much?

✏️ Q7: The two wrong answers are not even the same wrong answer. Look at the residual mean squares — 298.6 against 531. Why does leaving patch out make the denominator bigger, and so the F smaller?

✏️ Q8: anova(naive_model) also gives you a row for patch: F = 5.93, p < 0.001. Patches genuinely differ from one another. Hold onto that — it is the thing section 4 is about to make you pay for.

✏️ Q9: Where did the residual df come from in each model? Count them out from 80 observations, 4 treatments, 16 patches. What is each model claiming your sample size is?

Section 4 · The right way — patch as the error term

The treatment was applied to patches, so the treatment effect has to be tested against patch-to-patch variation, not quadrat-to-quadrat variation. Error(patch) tells aov() exactly that.

🔮 Predict first: the data has not changed and the treatment means have not changed. Will the p-value move a little, or a lot? Which direction?

▶ Type this yourself:

# ---- 4 The right way - patch as the error term -----------
nested_model <- aov(algae ~ treat + Error(patch), data = u_df)
summary(nested_model)

⚠️ Watch out! Error(patch) is enough here — patch IDs are unique 1–16 and each patch got exactly one treatment. Writing Error(treat:patch) gives identical numbers but warns Error() model is singular.

▶ Type this yourself — cross-check with afex, which is the better tool for unbalanced nested designs:

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

aov_car() prints a note that it is aggregating to one value per design cell — that is expected, and it is doing the same thing section 6 will make you do by hand.

✏️ Q10: Write down the F, both df, and the p for the nested model. Compare to Q6. Was your prediction right?

✏️ Q11: Now compare the treatment sums of squares in section 3 and section 4. They are the same number. So what actually changed between the two models?

✏️ Q12: The residual df for treat was 64 in the naive factorial and 76 with patch ignored entirely — and it is 12 in the nested model. Explain where the 12 comes from. (16 patches, 4 treatments.) Which of the three is your real sample size telling you?

✏️ Q13: Do nested_model and model_afx agree? They should — if they do not, one of them has a typo.

💡 Key idea: the numerator of the F-ratio did not move at all. The whole change is in the denominator — what you are willing to call noise. Treating subsamples as replicates shrinks the denominator, inflates F, and manufactures significance out of nothing. That is pseudoreplication, and it is the most common fatal flaw in ecological papers.

🧩 Block 3 — Post-hoc, and where the variance lives (after lecture Part 9)

Sections 5–7: compare the treatments, reproduce what the nested test is really doing, then put an effect size next to the p-value.

Section 5 · Post-hoc comparisons

▶ Type this yourself:

# ---- 5 Post-hoc comparisons ------------------------------
emm <- emmeans(nested_model, ~ treat)
emm

pairs(emm, adjust = "sidak")

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

⚠️ Watch out! emmeans() warns EMMs are biased unless design is perfectly balanced. Here the design is balanced — 4 patches per treatment, 5 quadrats per patch — so you can read past it. It fires on any aovlist object.

✏️ Q14: The treatment means run from 1.3% to 39.2% — a 30-fold spread. How many letters does the compact letter display use?

✏️ Q15: How can means that far apart fail to be distinguishable? Look at the SE column — what is the standard error on each mean, and where did it come from?

✏️ Q16: Read the last line of the cld() output carefully: “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.” Why is that second sentence there?

Section 6 · Where the variance actually is (you write this one)

This is the payoff. The lecture claimed the nested F-test is not really a test on 80 quadrats at all — it is a one-way ANOVA on the 16 patch means, wearing a disguise. Do it yourself and see how exactly it comes out.

🔮 Predict first: collapse the data to one mean per patch — 16 rows — and run an ordinary one-way ANOVA on those. What F and p will you get?

▶ Type this yourself:

# ---- 6 Where the variance actually is --------------------
patch_df <- u_df %>%
  group_by(treat, patch) %>%
  summarise(mean_algae = mean(algae), .groups = "drop")

patch_df

▶ Type this yourself:

patch_model <- aov(mean_algae ~ treat, data = patch_df)
summary(patch_model)

✏️ Q17: Compare this F, df, and p to section 4. Exactly the same, or just close? Was your prediction right?

✏️ Q18: The sums of squares are not the same as section 4 — they are exactly 1/5 of them. Why 1/5?

✏️ Q19: So what did the five quadrats per patch actually buy you? (They did not buy you df for the treatment test — look at Q12.)

Written for you in the skeleton — now break the balance. A storm scours all but one quadrat out of four of the patches. The filter() and the re-summarise are bookkeeping; the two aov() calls under them are the point. Run it and read it.

thin_df <- u_df %>%
  filter(!(patch %in% c(2, 6, 10, 14) & quad > 1))

thin_df %>% count(patch)

# the nested ANOVA on the thinned data
summary(aov(algae ~ treat + Error(patch), data = thin_df))

# the patch-mean shortcut, on the same thinned data
thin_patch_df <- thin_df %>%
  group_by(treat, patch) %>%
  summarise(mean_algae = mean(algae), .groups = "drop")

summary(aov(mean_algae ~ treat, data = thin_patch_df))

✏️ Q19b: On the balanced data those two gave you the same F and p. Do they still? Write both down. Which one is more optimistic?

✏️ Q19c: A patch mean built from 1 quadrat is noisier than one built from 5 — its variance is \(\sigma^2_{patch} + \sigma^2_{within}/n_j\). The one-way ANOVA on the patch means does not know that, and weights all 16 equally. What assumption of that ANOVA is now violated by construction?

✏️ Q19d: To weight each patch by how precise it is, you would need \(\sigma^2_{within}\). Where did it go when you averaged? This is the real reason averaging is not a free shortcut.

✏️ Q20: Print patch_df and look at the four patch means inside the 0% treatment: roughly 34.2, 62.0, 2.2, 58.4. Would running more quadrats per patch have made this experiment better able to detect a treatment effect? What would have?

💡 Key idea: subsamples make each experimental unit’s value more precise. They do not make you more replicated. If you want power to detect a treatment effect, you need more patches — more places where the treatment was independently applied. Twenty quadrats in four patches is still n = 4.

But the shortcut is only free when the design is balanced. Once patches hold different numbers of quadrats, their means stop being equally precise, the averaging trick stops matching the nested test, and you have thrown away the one quantity — \(\sigma^2_{within}\) — you would need to weight them properly. That costs you real power, and it is exactly the problem a mixed model is built to solve.

Section 7 · Effect size — a null result, or an underpowered one?

Section 4 gave you p = 0.091 and you wrote it down as “not significant.” That is a statement about confidence, not about size. Before you report a non-significant result, you owe the reader the other number: how big was the effect you failed to resolve?

effectsize takes the aovlist from aov(..., Error(patch)) directly — and the Group column in the output tells you which error stratum it used.

▶ Type this yourself:

# ---- 7 Effect size ---------------------------------------
eta_squared(nested_model)

omega_squared(nested_model)

🔮 Predict first: the naive factorial in section 3 had p = 0.00000007 and the nested model p = 0.091 — seven orders of magnitude apart. Will their effect sizes be that far apart too?

▶ Type this yourself — the same effect size off the two wrong models from section 3:

eta_squared(naive_model)

eta_squared(ignore_model)

You should get roughly:

model                            partial eta^2      p
nested, Error(patch)   (right)       0.40         0.091
naive factorial        (wrong)       0.43         0.00000007
patch ignored entirely (wrong)       0.26         < 0.001

▶ Type this yourself — turn it into a word:

interpret_eta_squared(eta_squared(nested_model)$Eta2_partial,
                      rules = "field2013")

✏️ Q21: Write down partial η² for the nested model and for the naive factorial. Was your prediction right? Which moved more between the two models — the effect size or the p-value?

✏️ Q22: interpret_eta_squared() calls the nested effect large, and its p-value is 0.091. Those are not in conflict. Say in one sentence what each one is measuring.

✏️ Q23: Look at the confidence interval on the nested η². It runs from 0 to 1 — it includes zero, which is the same message as p > 0.05. Now look at the CI for the naive model, which does not include zero. Which of those two CIs is lying to you, and why?

💡 Key idea: pseudoreplication corrupts the p-value, not the effect size. The wrong model and the right model agree the treatment effect is large — 0.43 against 0.40. What they disagree about, by seven orders of magnitude, is how confident you get to be about it.

That changes what you are entitled to write in section 8. “Urchin density had no effect on algal cover” is not what these data say. “A large effect (partial η² = 0.40) that four patches per treatment could not resolve (F₃,₁₂ = 2.72, p = 0.091)” is. The first sentence blames the urchins; the second blames the design — and only the second one is honest.

Write-up — graded

These go into the bottom of your script, as # comments after ANSWER:.

Section 8 · Report the model

✏️ Q24: Fill in all three models:

                              F      df        p
naive factorial (treat+patch)
patch ignored entirely
nested, Error(patch)

✏️ Q25: Write a 2–3 sentence results paragraph for the correct model, using your own F, df, and p.

✏️ Q25b: Add the effect size to your Q25 paragraph — partial η², its CI, and one clause saying what it means. Your paragraph should now make clear this was an underpowered experiment, not a demonstration that urchins do not matter.

✏️ Q26: If you had reported either of the section 3 numbers instead, what would you have claimed? What is the name for that mistake?

✏️ Q27: Does this experiment show urchins have no effect on algae cover? Be careful here — separate “we did not detect an effect” from “there is no effect,” and say which one you are entitled to.

Section 9 · The final figure

One publication-quality figure that shows the nesting: treatment on the x-axis, the patch means as points, so a reader can see the among-patch spread the F-test is actually working with. Proper axis labels, a real theme, readable sizing — this one is graded on presentation. Export it with ggsave().

# ---- 9 The final figure ----------------------------------
# Build it from patch_df, then:
# ggsave("figures/14_urchin_patches.pdf", width = 7, height = 5)

✏️ Q28: Your figure shows why the test came out the way it did. Point at the specific feature of the figure that explains the non-significant p-value.

Other Formats

Activity: Nested ANOVA

In-class Activity 14: Nested ANOVA

Today’s Objectives

🧩 Block 1 — What makes a design nested (after lecture Parts 1–3)

Section 1 · The data

Section 2 · See the nesting

🧩 Block 2 — The wrong model, then the right one (after lecture Parts 4–6, 8)

Section 3 · The wrong way — two of them

Section 4 · The right way — patch as the error term

🧩 Block 3 — Post-hoc, and where the variance lives (after lecture Part 9)

Section 5 · Post-hoc comparisons

Section 6 · Where the variance actually is (you write this one)

Section 7 · Effect size — a null result, or an underpowered one?

Write-up — graded

Section 8 · Report the model

Section 9 · The final figure

Review and checkpoint