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, and 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 is that 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, and 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.

"Not seen by distance $x$" means the distance you had to travel exceeded $x$, so we want $P(X \ge x)$, the survival function. There are two routes to it, and it is worth seeing both.

Way 1: from first principles (integrate the density). Accumulate all the probability beyond $x$ by integrating the pdf from $x$ outward:

Way 2: from the CDF you already know (no integral). Once you have named the distribution, you can skip the calculus. The Exponential CDF is a standard result, $F(x) = P(X \le x) = 1 - e^{-\lambda x}$, and what we want is its complement:

Same answer, and Way 2 takes one line. That is the payoff for recognising the distribution instead of re-deriving it; Way 1 is really just the on-the-spot derivation of that CDF.

The general result to carry away is that for any exponential process the survival function is a pure exponential decay, $P(X \ge x) = e^{-\lambda x}$. The fraction still unseen falls by the constant factor $e^{-\lambda}$ per extra meter, so there is a "half-distance" $\ln(2)/\lambda$ at which half the orangutans remain unseen, exactly like a half-life.

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)$.

(a) The Beta has mean $\dfrac{\alpha}{\alpha+\beta} = \dfrac{0.03\phi}{0.03\phi + 0.97\phi} = 0.03$: the $\phi$ cancels, so every choice of $\phi$ leaves the mean at exactly the 0.03 of the plain assumption. $\phi$ sets only the spread. Small $\phi$ makes cameras wildly unequal, large $\phi$ makes them nearly identical. That is why the two simulations below share a mean of ~6 and differ only in their tails.

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, the same mean with 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. That is 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, and 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, since 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. "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).

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.