Activity: Multiple Regression
Several predictors
Worksheet: Multiple Regression
📥 Bird data 📥 Tools script (.R) 📥 Script skeleton (.R)
The lecture’s dataset: 56 remnant forest patches in Victoria, Australia, from Loyn (1987). Bird abundance is the response, and six patch characteristics are the candidate predictors.
The order matters today. You will fit the predictors one at a time, then build one clean two-predictor model and read it carefully — that is the worked example. Only then do you put all six in at once and watch what happens.
Part 1 · Setup and the data
▶ Run this:
# ── Part 1.1 · Setup and the data ───────────────────────────────────────
library(tidyverse)
library(car)
source("themes_functions/r_themes_and_functions.R")
bird_df <- read_csv("data/loyn_birds.csv") %>%
mutate(log_area = log10(area_ha),
log_dist = log10(dist_m),
log_dist_larger = log10(dist_larger_m))
bird_df✏️ Your turn: Why are area and the two distances logged? Look at the range of area_ha.
____________________________________________________________Part 2 · One predictor at a time
▶ Run this — three separate simple regressions:
# ── Part 2.1 · One predictor at a time ──────────────────────────────────
summary(lm(abundance ~ log_area, data = bird_df))
summary(lm(abundance ~ grazing, data = bird_df))
summary(lm(abundance ~ year_isolated, data = bird_df))✏️ Your turn: Record each one.
log_area: slope = __________ p = __________
grazing: slope = __________ p = __________
year_isolated: slope = __________ p = __________
How many are significant? __________Part 3 · Two predictors — the worked example
Area and grazing were the two strongest predictors on their own. Put both in one model. This is the whole idea of multiple regression in its simplest form: one response, two predictors, two slopes.
▶ Run this:
# ── Part 3.1 · Two predictors — the worked example ──────────────────────
reduced_model <- lm(abundance ~ log_area + grazing, data = bird_df)
summary(reduced_model)✏️ Your turn: Write the fitted equation out, with the numbers in it.
abundance = ______ + ______ × log_area + ______ × grazing✏️ Your turn: Read each slope as a sentence. The phrase you need is “holding the other predictor constant.”
log_area: ___________________________________________________
grazing: ___________________________________________________Notice that both slopes shrank compared with Part 2 — area from about 9.8 to 6.9, grazing from about −5.0 to −2.9. Part 4 is about why.
Part 4 · Now all six at once
▶ Run this:
# ── Part 4.1 · Now all six at once ──────────────────────────────────────
full_model <- lm(abundance ~ log_area + log_dist +
log_dist_larger + year_isolated +
grazing + altitude_m,
data = bird_df)
summary(full_model)✏️ Your turn: Which predictors are still significant, and what happened to grazing?
Still significant: _________________________________________
grazing slope: alone ______ two-predictor ______ here ______✏️ Your turn: All three grazing slopes came from the same 56 patches. Is any of them wrong? Say what each one is actually measuring.
____________________________________________________________
____________________________________________________________Part 5 · Why it happened
▶ Run this:
# ── Part 5.1 · Why it happened ──────────────────────────────────────────
bird_df %>%
select(log_area, log_dist, log_dist_larger,
year_isolated, grazing, altitude_m) %>%
cor() %>%
round(2)
vif(full_model)✏️ Your turn: Find the two most strongly correlated predictors, and say in one sentence what that correlation means on the ground.
Pair: __________ and __________ r = __________
On the ground: ______________________________________________✏️ Your turn: All the VIFs are small. Does that mean collinearity is not a problem here?
____________________________________________________________Part 6 · Which model do you keep?
▶ Run this — test whether the four extra predictors earned anything, then compare the two models on AIC:
# ── Part 6.1 · Which model do you keep? ─────────────────────────────────
anova(reduced_model, full_model)
tibble(
model = c("six predictors", "log_area + grazing"),
adj_r2 = c(summary(full_model)$adj.r.squared,
summary(reduced_model)$adj.r.squared),
AIC = c(AIC(full_model), AIC(reduced_model))
)✏️ Your turn: Do the two methods agree? Which model do you keep?
anova p = __________ lower AIC = __________
Model I keep: _______________________________________________Part 7 · Check the assumptions
▶ Run this:
# ── Part 7.1 · Check the assumptions ────────────────────────────────────
par(mfrow = c(2, 2))
plot(reduced_model)
par(mfrow = c(1, 1))
shapiro.test(residuals(reduced_model))✏️ Your turn: One line per assumption.
Linearity: __________________________________________________
Equal variance: _____________________________________________
Normality (Shapiro p = ______): _____________________________
Influential points: _________________________________________Part 8 · Grazing as a factor
Grazing has been a number, 1 to 5, which forces every step to cost the same. Now let it be a category so each level is estimated separately.
▶ Run this — make it a factor and look at the order R chose:
# ── Part 8.1 · Grazing as a factor ──────────────────────────────────────
bird_df <- bird_df %>%
mutate(grazing_class = factor(grazing_class))
levels(bird_df$grazing_class)
summary(lm(abundance ~ log_area + grazing_class, data = bird_df))✏️ Your turn: Which level is R using as the reference, and why did it pick that one?
Reference level: __________ Why: _________________________▶ Run this — put the least-grazed level first instead:
# ── Part 8.2 · Grazing as a factor ──────────────────────────────────────
bird_df <- bird_df %>%
mutate(grazing_class = fct_relevel(grazing_class,
"minimal", "light", "moderate",
"heavy", "intense"))
levels(bird_df$grazing_class)
grazing_model <- lm(abundance ~ log_area + grazing_class, data = bird_df)
summary(grazing_model)✏️ Your turn: Write down the four grazing coefficients and their p-values.
light: __________ p = __________
moderate: __________ p = __________
heavy: __________ p = __________
intense: __________ p = __________✏️ Your turn: The R² and F-statistic are identical in 8.1 and 8.2 — only the coefficients moved. Given that, what did fitting grazing as the number 1–5 (one slope of about −2.9 birds per step, back in Part 3) hide from you?
____________________________________________________________
____________________________________________________________Part 9 · The figure you hand in
▶ Run this:
# ── Part 9.1 · The figure you hand in ───────────────────────────────────
final_plot <- ggplot(bird_df, aes(log_area, abundance,
colour = grazing_class)) +
geom_point(size = 2.5) +
geom_smooth(method = "lm", se = FALSE) +
labs(
title = "Bird abundance rises with patch area and falls under intense grazing",
subtitle = "56 remnant forest patches, Victoria, Australia (Loyn 1987)",
x = "log10 patch area (ha)",
y = "Bird abundance",
colour = "Grazing"
) +
theme_regular(base_size = 12)
final_plot
ggsave(
"figures/bird_regression.pdf",
plot = final_plot,
width = 7,
height = 4.5,
units = "in"
)✏️ Your turn: Open the PDF and check the legend order runs minimal → intense, not alphabetically.
Review and checkpoint
At this point you can:
📤 What to turn in
11_multiple_regression.Rfigures/bird_regression.pdf- This worksheet, with your written answers
Homework
Homework 6 — Multiple Regression: When Predictors Tell the Same Story
Multiple regression takes two class days, so this one is due the following week.
Seventy-three North American grassland sites, six climate and geography predictors, and one response: the proportion of the plant community that is C3 grass. Two of the predictors are so tangled that dropping one reverses the sign of the other — so the honest conclusion is about what the data cannot separate, which is exactly the move you practised in class.
📥 Grassland data 📥 Homework skeleton (.R)
The homework has its own skeleton, separate from the one you used in class. Its header opens with the eight-step process for testing a multiple regression — the same order you worked through today — so use that as your checklist rather than starting from a blank script.
Data source: Quinn & Keough, Experimental Design and Data Analysis for Biologists, Box 6.1 (paruelo.csv); originally Paruelo & Lauenroth (1996), Ecological Applications.