Continuous distributions, linear models, and generalized linear models

Where the multiplier comes from. Don't memorize "exponentiate the coefficient" — derive it once, because the same move reads every log-link slope. Start from what the model literally says:

and exponentiate both sides to get onto the count scale:

Now ask the actual question: what happens to the expected count when you move one km further upstream, $\text{km} \to \text{km}+1$? Take the ratio of the new mean to the old so the intercept cancels:

The baseline location drops out — only $e^{\beta_1}$ survives. That is the general rule worth keeping: each one-unit step in a raw predictor multiplies the mean by $e^{\beta_1}$, and a step of $d$ units multiplies it by $e^{\beta_1 d} = (e^{\beta_1})^d$. Additive on the log scale ($\beta_1$ per km), multiplicative on the count scale ($e^{\beta_1}$ per km).

Only now do the numbers enter. With $\hat\beta_1 = -0.05$, $e^{-0.05} = 0.951$, so each additional km upstream drops expected redd counts by ~4.9%. Over 14 km the factor compounds to $(e^{-0.05})^{14} = e^{-0.70} = 0.50$ — half as many redds as at the stream mouth.

The half-distance follows from the same rule: set $e^{\beta_1 d} = \tfrac12$ and solve, $d = \ln(2)/|\beta_1| = \ln(2)/0.05 = 13.9$ km. This "half-life" language is intuitive for biologists who already think in exponential decay.

Same log link as Example 1, so the same rule applies: one Beaufort unit multiplies the expected detection rate by $e^{\hat\beta} = e^{-0.41} = 0.664$ — a ~34% drop per unit. The new twist is a jump of more than one unit, and here the $(e^\beta)^d$ half of the rule earns its keep.

Going from sea state 1 to sea state 4 is a step of $d = 3$ units. Compounding the per-unit factor three times:

The two routes are the same arithmetic: multiply the per-unit factor three times, $0.664^3 = 0.293$, or add $3\beta = -1.23$ on the log scale and exponentiate once, $e^{-1.23} = 0.293$. Either way it's a 71% drop in detections across those three sea-state units — a strong "detection-availability" effect that dominates survey logistics: a slightly rougher sea, and most of your animals stop being detectable.

With a factor predictor R codes the level as a 0/1 dummy, so there is no "per unit" to compound — the coefficient jumps you straight from one group to the other. Write the two group means the model implies:

because the dummy is $0$ for control and $1$ for fertilizer. Exponentiate and take the ratio, and the intercept (the control level) cancels:

So the general rule for a categorical log-link coefficient: $e^{\beta}$ is exactly the ratio of that level's mean to the reference level's mean — the fold change itself, no compounding involved. Plug in $\hat\beta_1 = 0.85$: $e^{0.85} = 2.34$. The fertilizer treatment more than doubles the expected caterpillar count relative to control. This is the "fold change" language ecologists already speak — they just have to remember that the raw coefficient is its logarithm.

The interaction changes the slope in temp, so the trick is to collect the temp terms before exponentiating. Write the linear predictor with the trawl dummy set to its two values. For gillnet (reference, dummy $= 0$) only the main effect survives; for trawl (dummy $= 1$) the interaction term switches on:

Group the coefficients that multiply temp: for gillnet that is just $\beta_{\text{temp}}$; for trawl it is $\beta_{\text{temp}} + \beta_{\text{int}}$. That grouped quantity is the log-scale slope, and the per-degree multiplier is its exponential — the general rule being the temperature effect for any gear is $e^{(\text{sum of every coefficient attached to temp for that gear})}$.

Now the numbers:

(a) Gillnet slope $= \beta_{\text{temp}} = -0.05$, so $e^{-0.05} = 0.951$ → a 4.9% drop in catch per degree of warming.

(b) Trawl slope $= \beta_{\text{temp}} + \beta_{\text{int}} = -0.05 + 0.08 = +0.03$, so $e^{0.03} = 1.030$ → a 3% increase per degree.

The signs are opposite: warming thins the gillnet catch but lifts the trawl catch. That reversal lives entirely in the interaction term — it is invisible in either main effect alone.

Why $e^\beta$ is an odds ratio. The logit link says the model is linear not in $p$ but in the log-odds. Start from what it literally states:

Exponentiate both sides — the left side is now the odds itself:

Ask the same ratio question we used for the log link, but now it lands on the odds. Bump canopy by one point and divide new odds by old:

The baseline odds cancel — so the general rule is each one-unit step in $x$ multiplies the odds by $e^{\beta_1}$, and a $d$-unit step multiplies them by $(e^{\beta_1})^d = e^{\beta_1 d}$. Identical structure to the log link, except the thing being multiplied is odds, not the mean. Plug in $\hat\beta_1 = 0.06$: $e^{0.06} = 1.062$, so each canopy point multiplies occupancy odds by 1.062; a 10-point rise multiplies them by $e^{0.6} = 1.82$.

Connecting the scales: why the probability jump is not constant. The odds ratio $e^{\beta_1}$ is the same everywhere on the curve, but the probability change it buys depends on where you start, because you have to convert odds back to $p$ via $p = \text{odds}/(1+\text{odds})$. Take a site at $p = 0.3$, i.e. odds $= 0.3/0.7 = 0.43$. A 10-point canopy increase moves the odds to $0.43 \cdot 1.82 = 0.78$, which is $p = 0.78/1.78 = 0.44$ — a 14 percentage-point jump. The same $\times 1.82$ applied to a site already at $p = 0.95$ would barely move it. A constant odds ratio, a baseline-dependent probability change: that is the logit link's signature, and the reason you must state the baseline when you report an effect.

As in Example 3 the predictor is a factor, so the coefficient carries you between two groups rather than per unit. The three "scales" are just three stops as you undo the logit link — read them left to right, each one algebra away from the last. Write the two group means the model implies:

with females the reference (dummy $= 0$) and males the dummy $= 1$.

• Log-odds (subtract): the two equations differ by exactly $\beta$, so males have log-odds $0.4$ lower than females. This is the raw coefficient.
• Odds ratio (exponentiate): subtracting then exponentiating gives $\dfrac{\text{odds}_m}{\text{odds}_f} = e^{\beta} = e^{-0.4} = 0.67$ — exactly the categorical-coefficient rule from Example 3, now on the odds scale. Males' odds of fledging are 67% of females'.
• Probability (convert): the odds ratio alone doesn't give a probability — you need a baseline. If females fledge with $p_f = 0.6$, their odds are $0.6/0.4 = 1.5$, so males' odds are $0.67 \cdot 1.5 = 1.0$, i.e. $p_m = 1.0/(1+1.0) = 0.5$. A 10-percentage-point gap.

Same effect, three views: a fixed log-odds difference and a fixed odds ratio, but a probability gap that only exists once you name where the females sit.

Make "steepness" precise by differentiating the curve. The linear predictor is $\eta = \beta_0 + \beta_1 x$, and $p(\eta) = e^\eta/(1+e^\eta)$. The probability change per unit $x$ is, by the chain rule,

(The middle factor is the tidy fact that the logistic curve's slope is $\dfrac{dp}{d\eta} = p(1-p)$.) So a fixed coefficient $\beta_1$ produces a probability change that is scaled by $p(1-p)$ — and that factor is not constant: it peaks at $p = 0.5$, where $p(1-p) = 0.25$, and collapses toward $0$ as $p \to 0$ or $p \to 1$. That is exactly why one unit of $x$ moves $p$ a lot in the middle and barely at the ends.

The peak value $0.25$ is also where the divide-by-4 rule comes from: at the steepest point the probability change per unit $x$ is $\beta_1 \cdot 0.25 = \beta_1/4$, an upper bound on how fast the curve can move. Biologically this is the right shape — a 10-degree warming shifts survival sharply when survival is near 50%, but hardly touches it when survival is already 99% and the curve has gone flat.

Derive the survival function — two ways. "Not yet detected at distance $x$" means the detection distance $X$ exceeds $x$, so we want the survival function $S(x) = P(X \ge x)$. There are two routes to it, and it is worth seeing both.

Way 1 — from first principles (integrate the density). The Exponential$(\lambda)$ density is $f(t) = \lambda e^{-\lambda t}$ for $t \ge 0$. Accumulate all the probability beyond $x$ by integrating from $x$ outward:

Way 2 — from the CDF you already know (no integral). Once you recognize that a "distance to first detection under a constant per-meter detection hazard" is an Exponential random variable, you can skip the calculus. The Exponential CDF is a standard result, $F(x) = P(X \le x) = 1 - e^{-\lambda x}$, and the survival function is simply its complement:

Same answer, and Way 2 takes one line — the payoff for naming the distribution instead of re-integrating. (Way 1 is really just the on-the-spot derivation of that CDF; the two are the same calculation, done from different starting points.)

The general rule to carry away: for any exponential process the survival function is a pure exponential decay $e^{-\lambda x}$, so the fraction still undetected drops by the constant factor $e^{-\lambda}$ per meter — and, exactly as in the log-link redd example, there is a "detection half-distance" $\ln(2)/\lambda$ at which half the animals remain undetected. Here $\ln(2)/0.02 = 34.7$ m.

Now plug in $\lambda = 0.02$ per meter and $x = 50$ m, so $\lambda x = 1$:

A 50-m-distant orangutan has about a 37% chance of being missed — just past the half-distance, so a bit more than half are still detectable.

Mean — one step at a time.

$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$, here with $a=\sigma,\ b=\mu$)

$\quad = \mu + \sigma\cdot 0$   (because $E[Z]=0$ for $Z\sim\mathcal N(0,1)$)

$\quad = \mu.$

Variance — one step at a time.

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

$\quad = \mathrm{Var}(\sigma Z)$   (adding a constant shifts location but not 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)$, here $a=\sigma$)

$\quad = \sigma^2\cdot 1$   (because $\mathrm{Var}(Z)=1$ for $Z\sim\mathcal N(0,1)$)

$\quad = \sigma^2.$

Shape. The distribution stays Normal because a linear transformation of a Gaussian is Gaussian — the moment-generating function of $\mu+\sigma Z$, namely $e^{\mu t + \sigma^2 t^2/2}$, is exactly the MGF of a $\mathcal N(\mu,\sigma^2)$. So $Y\sim\mathcal N(\mu,\sigma^2)$.

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 — same mean
sd(bin); sd(betbin) # 2.5 vs 5.6 — about 5x the variance

(c) With $\phi$ this small the Beta piles $p_i$ up near 0 or 1, so most cameras are near-blind while a few are hot. A camera that draws a tiny $p_i$ contributes zero on all ten of its nights, so all-zero deployments become common (0.166 vs 0.003) — even though the mean is unchanged at ~6. That is overdispersion: same mean, fatter tails.

f_lm <- lm(log(Count + 1) ~ logED + Year, data = d)
f_p  <- glm(Count ~ logED + Year, family = poisson,      data = d)
f_qp <- glm(Count ~ logED + Year, family = quasipoisson, data = d)
f_nb <- MASS::glm.nb(Count ~ logED + Year, data = d)

# overdispersion diagnostic for the Poisson
f_p$deviance / f_p$df.residual              # 2296.6 / 52 = 44.2
summary(f_qp)$dispersion                     # 51.0
f_nb$deviance / f_nb$df.residual             # 48.8 / 52 = 0.94  (theta = 0.73)
2^coef(f_nb)["logED"]                         # 1.86

Overdispersion. The Poisson residual deviance is 44× its df — massive overdispersion, so its SE on logED (0.013) is far too small. The quasipoisson rescales that SE by $\sqrt{51}\approx 7.1$ to 0.095; the negative binomial (a proper likelihood, $\theta = 0.73$) gives 0.081 and a residual-deviance ratio near 1 (0.94) — a well-calibrated fit. Trust fits 3 and 4, not the naive Poisson.

Back-transform — where the ×1.86 comes from. Don't memorize the recipe; derive it, because the same three steps turn any log-link slope into a sentence about biology. Start from what the model literally says. With a log link and a log-transformed predictor,

and exponentiating both sides puts it on the count scale:

Because we logged the predictor, eDNA enters as a power law with exponent $\beta_1$ — and that is what makes the effect multiplicative. Now ask the actual question: what happens to the expected count when the eDNA signal doubles, $\text{eDNA} \to 2\,\text{eDNA}$, with Year held fixed? Take the ratio of the new expected count to the old, so everything that didn't change cancels:

The intercept, the Year term, and even the baseline eDNA level all cancel — only the factor you multiplied the predictor by survives, raised to $\beta_1$. That is the general result, and it is worth keeping: multiply a log-transformed predictor by any factor $c$, and the expected response multiplies by $c^{\beta_1}$ (doubling is just $c = 2$; a 50% rise is $c = 1.5$, giving $1.5^{0.897}=1.44$, a 44% increase). It is the same fact seen from the link scale — doubling adds $\log 2$ to $\log(\text{eDNA})$, hence adds $\beta_1\log 2$ to the linear predictor, and $e^{\beta_1\log 2}=2^{\beta_1}$ back on the response scale. Additive on the link, multiplicative on the response.

Only now do the numbers enter. Plug in the NB estimate $\hat\beta_1 = 0.897$:

So a doubling of the flow-corrected eDNA signal is associated with roughly an 86% increase (×1.86) in the expected daily adult sockeye count — matching the graded Problem 5. "The slope was 0.90 (SE 0.08)" is not a finding; "doubling the eDNA signal is associated with 86% more sockeye" is the sentence that goes in the paper.

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 survival probability per class -- identical across the two forms
plogis(coef(f1)[1] + c(0, coef(f1)[-1]))    # crew, first, second, third
tapply(long$survived, long$class, mean)      # exactly the same numbers
# titanic_prop has 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)

Fitted probabilities per class come out crew 0.240, first 0.625, second 0.414, third 0.252 — identical to the raw proportions (a single categorical predictor gives a saturated model, so each group's MLE is just its empirical proportion).

(c) Why the two forms agree. Write $y_i\in\{0,1\}$ for each passenger and group them by class. For a group with $n$ passengers and $y=\sum y_i$ survivors, the two likelihoods are

They differ only by the factor $\binom{n}{y}$ — the number of ways to arrange $y$ survivors among $n$ people. Taking logs, the two log-likelihoods differ by the additive constant $\log\binom{n}{y}$, which does not contain $p$. Because the MLE and its standard errors depend only on how the log-likelihood changes with $p$ (its derivatives in $p$), a term free of $p$ drops out: the estimated coefficients, fitted probabilities, and SEs are identical. (The coefficient tables can still look different if the two data frames sort their factor levels differently and thus pick different reference classes — align the references and the numbers match exactly.)

What the constant does shift is the absolute log-likelihood, and therefore the AIC. Aggregating the Titanic passengers into the class×sex×age cells of titanic_prop.csv gives $\sum_{\text{groups}}\log\binom{n}{y} = 1025.44$, so the Binomial log-likelihood sits exactly that far above the Bernoulli one ($-268.84$ vs $-1294.28$) and the AICs differ by $2 \times 1025.44 = 2050.87$.

The residual deviance is a different story, and it is worth getting right. Deviance is $2(\ell_{\text{saturated}} - \ell_{\text{model}})$, and $\log\binom{n}{y}$ sits in both terms — so it cancels, and the deviance does not shift by that constant. The two deviances still differ (Bernoulli 2588.56 vs Binomial 491.06, a gap of 2097.49 — which is neither 1025.44 nor 2050.87), but for another reason entirely: the saturated reference model changes. The Bernoulli form's saturated model gets one parameter per passenger and so fits every row perfectly; the Binomial form's gets only one per cell. The gap is exactly the Bernoulli deviance of that cell-saturated model — the individual-level noise that aggregation simply stops measuring.

deviance(f1)                                                             # 2588.555  Bernoulli
deviance(f2)                                                             #  491.061  Binomial
deviance(glm(survived ~ class*sex*age, family = binomial, data = long))  # 2097.495  <- the gap

So never compare the deviance or AIC of a Bernoulli (one-row-per-individual) fit against an aggregated Binomial fit. The AIC gap is bookkeeping; the deviance gap is a change of yardstick. Neither is evidence about model fit. Within one data form, both are still perfectly good for comparing nested models.

library(effects); library(emmeans)
long <- read.csv("data/titanic_long.csv", stringsAsFactors = TRUE)
fit  <- glm(survived ~ class * sex + age, family = binomial, data = long)
plot(allEffects(fit))                                  # fitted probabilities
emmeans(fit, pairwise ~ class | sex, type = "response")  # class contrasts within each sex

(a) Largest pairwise odds ratio. The emmeans pairwise output (odds-ratio scale) gives, among females, first-vs-third class as the largest contrast:

# sex = female:
#  contrast        odds.ratio     SE  z.ratio  p.value
#  first / third        48.7    25.7    7.35   <.0001
# sex = male:
#  first / third         2.78    0.55   5.14    <.0001

First-class women had about 49× the odds of survival of third-class women — the sharpest class gradient in the model. Among men the same first-vs-third contrast is only ~2.8×: the "women and children first" protocol compressed the class effect for men.

(b) Why not TukeyHSD on an aov? TukeyHSD operates on an aov object, which assumes Normal errors, constant variance, and an identity link. A survival GLM has Bernoulli errors (variance $p(1-p)$, not constant) and a logit link, so Tukey on an aov would use the wrong likelihood, ignore the link, and ignore the interaction. emmeans on the fitted GLM forms the contrasts on the correct log-odds scale (optionally back-transformed to the response scale) with proper delta-method SEs and a Tukey multiplicity adjustment.

(c) Female 1st-vs-3rd, adult passenger.

predict(fit,
        newdata = data.frame(class = c("first","third"), sex = "female", age = "adult"),
        type = "response")
#     1st        3rd
#  0.972      0.419

An adult first-class woman is modeled at $p = 0.972$ versus $p = 0.419$ for an adult third-class woman — a probability difference of about 0.55 (55 percentage points). On the link scale that is a log-odds gap of $3.56 - (-0.33) = 3.89$, i.e. an odds ratio of $e^{3.89} \approx 49$, matching part (a). The odds ratio and the probability difference describe the same effect on two scales; report both, because a huge odds ratio near the probability ceiling still corresponds to a bounded (though here still large) probability gap.