Activity: GLM & Logistic Regression
Generalized linear models
In-class Activity 15: Generalized Linear Models
Today’s Objectives
- Fit the week-9 ANOVA on count data and watch its assumptions fail
- Say what a GLM’s three components are — distribution, predictors, link function
- Check that a Gaussian GLM reproduces
lm()exactly - Fit a Poisson GLM, and find out why the obvious choice for counts still fails
- Measure overdispersion, and say what it does to your p-values
- Fit a negative binomial GLM and confirm it fixed the problem
- Run the same escalation on a continuous predictor, on your own
- Fit a logistic regression on a binary response and interpret an odds ratio
📥 Script skeleton (.R) 📥 Polis lizard data
- Save the skeleton as
scripts/15_glm_logistic_regression.R, andpolis.csvin yourdata/folder. - The Galapagos data comes from the
farawaypackage — no file needed.
- 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 two data reads — and nothing else. Every line of statistics you type yourself.
- Section 6 has no code here either. By then you will have run the same three steps three times, so you do them once more on your own.
- 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, 9, and 10 — is graded. Those answers go into the script as
#comments. - Hand in
scripts/15_glm_logistic_regression.R. Run it top to bottom (Ctrl/Cmd + Shift + Enter) first — it must run with no errors.
Today is a chain of three models where each one fails for a different reason, and the failure is the content. Before each model, say out loud whether you expect it to work — then look at the diagnostic, not the p-value.
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 10, which is the graded deliverable.
Where this fits
For the last three weeks every model has been a linear model — one-way ANOVA, factorial, nested. All of them assume the residuals are normal with constant variance. That is a fine assumption for body mass or algae cover.
It is a bad assumption for a count. A count is a whole number, it cannot go below zero, and its variance usually grows with its mean. Today is what you do instead.
A generalized linear model keeps the linear part and swaps out two things:
| Component | What it is | Example |
|---|---|---|
| Random | the response’s distribution | Poisson for counts |
| Systematic | the predictors | size_cat, area |
| Link | connects the mean to the predictors | log for counts, logit for binary |
| Response looks like | Family |
|---|---|
| continuous, symmetric | gaussian — this is just lm() |
| counts | poisson |
| counts, variance > mean | negative binomial (glm.nb) |
| binary — present/absent | binomial(link = "logit") |
| positive, continuous, skewed | Gamma |
🧩 Block 1 — The assumptions fail (after lecture Part 1)
Sections 1–2: meet the count data, then fit last week’s model to it on purpose.
Section 1 · The data
Written for you in the skeleton — run it and read it.
# ---- Libraries -------------------------------------------
library(janitor)
library(car) # Anova()
library(MASS) # glm.nb()
library(performance) # check_overdispersion()
library(DHARMa) # simulateResiduals()
library(emmeans)
library(broom)
library(parameters) # model_parameters()
library(pscl) # pR2()
library(ResourceSelection) # hoslem.test()
library(faraway) # the gala data
library(tidyverse)
options(scipen = 999)The gala data: 30 Galapagos islands, used to test MacArthur–Wilson island biogeography.
| Column | What it is |
|---|---|
spp |
number of plant species — a count |
endemics |
number of endemic species — a count |
area |
island area, km² |
elevation |
maximum elevation, m |
nearest |
distance to the nearest island, km |
# ---- 1 The data ------------------------------------------
g_df <- gala %>%
clean_names() %>%
rename(spp = species) %>%
mutate(
size_cat = case_when(
area < 1 ~ "small",
area >= 1 & area < 100 ~ "medium",
area >= 100 ~ "large"
),
size_cat = factor(size_cat, levels = c("small", "medium", "large"))
) %>%
filter(area < 3000)
head(g_df)
count(g_df, size_cat)filter(area < 3000) drops Isabela, which is about 20× larger than anything else and flattens every plot. size_cat bins area into three groups so we start from a categorical predictor — the same shape of question as week 9’s ANOVA.
✏️ Q1: spp is a count. Name two things that are true of counts but not true of a normal distribution.
✏️ Q2: How many islands are in each size category? Is this balanced?
Section 2 · Why last week’s ANOVA will not do
Fit it the week-9 way, on purpose, and let the assumptions fail.
▶ Type this yourself:
# ---- 2 Why last week's ANOVA will not do -----------------
ggplot(g_df, aes(x = spp)) +
geom_histogram(binwidth = 25)
ggplot(g_df, aes(x = size_cat, y = spp)) +
geom_boxplot()🔮 Predict first: for a normal distribution the variance has nothing to do with the mean. Before you run the next chunk — guess the mean of
spp, then guess its variance.
▶ Type this yourself:
g_df %>%
summarise(
mean_spp = mean(spp),
var_spp = var(spp),
ratio = var_spp / mean_spp
)▶ Type this yourself:
lm_model <- lm(spp ~ size_cat, data = g_df)
summary(lm_model)
Anova(lm_model, type = 3)
shapiro.test(residuals(lm_model))
leveneTest(spp ~ size_cat, data = g_df)✏️ Q3: Write down the mean, the variance, and both assumption p-values. Did either assumption pass?
✏️ Q4: Anova() gave size_cat a p-value of about 0.0000035. Did this model earn the right to that number?
✏️ Q5: The variance of spp is more than 100× its mean. A normal distribution lets the variance be any number it likes, so that ratio on its own is not a violation. Look at the boxplot instead — which assumption is the ratio actually breaking?
✏️ Q6: Look at the histogram. If you fitted a normal curve to the “small” islands and asked for a 95% interval, what impossible value would it predict?
💡 Key idea: this is the moment the ANOVA toolkit runs out. Nothing is wrong with the question — island size probably does drive species richness — but the model being used to answer it is making promises the data does not keep. The rest of today is three attempts to fix that.
🧩 Block 2 — Three families, one question (after lecture Parts 2–4)
Sections 3–5: Gaussian (the sanity check), Poisson (fails), negative binomial (works).
Section 3 · A Gaussian GLM is the same linear model
Before trusting glm() with an unfamiliar distribution, check that it reproduces the one you already know.
🔮 Predict first: a Gaussian GLM with an identity link is a linear model. Will the coefficients be identical to
lm_model, or just similar?
▶ Type this yourself:
# ---- 3 A Gaussian GLM is the same linear model -----------
gauss_model <- glm(spp ~ size_cat, data = g_df,
family = gaussian(link = "identity"))
summary(gauss_model)
Anova(gauss_model, type = "III", test = "F")✏️ Q7: Compare the coefficients and the size_cat p-value to section 2. Identical, or just close?
✏️ Q8: So if a Gaussian GLM is just lm() again, what is glm() actually adding? (What did you have to name in the glm() call that lm() never asked for?)
💡 Key idea:
lm()is the special case ofglm()where the family is Gaussian and the link is the identity. Everything that follows changes exactly those two arguments. You are not learning a new model — you are learning two knobs on the one you already have.
Section 4 · Poisson GLM — the obvious choice, and it fails
Poisson is the distribution for counts. It has only one parameter, which means its mean and its variance are forced to be the same number. Remember the ratio you computed in Q3.
🔮 Predict first: the Poisson distribution insists variance = mean. Your data’s ratio was over 100. Will the overdispersion check come back clean?
▶ Type this yourself:
# ---- 4 Poisson GLM - the obvious choice, and it fails ----
poiss_model <- glm(spp ~ size_cat, data = g_df,
family = poisson(link = "log"))
summary(poiss_model)
Anova(poiss_model, type = "III", test = "LR")A Poisson GLM cannot use an F-test, so we use a likelihood ratio test instead — it compares the fit of the full model against a model without the predictor.
▶ Type this yourself:
check_overdispersion(poiss_model)
sim_poiss <- simulateResiduals(poiss_model, n = 1000)
plot(sim_poiss)
testDispersion(sim_poiss)⚠️ Watch out!
shapiro.test(residuals(poiss_model))will run, and it will give you a number. The number is meaningless — normality of raw residuals is not a Poisson assumption. R will not stop you from running a test that does not apply. This is why you check the diagnostic that matches the model, not the one you remember.
✏️ Q9: Why can a Poisson GLM not use an F-test? (What does an F-test need in its denominator that this model never estimates?)
✏️ Q10: The dispersion ratio should be about 1. What did you get? Given Q3, is that a surprise?
✏️ Q11: Overdispersion makes standard errors too small, which makes test statistics too big. Write down the size_cat LR chi-square here — you will compare it to section 5’s.
✏️ Q12: Look at the DHARMa Q-Q plot. Is the problem subtle?
💡 Key idea: Poisson was the right instinct and still the wrong model. “It’s count data, so use Poisson” is not enough — Poisson makes a second, much stronger claim, that the variance equals the mean. Ecological counts almost never oblige, because real populations cluster.
Section 5 · Negative binomial — the one that works
The negative binomial adds a second parameter, theta, which lets the variance exceed the mean by however much it needs to.
▶ Type this yourself:
# ---- 5 Negative binomial - the one that works ------------
nb_model <- glm.nb(spp ~ size_cat, data = g_df)
summary(nb_model)
Anova(nb_model, type = "III", test = "LR")▶ Type this yourself:
check_overdispersion(nb_model)
sim_nb <- simulateResiduals(nb_model, n = 1000)
plot(sim_nb)
AIC(poiss_model, nb_model)▶ Type this yourself:
nb_emm <- emmeans(nb_model, ~ size_cat, type = "response")
nb_emm
pairs(nb_emm)
multcomp::cld(nb_emm, Letters = letters)✏️ Q13: The NB dispersion ratio and its p-value — is the overdispersion gone?
✏️ Q14: Compare the size_cat LR chi-square to the Poisson one from Q11. By what factor was the Poisson version inflated?
✏️ Q15: The AIC gap between the two models is enormous. Same data, same predictor, same link — so where did all that fit come from?
✏️ Q16: type = "response" back-transforms the means from the log scale. Look at the confidence limits: why are they no longer symmetric around each mean?
✏️ Q17: Read the letter display. Are all three size categories distinguishable?
💡 Key idea: you fitted three models to one question and got three different p-values, all of them tiny. The number was never the problem. What changed across the three is whether you are entitled to it — and only the third model’s diagnostics say yes.
🧩 Block 3 — On your own, then a binary response (after lecture Part 5)
Section 6 is yours with no code at all. Section 7 is a new kind of response.
Section 6 · Do it again with a continuous predictor (you write this one)
No code on this page. You have now run the same three steps three times, so run them once more without the script in front of you.
Binning area into three categories threw information away. Use area directly.
Question: how does plant species richness change with island area?
Work through:
- Plot
sppagainstarea. lm(spp ~ area)— and check the two assumptions.- The Poisson version, with an overdispersion check.
- The negative binomial version, with an overdispersion check, DHARMa, and the LR test.
- One plot of
sppagainstareawith the NB fit drawn on it.geom_smooth(method = "glm.nb")will do it. - Decide which of the three you would report.
💡 Hint: the only things that change from sections 2–5 are the formula and the data column. If you get stuck, go back and read what you typed, not this page.
✏️ Q18: Which model did you end up reporting, and what was the single piece of output that decided it?
✏️ Q19: Did binning area into three groups in section 1 cost you anything? Compare what size_cat told you to what area tells you.
✏️ Q20: Look at your NB fit on the plot. It is a curve, not a line. Where does the curve come from — the distribution, or the link function?
Section 7 · Logistic regression — a binary response
A count was one kind of non-normal response. A yes/no is another. Logistic regression models the probability of a 1, using the logit link to stretch probability (stuck between 0 and 1) onto a scale that runs from −∞ to +∞:
\[\log\left(\frac{\pi}{1-\pi}\right) = \beta_0 + \beta_1 x\]
Polis et al. (1998) recorded whether the lizard Uta was present on 19 islands in the Gulf of California, against each island’s perimeter-to-area ratio — a measure of how exposed the island’s shoreline is. These are the real data from the paper, as reproduced in Quinn & Keough Box 13-1.
Written for you in the skeleton — the "." missing-value codes need handling, and Peninsula has no paratio:
# ---- 7 Logistic regression - a binary response -----------
polis_df <- read_csv("data/polis.csv", na = c(".", "")) %>%
clean_names() %>%
drop_na(paratio, pa)
polis_df
count(polis_df, uta)▶ Type this yourself:
ggplot(polis_df, aes(x = paratio, y = pa)) +
geom_point()🔮 Predict first: exposed shorelines are hard places to be a lizard. Will the slope be positive or negative?
▶ Type this yourself:
lizard_model <- glm(pa ~ paratio, data = polis_df,
family = binomial(link = "logit"))
summary(lizard_model)▶ Type this yourself — the likelihood ratio test, against a model with no predictor:
null_model <- glm(pa ~ 1, data = polis_df, family = binomial(link = "logit"))
anova(null_model, lizard_model, test = "Chisq")▶ Type this yourself — the odds ratio, three ways:
exp(coef(lizard_model))
exp(confint(lizard_model, "paratio"))
model_parameters(lizard_model, exponentiate = TRUE)
tidy(lizard_model, exponentiate = TRUE, conf.int = TRUE)▶ Type this yourself — goodness of fit:
r2(lizard_model)
pR2(lizard_model)
hoslem.test(lizard_model$y, fitted(lizard_model), g = 10)▶ Type this yourself — the fitted probability curve:
pred_df <- data.frame(
paratio = seq(min(polis_df$paratio), max(polis_df$paratio), length.out = 100)
)
pred_df$prob <- predict(lizard_model, newdata = pred_df, type = "response")
ggplot() +
geom_point(data = polis_df, aes(x = paratio, y = pa)) +
geom_line(data = pred_df, aes(x = paratio, y = prob))✏️ Q21: Report the slope and the odds ratio with its CI. Does the CI include 1, and why does that matter more than whether it includes 0?
✏️ Q22: Say what the odds ratio means about Uta and shoreline exposure — in lizards, not in log-odds. “For every one-unit increase in perimeter-to-area ratio, the odds of finding Uta …”
✏️ Q23: For the Hosmer–Lemeshow test you want a non-significant p-value. Why is that backwards from every other test you have run this semester?
✏️ Q24: r2() reports a pseudo-R². Why “pseudo” — what can’t you compute here that you could in a regression?
✏️ Q25: Look at your curve. At roughly what perimeter-to-area ratio is Uta a coin flip?
💡 Key idea: the link function is doing the work in both halves of today. A
loglink keeps predicted counts above zero; alogitlink keeps predicted probabilities between 0 and 1. In both cases the linear part of the model is unchanged — you are still fitting a straight line, just on a scale where a straight line cannot predict something impossible.
Write-up — graded
These go into the bottom of your script, as # comments after ANSWER:.
Section 8 · Report the count model
✏️ Q26: Fill in all three models of spp ~ size_cat:
model test statistic p dispersion
lm / Gaussian F -
Poisson LR chisq
negative binomial LR chisq✏️ Q27: Two or three sentences reporting the negative binomial result, with your own statistic, df, and p. Say which model you fitted and why that one.
✏️ Q28: Someone reports the Poisson p-value for size_cat. What have they done wrong, and in which direction is their answer wrong?
✏️ Q29: Section 2’s ANOVA also gave a significant size_cat effect. So did the whole exercise change the conclusion, or only the honesty of it? Explain — and say why that still matters.
Section 9 · Report the logistic model
✏️ Q30: Slope, odds ratio and its 95% CI, the LR test result, and a pseudo-R².
✏️ Q31: One sentence an ecologist would actually write about perimeter-to-area ratio and the chance of finding Uta.
Section 10 · The final figure
One publication-quality figure. Pick either the NB fit of spp on area from section 6, or the logistic curve from section 7. Show the raw data under the fitted line. Proper axis labels, a real theme, readable sizing — this one is graded on presentation. Export it with ggsave().
# ---- 10 The final figure ---------------------------------
# ... your plot here, then:
# ggsave("figures/15_glm_fit.pdf", width = 7, height = 5)✏️ Q32: Why would a straight line from lm() have been the wrong shape for the figure you chose? Name the impossible value it would have predicted.
Review and checkpoint
At this point you should be able to:
Run your whole script top to bottom with Ctrl/Cmd + Shift + Enter. It must run clean.
💡 Key idea: the chain you ran today — fit it, check the diagnostic that matches the model, change one assumption, refit — is the whole method. Next week the thing you change is not the distribution but the error structure, with random effects.