Activity: Linear Mixed Models

Random effects for nested data

Hands-on activity: refitting the sea urchin nested design with lme4::lmer(), choosing between the F-test and the Chi-square test when they disagree, reading variance components and the ICC, and proving the mixed model reproduces the traditional nested ANOVA.

Author

Bill Perry

In-class Activity 16: Linear Mixed Models

Today’s Objectives

Tell a fixed effect from a random effect, and say why patch is random
Write random-effects notation for nested and for crossed designs
Fit a nested design as a mixed model with lme4::lmer()
Run the F-test and the Chi-square test, and find that they disagree
Choose between them, and defend the choice
Read variance components and the ICC — information the nested ANOVA never gave you
Reproduce the patch variance from week 9’s mean squares, and say why REML does not need to
Run DHARMa diagnostics on a mixed model
Prove lmer() reproduces aov(... + Error()) on a balanced design

📂 Files for today

📥 Urchin data 📥 Script skeleton (.R)

Save the skeleton as scripts/16_linear_mixed_models.R, and andrew.csv in your data/ folder.
Same data as Activity 14 — if you still have that script, keep it open.

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 and the data read — and nothing else. Every line of statistics you type yourself.
Section 6 has no code here either. It is a direct comparison back to Activity 14, and you have already written every line of it once.
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 7 and 8 — is graded. Those answers go into the script as # comments.
Hand in scripts/16_linear_mixed_models.R. Run it top to bottom (Ctrl/Cmd + Shift + Enter) first — it must run with no errors.

🔮 Predict before you run

Today’s trap is not a broken assumption — it is two legitimate tests that give you different answers on either side of α = 0.05. Before you run each one, say which way you expect it to come out.

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 8, which is the graded deliverable.

Where this fits

In Activity 14 you analysed the sea urchin experiment with aov(algae ~ treat + Error(patch)). It gave the right answer — F = 2.72 on 3 and 12 df, p = 0.091 — by forcing the treatment to be tested against patch-to-patch variation.

Error() is the 1960s way to get the right denominator. It works, and it is brittle: it needs a balanced design, it cannot handle more than one grouping cleanly, and it tells you nothing about how much variation lives among patches.

A mixed model does the same job by naming patch as a random effect. Same question, better machinery — and one new piece of output that the nested ANOVA simply could not produce.

	Fixed effect	Random effect
What it is	levels you chose	levels drawn from a population
What you want	the difference between them	how much they vary
Here	`treat` — the four densities	`patch` — 16 bits of reef
Generalizes to	these four densities	reef patches in general

Random-effects notation:

Code	Meaning
`(1\|group)`	a random intercept for each group
`(slope\|group)`	random intercept and random slope
`(1\|outer/inner)`	`inner` nested inside `outer`
`(1\|a) + (1\|b)`	crossed, not nested

🧩 Block 1 — Fixed, random, and the model (after lecture Parts 1–2)

Sections 1–2: the data again, then the mixed model.

Section 1 · The sea urchin data, again

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

# ---- Libraries -------------------------------------------
library(janitor)
library(car)          # Anova()
library(lme4)         # lmer()
library(lmerTest)     # Satterthwaite df and p-values
library(emmeans)
library(performance)  # icc()
library(DHARMa)
library(tidyverse)

options(scipen = 999)
options(contrasts = c("contr.sum", "contr.poly"))  # for Anova(type = 3)

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

head(u_df)

✏️ Q1: In Activity 14 you got p = 0.091 and stopped there. Two weeks on, what question about those 16 patches is still completely unanswered?

Section 2 · Fit the mixed model

▶ Type this yourself:

# ---- 2 Fit the mixed model -------------------------------
mixed_model <- lmer(algae ~ treat + (1|patch), data = u_df,
                    control = lmerControl(calc.derivs = FALSE))

summary(mixed_model)

control = lmerControl(calc.derivs = FALSE) just suppresses a convergence check that is noisy on small datasets. It does not change the fit.

✏️ Q2: Find the Random effects block in the summary. There are two rows. What is each one, and what are the two numbers?

✏️ Q3: Why is patch the random effect and treat the fixed one? What would it mean to call patch fixed — what would you be claiming you wanted to know?

✏️ Q4: The summary reports Number of obs: 80, groups: patch, 16. Which of those two numbers is your real replication for testing treat?

🧩 Block 2 — Two tests that disagree (after lecture Part 3)

Section 3: the same model, three tests, two different answers.

Section 3 · Two tests, two answers

Mixed models have an awkward problem: there is no single agreed way to count the denominator degrees of freedom. Different approximations give different p-values.

🔮 Predict first: Activity 14’s nested ANOVA gave p = 0.091. Will the mixed model’s F-test land near that, or somewhere else?

▶ Type this yourself — the F-test, with Kenward-Roger df:

# ---- 3 Two tests, two answers ----------------------------
Anova(mixed_model, type = 3, test.statistic = "F")

💡 Read the header R prints: “Kenward-Roger df”. car::Anova() has no ddf argument — it always uses Kenward-Roger when you ask for test.statistic = "F".

▶ Type this yourself — the Satterthwaite version, which is what lmerTest adds:

anova(mixed_model)

🔮 Predict first: before the next chunk — the Chi-square test is the same model again. Same answer, or different?

▶ Type this yourself — the Chi-square test:

Anova(mixed_model, type = 3, test.statistic = "Chisq")

✏️ Q5: Write down all three p-values. Two agree and one does not. Which is the odd one out, and which side of 0.05 does it fall on?

✏️ Q6: Chi-square = F × numerator df, and the Chi-square test assumes infinite denominator df. How many do you actually have? Is 12 infinite?

✏️ Q7: The rule of thumb: under ~100 observations or fewer than ~20 random-effect levels, use the F-test; over ~500 observations and more than ~50 levels, Chi-square is fine; in between, the F-test is safer. You have 80 observations and 16 patches. Which does the rule pick?

✏️ Q8: Here is the uncomfortable part. If you had run only the Chi-square test, you would have reported a significant treatment effect and never known. What does that tell you about running one test and stopping?

✏️ Q9: Could you defend the F-test to a reviewer who wanted the Chi-square? Write the one sentence you would say.

💡 Key idea: this is not a case of one test being broken. Both are legitimate approximations; they differ because one of them pretends you have more information than you do. When two defensible analyses straddle α = 0.05, the honest move is to say so and justify your choice — not to pick the one you like and delete the other.

🧩 Block 3 — What the nested ANOVA couldn’t tell you (after lecture Parts 4–6)

Sections 4–5: variance components, then diagnostics.

Section 4 · Variance components — the new information

Everything so far has reproduced Activity 14. This is what Error(patch) could never give you.

🔮 Predict first: look back at Activity 14’s patch means — within the 0% treatment alone they ran from 2.2% to 62%. Before you run this: is more of the variation among patches, or among quadrats within a patch?

▶ Type this yourself:

# ---- 4 Variance components - the new information ---------
VarCorr(mixed_model)

icc(mixed_model)

✏️ Q10: Write down the patch variance and the residual variance. Which is bigger?

✏️ Q11: The ICC (intraclass correlation) is the share of total variance sitting among patches. What is yours? Say what it means in one plain sentence about this reef.

✏️ Q12: What would an ICC near zero mean about the patches? Near one?

✏️ Q13: The nested ANOVA in Activity 14 gave you a p-value and nothing else. Name the practical decision you could make with the ICC that you could not make with the p-value alone. (Think about designing the next experiment.)

▶ Type this yourself — where that patch variance came from, the old way:

a_tab <- summary(aov(algae ~ treat + treat:patch, data = u_df))[[1]]
a_tab

MS_patch  <- a_tab[2, "Mean Sq"]
MS_within <- a_tab[3, "Mean Sq"]

(MS_patch - MS_within) / 5

🔮 Predict first: lmer() never computed a sum of squares. That line above is pure arithmetic on week 9’s ANOVA table. Will it land on the same number VarCorr() gave you, or just near it?

✏️ Q13b: Did it match? To how many decimal places?

✏️ Q13c: That / 5 is the number of quadrats per patch. Suppose a storm had cut four patches down to one quadrat each. What number would you divide by then — and does lmer() have to answer that question at all?

✏️ Q13d: MS_patch could come out smaller than MS_within by chance, which makes that subtraction negative. A variance cannot be negative. What do you think lmer() reports instead, and what warning have you been told to expect when it does?

▶ Type this yourself — the treatment means from the model:

emm <- emmeans(mixed_model, ~ treat)
emm

pairs(emm, adjust = "sidak")

multcomp::cld(emm, Letters = letters)

✏️ Q14: Look hard at the confidence limits on the 33% and 100% means. They go below zero. Algae cover is a percentage — it cannot be negative.

✏️ Q15: You spent all of last week on models that refuse to predict impossible values. So what is this model assuming about algae that is not true, and what would you do about it? (You do not have to fix it today — name the fix.)

💡 Key idea: the mixed model fixed the error structure and left the distribution alone. Last week you fixed the distribution and left the error structure alone. Real data often needs both at once — that is a generalized linear mixed model, glmer(), and you now have both halves of it.

Section 5 · Diagnostics

▶ Type this yourself:

# ---- 5 Diagnostics ---------------------------------------
set.seed(123)
sim_mixed <- simulateResiduals(mixed_model, plot = FALSE)

plotQQunif(sim_mixed)
plotResiduals(sim_mixed)
testDispersion(sim_mixed)

▶ Type this yourself:

plot(mixed_model)

resid_df <- u_df %>%
  mutate(.resid = residuals(mixed_model))

ggplot(resid_df, aes(x = patch, y = .resid)) +
  geom_boxplot()

✏️ Q16: Do the DHARMa diagnostics look acceptable?

✏️ Q17: A mixed model also assumes the random effects are normally distributed. You have 16 patches. Is that enough to check that assumption honestly? What would you say in a paper?

Section 6 · Does it agree with Activity 14? (you write this one)

No code on this page. You wrote every line of this once already, two weeks ago.

The claim: on a balanced design, lmer(algae ~ treat + (1|patch)) and aov(algae ~ treat + Error(patch)) are the same test. Check it.

Work through:

Refit Activity 14’s nested model: aov() with + Error(patch) on u_df.
Put its F, df, and p next to the F-test from section 3.
Collapse to the 16 patch means and run a plain one-way aov() — as you did in Activity 14 — and add that row too.
Three rows, one table. Decide which numbers match and which do not.

✏️ Q18: Do all three agree? Exactly, or just closely?

✏️ Q19: If the F-test is identical to a method from the 1960s, what did you actually gain by fitting a mixed model? Name two things — section 4 is one of them.

✏️ Q20: This design is perfectly balanced: 4 patches per treatment, 5 quadrats per patch. Where would aov(Error()) and lmer() start to come apart? (What if one patch had been lost to a storm?)

💡 Key idea: the mixed model did not give you a better p-value — it gave you the same one, plus variance components, plus a path to unbalanced data, crossed random effects, and non-normal responses. That is why it replaced Error(), not because the old answer was wrong.

Homework

Homework 9 — Linear Mixed Models

The rusty crayfish growth experiment — but now lake is a random effect rather than a fixed factor. Same data as Homework 8, deliberately: the comparison between the two analyses is the assignment.

Data: sargent_lodge_crayfish.csv — Sargent & Lodge (2014), Evolutionary Applications 7: 949–961.

Start a new script for it — do not add it to 16_linear_mixed_models.R.

Write-up — graded

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

Section 7 · Report the model

✏️ Q21: Fill in the table:

method                        F / Chisq   df      p
lmer, F (Kenward-Roger)
lmer, F (Satterthwaite)
lmer, chi-square
aov with Error(patch)

✏️ Q22: Variance components: patch, residual, and the ICC.

✏️ Q23: A 2–3 sentence results paragraph. Report the fixed effect with the test you chose and the variance component — both belong in the sentence.

✏️ Q24: Does this experiment show urchins have no effect on algae cover? Same question as Activity 14 — has your answer changed now that you can see the variance components?

✏️ Q25: You now have three ways to analyse a nested design: aov with Error(), lmer with a random intercept, and a one-way ANOVA on the unit means. Which would you reach for first on a design of your own, and why?

Section 8 · The final figure

One publication-quality figure showing both levels of variation: the treatment means with their intervals from the model, and the individual patch means behind them. Proper axis labels, a real theme, readable sizing — this one is graded on presentation. Export it with ggsave().

# ---- 8 The final figure ----------------------------------
# ... your plot here, then:
# ggsave("figures/16_urchin_mixed.pdf", width = 7, height = 5)

✏️ Q26: Your figure shows two kinds of spread — among patches within a treatment, and among treatments. Which one is bigger, and which one does the p-value care about?

Other Formats

Activity: Linear Mixed Models

In-class Activity 16: Linear Mixed Models

Today’s Objectives

Where this fits

🧩 Block 1 — Fixed, random, and the model (after lecture Parts 1–2)

Section 1 · The sea urchin data, again

Section 2 · Fit the mixed model

🧩 Block 2 — Two tests that disagree (after lecture Part 3)

Section 3 · Two tests, two answers

🧩 Block 3 — What the nested ANOVA couldn’t tell you (after lecture Parts 4–6)

Section 4 · Variance components — the new information

Section 5 · Diagnostics

Section 6 · Does it agree with Activity 14? (you write this one)

Homework

Write-up — graded

Section 7 · Report the model

Section 8 · The final figure

Review and checkpoint