Worked solutions — Linear Models, Links, and GLMs

Exercise M1 — Orangutan line transects (Exponential)

Setup. Detection function $f(x) = \lambda e^{-\lambda x}$ with $\lambda = 0.02$ per meter.

The survival function (probability of not yet being detected by distance $x$) is

At $x = 50$:

lambda <- 0.02
exp(-lambda * 50) # 0.3679

Biological interpretation. A 37% probability of missing an orangutan that is 50 m off the trackline is huge. This is exactly why distance-sampling analyses must correct raw counts by a fitted detection function instead of treating raw counts as abundance estimates. The "effective strip half-width" — the integral of $S(x)$ — gives you a single number to multiply transect length by to get the area effectively searched. For $\lambda = 0.02$, that effective half-width is $1/\lambda = 50$ m, even though the truncation distance might be 100 m or more.

Exercise M2 — Linear transformations of Normals

Let $Z \sim \mathcal{N}(0,1)$ and $Y = \mu + \sigma Z$. Show each step and name the rule used.

Mean.

$E[Y] = E[\mu + \sigma Z]$   (substitute $Y=\mu+\sigma Z$)

$\quad = \mu + \sigma\,E[Z]$   (by linearity of expectation, $E[aX+b]=a\,E[X]+b$, with $a=\sigma,\ b=\mu$)

$\quad = \mu + \sigma\cdot 0 = \mu$   (because $E[Z]=0$ for a standard Normal).

Variance.

$\mathrm{Var}(Y) = \mathrm{Var}(\mu + \sigma Z)$   (substitute $Y=\mu+\sigma Z$)

$\quad = \mathrm{Var}(\sigma Z)$   (adding the constant $\mu$ does not change spread: $\mathrm{Var}(X+b)=\mathrm{Var}(X)$)

$\quad = \sigma^2\,\mathrm{Var}(Z)$   (by the scaling rule $\mathrm{Var}(aX+b)=a^2\,\mathrm{Var}(X)$, with $a=\sigma$)

$\quad = \sigma^2\cdot 1 = \sigma^2$   (because $\mathrm{Var}(Z)=1$ for a standard Normal).

The shape stays Gaussian because the characteristic function (or MGF) of $\mu + \sigma Z$ is still that of a Normal — a linear transformation of a Normal random variable is Normal.

set.seed(1)
z <- rnorm(1e6)
y <- 10 + 3 * z
mean(y); sd(y) # ~ 10 and ~ 3
hist(y, breaks=80) # symmetric bell

Biological interpretation. This is the mathematical reason every regression diagnostic plot puts standardized residuals on a $\mathcal{N}(0,1)$ axis. Any Normal can be rescaled to a standard Normal — so all the percentile thresholds (the $\pm 1.96$ for a 95% interval, the $\pm 2.58$ for 99%) are universal once you standardize. Students who internalize this also understand why standardizing predictors is so often useful: it makes effect-size comparisons across covariates dimensionally meaningful.

Exercise M3 — Beta-Binomial trail cameras

(a) Mean of Beta$(0.03, 0.97)$: $\alpha/(\alpha + \beta) = 0.03/1 = 0.03$. Same expected detection probability as the plain binomial.

(b) Simulation. Note the structure: one $p_i$ per camera, reused across that camera's 10 nights. The clustering is what creates the overdispersion.

set.seed(1); reps <- 5000
K <- 20; J <- 10 # 20 cameras x 10 nights = 200 camera-nights

bin <- rbinom(reps, K*J, 0.03)
betbin <- replicate(reps, {
 p <- rbeta(K, 0.03, 0.97) # ONE p per camera...
 sum(rbinom(K, J, p)) # ...reused across its 10 nights
})
mean(bin == 0); mean(betbin == 0) # 0.003 vs 0.166
mean(bin); mean(betbin) # 6.0 vs 6.0
sd(bin); sd(betbin) # 2.5 vs 5.6

(c) Why so much "all-zero" mass? With $\alpha + \beta = 1$ the Beta is U-shaped: most cameras draw a $p_i$ near 0, a few near 1. A near-blind camera returns zero on all ten of its nights, so whole deployments come back empty far more often (0.166 vs 0.003) while the mean is unchanged at ~6. The reuse of $p_i$ within a camera is essential. If every camera-night drew its own independent $p$, the mixture would integrate out — each night marginally Bernoulli(0.03) and independent, giving exactly Binomial(200, 0.03), SD 2.4, no overdispersion at all. Heterogeneity alone does not overdisperse; heterogeneity shared across repeated observations does.

Biological interpretation. This is the camera-trap version of what Day 1 taught with mixtures of Poissons: heterogeneity among sampling units inflates variance even when the marginal mean is correct. If your wildlife survey is "we deployed 200 cameras for one night each," the right model is almost never a plain binomial — cameras differ in placement, lure, weather, and species behaviour near each unit. Beta-binomial (or random-effect logistic) corrects the SEs that a plain binomial would underestimate.

Exercise A1 — Sockeye eDNA: choosing the distribution

library(MASS)
sock <- read.csv("data/sockeye_adult.csv", stringsAsFactors = FALSE)
table(sock$Sockeyetype)                          # adult level is "Sockeye_Adults"
d <- subset(sock, Sockeyetype == "Sockeye_Adults")
d$Year <- factor(sub(".*/", "", d$Date))         # "15" or "16"
d <- d[!is.na(d$Count) & !is.na(d$Qcorr_qPCR), ] # complete cases (n = 55)
min_pos <- min(d$Qcorr_qPCR[d$Qcorr_qPCR > 0])    # 1e-05
d$logED <- log(d$Qcorr_qPCR + min_pos)

# 1. Log-normal regression on log(count+1)
f_lm <- lm(log(Count + 1) ~ logED + Year, data = d)
summary(f_lm)                                    # logED slope 0.570 (SE 0.047)

# 2. Plain Poisson
f_p <- glm(Count ~ logED + Year, family = poisson, data = d)
f_p$deviance / f_p$df.residual                   # 2296.6 / 52 = 44.2  (overdispersed!)

# 3. Quasipoisson
f_qp <- glm(Count ~ logED + Year, family = quasipoisson, data = d)
summary(f_qp)$dispersion                          # 51.0

# 4. Negative binomial
f_nb <- glm.nb(Count ~ logED + Year, data = d)
f_nb$deviance / f_nb$df.residual                  # 48.8 / 52 = 0.94  (theta = 0.73)
2^coef(f_nb)["logED"]                             # 1.86

Diagnostic check. The plain Poisson is grossly overdispersed: residual deviance / df $\approx 44$ (any value $\gg 1$ signals overdispersion), so its SE on logED (0.013) is far too small. The quasipoisson keeps the same point estimate but rescales the SEs by $\sqrt{\hat\phi} = \sqrt{51}\approx 7.1$, giving 0.095. The negative binomial provides a proper likelihood (so you can compare via AIC) and a dispersion parameter $\theta = 0.73$; its residual-deviance ratio (0.94) is near 1, a well-calibrated fit. Trust fits 3 and 4, not the naive Poisson.

Coefficient back-transform. The eDNA slope from the NB fit is $\hat\beta = 0.897$ on the $\log(\text{eDNA})$ predictor. A "doubling of the flow-corrected eDNA signal" corresponds to a one-unit change on the $\log_2$ scale, i.e. $\ln 2 = 0.693$ on the natural-log scale. That doubling adds $0.897 \times 0.693 = 0.622$ to the linear predictor; because the link is the log, the additive shift becomes a multiplicative effect on the expected count:

A doubling of the flow-corrected eDNA signal is associated with roughly an 86% increase (×1.86) in expected daily adult sockeye count — the same value the graded Problem 5 reports. Additive on the link scale, multiplicative on the response scale: that is the whole point of the log link, and it is the sentence students should be able to produce unprompted.

Biological interpretation. eDNA is shed proportionally to biomass (per unit volume of water). The fitted slope $\hat\beta = 0.897$ is a little below 1 on the $\log(\text{eDNA})$ scale, indicating a slightly sublinear eDNA-to-count relationship — possibly because fish in dense schools shed less per fish (boundary-layer effects), or because qPCR saturates at high concentrations. A slope above 1 would suggest the opposite — eDNA "amplifying" the count signal, perhaps because larger fish (which dominate the run-timing tails) shed disproportionately more eDNA per individual. The ×1.86 number is what goes in your abstract; the SE matters only insofar as it tells you whether the slope is distinguishable from 0 or from 1.

Exercise A2 — Titanic: long vs aggregated form

long <- read.csv("data/titanic_long.csv")
prop <- read.csv("data/titanic_prop.csv")

f1 <- glm(survived ~ class, family = binomial, data = long)
f2 <- glm(cbind(Yes, No) ~ Class, family = binomial, data = prop)

# Fitted probability per class (crew, first, second, third)
lp <- coef(f1)[1] + c(0, coef(f1)[-1])
plogis(lp)                          # fitted p per class
tapply(long$survived, long$class, mean)   # raw class proportions = same
# titanic_prop is one row per Class x Sex x Age, so aggregate to Class first:
tapply(prop$Yes, prop$Class, sum) / tapply(prop$Yes + prop$No, prop$Class, sum)

Output (each class's absolute log-odds, odds, fitted $\hat p$, and the raw proportion):

Why the two fits agree (part c). Group the passengers by class. For one group of $n$ passengers with $y=\sum_i y_i$ survivors, the one-row-per-person (Bernoulli) likelihood and the aggregated (Binomial) likelihood are

They differ only by the binomial coefficient $\binom{n}{y}$ — the number of orderings of $y$ survivors among $n$ people (e.g. the 285 second-class passengers had $y=118$ survivors). On the log scale that is an additive constant $\log\binom{n}{y}$ that contains no $p$. Since the MLE, its SEs, and every likelihood-ratio test depend only on how the log-likelihood varies with $p$, this constant cancels: the coefficients, fitted probabilities, and SEs are identical. (The coefficient tables can also look different simply because long and prop sort their factor levels differently and pick different reference classes; align the references and the estimates match exactly.)

What the constant shifts — and what it doesn't. This is the step most people garble, so be precise with students. The constant shifts the log-likelihood and the AIC. For these data $\sum_{\text{groups}}\log\binom{n}{y} = 1025.44$ over the class×sex×age cells of titanic_prop.csv, and the Binomial log-likelihood sits exactly that far above the Bernoulli one ($-268.84$ vs $-1294.28$), so the AICs differ by $2 \times 1025.44 = 2050.87$.

It does not shift the residual deviance. Deviance is $2(\ell_{\text{saturated}} - \ell_{\text{model}})$, and $\log\binom{n}{y}$ appears in both terms, so it cancels out entirely. The two deviances do differ — Bernoulli 2588.56 vs Binomial 491.06 — but that gap is 2097.49, which is neither 1025.44 nor 2050.87. Its real source is the change of saturated reference: the Bernoulli form's saturated model has one parameter per passenger and fits every row perfectly ($\ell_{\text{sat}} = 0$), while the Binomial form's has one per cell. The gap is exactly the Bernoulli deviance of the cell-saturated model — the individual-level variation that aggregation stops trying to explain.

n <- prop$Yes + prop$No
sum(lchoose(n, prop$Yes))                 # 1025.437  the constant

logLik(f1); logLik(f2)                    # -1294.278  vs  -268.841   (shifted by 1025.437)
AIC(f1);    AIC(f2)                       #  2596.555  vs   545.682   (shifted by 2050.873)
deviance(f1); deviance(f2)                #  2588.555  vs   491.061   (gap 2097.495 -- neither!)

# the deviance gap is the change of saturated reference:
deviance(glm(survived ~ class*sex*age, family = binomial, data = long))   # 2097.495

The moral is unchanged and worth stating flatly: never compare the deviance or AIC of a Bernoulli fit against an aggregated Binomial fit. The AIC gap is bookkeeping; the deviance gap is a change of yardstick. Neither carries information about model quality. Within a single data form, both remain valid for comparing nested models.

Biological interpretation. The class coefficients trace a steep social gradient: 1st-class odds of survival are about 5× those of 3rd class ($1.66/0.34 \approx 4.9$), or equivalently a log-odds gap of about 1.60 ($+0.51 - (-1.09)$). Translating to probabilities, 1st-class passengers survived at 62.5% vs 25.2% in 3rd — a 37 percentage-point gap. The model recovers the raw proportions exactly because with a single categorical predictor and no other constraints, the MLE just is the empirical proportion in each group. This is a useful sanity check: any time your saturated logistic model disagrees with the empirical fractions, something is wrong with your code.

Exercise A3 — Titanic with interactions

library(effects); library(emmeans)
long <- read.csv("data/titanic_long.csv", stringsAsFactors = TRUE)  # age is a factor: adult/child
f3 <- glm(survived ~ class * sex + age, family = binomial, data = long)
summary(f3)
plot(allEffects(f3))
emmeans(f3, pairwise ~ class | sex, type = "response")

(a) Largest pairwise odds ratio. The female 1st-class vs female 3rd-class contrast is by far the largest: OR $\approx 48.7$ (log-odds difference $\approx 3.89$, $p < 0.0001$). Averaged over age, first-class women survived at $p \approx 0.98$ and third-class women at $p \approx 0.55$. For contrast, the same first-vs-third contrast among men is only OR $\approx 2.8$ — the "women and children first" protocol compressed the class gradient for men.

(b) Why not TukeyHSD on an aov? TukeyHSD treats the model as a Gaussian linear model with constant variance and an identity link. On a binary outcome, the variance is $p(1-p)$ — heteroskedastic by construction — and the response is on the wrong scale (0/1) for Normal-theory inference. emmeans applies the correct binomial likelihood and lets you choose whether to report contrasts on the link scale (log-odds) or the response scale (probability), with the correct delta-method standard errors.

(c) Female 1st-vs-3rd, adult passenger, probability difference. In this dataset age is a two-level factor (adult/child), so we predict for an adult female:

predict(f3, newdata = data.frame(class = c("first","third"), sex = "female", age = "adult"),
 type = "response")
#  first  third
#  0.972  0.419

An adult first-class woman was modeled at $p = 0.972$ versus $p = 0.419$ for an adult third-class woman — about 55 percentage points higher. On the link scale that is a log-odds gap of $3.56 - (-0.33) = 3.89$, i.e. $e^{3.89}\approx 49$, reconciling with the odds ratio in part (a): the OR and the probability difference are the same effect read on two scales.

Biological interpretation. The interaction encodes the famous "women and children first" effect being selectively applied. The protective effect of being female (the main sex coefficient) is amplified for higher-class passengers and partially erased for 3rd-class — because the structural barriers (gates, distance to boats, language) overwhelmed the sex-prioritized evacuation in 3rd class. Reading the model on the odds-ratio scale tells you the social pattern; reading it on the probability scale tells you the human cost.

Notes on the log-link examples (L1–L8)

Notes on the logit-link examples (B1–B8)

Try-it exercises (log + logit)

Log-link "Try" answers

• Try 1. $e^{0.12 \cdot 5} = e^{0.6} = 1.82$. Five degrees of warming multiplies expected counts by 1.82.
• Try 2. $\beta = \ln(1.5) = 0.405$ on the $\log_2$ predictor.
• Try 3. Point estimate $e^{0.40} = 1.49$; CI $(e^{0.18}, e^{0.62}) = (1.20, 1.86)$. Asymmetric on the response scale.

Logit-link "Try" answers

• Try L1. $e^{0.69} = 2.00$. The odds double.
• Try L2. $p = $ plogis(c(-1, 2, 5)) = (0.269, 0.881, 0.993).
• Try L3. $\text{logit}(0.75) = \ln 3 = 1.099$; $x = (1.099 + 0.5)/0.03 = 53.3$.
• Try L4. Marginal effect at $p = 0.2$: $0.5 \cdot 0.2 \cdot 0.8 = 0.08$, or 8 percentage points per unit $x$.