Continuous distributions, linear models, and GLMs

Which distribution for a right-skewed body-size sample?

You record body weight (g) for every fish kept in a 30-day creel survey on a coastal river. Most fish are sub-legal-size cutthroat (small), with a long right tail of large hatchery steelhead retained in the same fishery. The empirical distribution is strictly positive and strongly right-skewed. (The same reasoning applies to tree DBH, seed mass, or any size variable built up by multiplicative growth.)

1. Which of {Normal, Lognormal, Gamma, Exponential, Beta} would you choose, and what is the data-generating story that justifies it in one sentence?
2. If $\log(\text{weight}) \sim \mathcal{N}(5.0,\,0.8^2)$, what is the expected weight on the raw scale? (Recall the lognormal mean is $e^{\mu + \sigma^2/2}$.)
3. What is the median weight? Why does it differ from the mean, and which is larger?
4. If you mistakenly fit a Normal to these weights and reported a 95% prediction interval centered at the mean, why would the lower endpoint be biological nonsense?

Length–weight allometry and the power law

For 200 Chinook salmon you have fork length (cm) and weight (g). Lengths range 30–110 cm; weights range 250 g to 18 kg. Two analysts fit:

• Analyst A: lm(weight ~ length)
• Analyst B: lm(log(weight) ~ log(length))

1. Describe what Analyst A's residuals-vs-fitted plot would look like, and name the assumption it violates.
2. Analyst B finds $\hat\beta_1 = 3.05$. Connect this slope to the standard allometric equation $W = a L^b$ and interpret $b$ biologically. Why is $b = 3$ (the "cube law") the natural null?
3. What multiplicative effect on weight corresponds to a 25% increase in length? To a doubling of length?
4. Analyst B back-transforms $\hat\beta_0 + \hat\beta_1 \log(L)$ by exponentiating to predict weight at a given length. What does that back-transformed value give you — the mean or the median weight — and what correction turns it into the other one?

Line-transect detection as an exponential process

On a line transect the probability of detecting an animal declines with perpendicular distance $x$. Model the detection process as Exponential with rate $\lambda = 0.02$ per meter, so that the chance of not yet having detected an animal at distance $x$ follows the survival function of the Exponential.

1. Write the survival function $S(x) = P(X \ge x)$ for an Exponential$(\lambda)$ and state its mean.
2. Compute the probability of missing an animal that is 50 m off the transect line.
3. At what perpendicular distance does the detection process reach its "half-life" (the distance at which $S(x) = 0.5$)? Show that this equals $\ln 2 / \lambda$.

Are the sockeye counts Poisson? Diagnose overdispersion

Load the adult sockeye counts and look at them before fitting any predictors. This problem is about the marginal count distribution.

d <- read.csv("data/sockeye_adult.csv", stringsAsFactors = FALSE)
d <- subset(d, Sockeyetype == "Sockeye_Adults")
d <- d[!is.na(d$Count), ]
mean(d$Count); var(d$Count)

1. Compute the mean and variance of Count. What is the dispersion ratio (variance / mean)? Does a Poisson — which forces variance = mean — look plausible?
2. Solve for the negative-binomial size parameter $k$ that reproduces the empirical variance under $\text{Var} = \mu + \mu^2/k$. (Use the empirical mean for $\mu$.)
3. Simulate 5000 draws from both rpois(mu) and rnbinom(mu, size = k) with your $\mu$ and $k$. Compare $P(\text{count} = 0)$ and the probability of a large "bumper run" (e.g. count > 400) between the two.
4. Name a biological mechanism that would generate this overdispersion in daily salmon counts, and say in one sentence what would go wrong if you reported Poisson-based confidence intervals for the mean daily count.

Sockeye eDNA → daily count: pick the right model and read the slope

Do adult sockeye numbers track the flow-corrected eDNA signal? Use the same subset as Problem 4, keep only rows where both the count and the eDNA are observed, and build a Year factor from the date.

d <- read.csv("data/sockeye_adult.csv", stringsAsFactors = FALSE)
d <- subset(d, Sockeyetype == "Sockeye_Adults")
d$Year <- factor(sub(".*/", "", d$Date))              # "15" or "16"
d <- d[!is.na(d$Count) & !is.na(d$Qcorr_qPCR), ]
min_pos <- min(d$Qcorr_qPCR[d$Qcorr_qPCR > 0])          # small constant so log() is finite
d$logED <- log(d$Qcorr_qPCR + min_pos)

Fit, in sequence, and compare:

1. lm(log(Count + 1) ~ logED + Year, data = d) — the old-school log-transform approach.
2. glm(Count ~ logED + Year, family = poisson, data = d).
3. glm(Count ~ logED + Year, family = quasipoisson, data = d).
4. MASS::glm.nb(Count ~ logED + Year, data = d).

1. For the Poisson fit, report residual deviance / residual df. Is overdispersion present? What does the quasipoisson dispersion estimate tell you, and how does the NB residual-deviance ratio compare?
2. Compare the standard error on the logED coefficient across the four fits. Which models give realistically wide confidence intervals and which give unrealistically narrow ones? Why?
3. Take the log-link slope on logED from your preferred model. Because the predictor is $\log(\text{eDNA})$, a doubling of the flow-corrected eDNA signal adds $\hat\beta \cdot \ln 2$ to the linear predictor, and because the link is the log, that additive shift multiplies the expected count by $2^{\hat\beta}$. Compute that multiplier and finish the sentence: "a doubling of the flow-corrected eDNA signal is associated with a ×___ change in the expected daily adult count."

Titanic survival: individual-level vs aggregated binomial

The same logistic regression can be fit from one-row-per-passenger data or from aggregated (successes, failures) counts. Show that they agree.

long <- read.csv("data/titanic_long.csv", stringsAsFactors = TRUE)
prop <- read.csv("data/titanic_prop.csv", stringsAsFactors = TRUE)
f1 <- glm(survived ~ class, family = binomial, data = long)
f2 <- glm(cbind(Yes, No) ~ Class, family = binomial, data = prop)

1. Fit both models. The fitted survival probability for each class should match between the two fits and should equal the raw survival proportion in that class. Verify this with predict(f1, type = "response") and tapply(long$survived, long$class, mean). Why must the individual-level (Bernoulli) and aggregated (Binomial) likelihoods give the same fitted probabilities?
2. The raw coefficient tables from f1 and f2 do not look identical. Explain why — think about which class each fit uses as its reference level. (Hint: levels(long$class) vs levels(prop$Class).)
3. Take the intercept and the first-class coefficient from f1 and convert them, by hand, from log-odds to odds to probability. Confirm your first-class probability against the fitted value and against the raw first-class survival proportion.
4. Express the first-class-vs-reference coefficient as an odds ratio and state, in plain English, how many times greater the odds of survival were for that class relative to the reference.

Titanic: does the class effect depend on sex?

Fit an interaction model and interpret it on the odds and probability scales.

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

1. From the emmeans pairwise output (on the odds-ratio scale), which class contrast within a sex is the largest? Report the odds ratio and say what it means in words.
2. The model above already includes age (a two-level factor, adult vs child) as a main effect. Express the agechild coefficient as an odds ratio and interpret it, holding class and sex fixed.
3. Using emmeans(fit, ~ class | sex, type = "response"), read off the fitted survival probability for first-class females and third-class females, and report the probability difference between them. Contrast this with what you would conclude from the odds ratio alone.
4. Why is a TukeyHSD on an aov object the wrong tool for these class contrasts, whereas emmeans on the fitted GLM is correct?

Reading a Poisson coefficient table as biology

From 120 ten-minute point counts in oak woodland you fit

glm(warblers ~ canopy_oak + elev_km, family = poisson, data = pc)

and obtain intercept $\hat\beta_0 = +0.10$, $\hat\beta_{\text{oak}} = +0.020$ per percent canopy oak, and $\hat\beta_{\text{elev}} = -0.55$ per km elevation.

1. What is the predicted mean count at oak = 0% and elev = 0 km?
2. What is the multiplicative effect on expected count of going from 10% to 60% oak canopy?
3. What is the elevation "half-decay" — the increase in elevation at which the expected count is halved? (Solve $e^{\hat\beta_{\text{elev}} \cdot \Delta} = 0.5$.)
4. Write the one-sentence biology summary you would put in a manuscript, with both effects on the multiplicative scale.

When the sign of an effect flips between groups

Macroinvertebrate counts per Surber sample from 90 collections are modeled as Poisson with predictors discharge (m³/s, continuous) and season (factor: spring reference, fall):

glm(macros ~ discharge * season, family = poisson, data = samps)

Coefficients: intercept = 3.50, discharge = $-0.30$, seasonFall = $+0.40$, discharge:seasonFall = $+0.50$.

1. What is the discharge slope in spring, on both the log scale and as a multiplicative (per-m³/s) effect?
2. What is the discharge slope in fall? Note that the two seasons respond in opposite directions — a finding you cannot get from a single pooled slope.
3. What is the expected count at discharge = 0 in the fall?
4. Solve for the discharge at which the spring and fall fitted counts are equal. Is that crossover inside the plausible range of discharges, and what does that imply about which season has higher counts across the observed gradient?

Turning counts into rates with an offset

Across 200 longline sets, the number of hooks varies from 500 to 5000. You record bycatch counts of yelloweye rockfish per set. Two analysts:

• Analyst A: glm(bycatch ~ depth_m, family = poisson)
• Analyst B: glm(bycatch ~ depth_m + offset(log(hooks)), family = poisson)

1. What does Analyst B's intercept represent that Analyst A's does not?
2. Why must the effort term enter as offset(log(hooks)) on the log scale, rather than by dividing the response by hooks?
3. If Analyst B's depth slope is $+0.012$ per m, write a one-sentence biological summary of the per-hook bycatch rate, and give the compounded multiplier over a 100 m depth change.
4. If shallow sets tended to use fewer hooks, in which direction would Analyst A's depth slope be biased, and why?

Occupancy thresholds and choosing a management lever

Presence/absence of hermit thrush at 300 sites is modeled as glm(present ~ basal_area + shrub_cover, family = binomial), giving intercept $-3.0$, basal-area slope $+0.08$ per m²/ha, and shrub slope $+0.04$ per percent shrub.

1. For a site with shrub cover = 25%, at what basal area does predicted occupancy probability equal 0.5? (Solve the linear predictor = 0.)
2. For a site with shrub cover = 50%, what is the new "threshold" basal area? Interpret the shift.
3. Compute the marginal effect of basal area on probability at the $p = 0.5$ inflection, using $\beta\,p(1-p)$. Confirm it against the divide-by-4 rule ($\beta/4$).
4. At a low-suitability site currently at $p = 0.3$, a manager can either raise basal area by 5 m²/ha or raise shrub cover by 10 percentage points. Using the marginal-effect approximation $\beta\,p(1-p)$ (or an exact plogis computation), which intervention yields the larger predicted gain in occupancy probability?

From a log-scale contrast to a percent change with uncertainty

Eelgrass above-ground biomass is regressed on sediment type (a factor with medium as the reference level) using a Gamma GLM with a log link. The "fine sediment vs medium sediment" contrast returns a coefficient $\hat\beta = -0.65$ with a 95% Wald CI of $(-0.95,\,-0.35)$ on the log scale.

1. Write out the full model explicitly in equations, not just prose: the response distribution, the link function, and the linear predictor with all terms. Define your indicator variable(s) and state the reference level.
2. Back-transform the point estimate. What multiplicative effect on biomass does it represent?
3. Back-transform the confidence interval. Is the resulting interval symmetric around the point estimate on the multiplicative scale?
4. Translate the point estimate and CI into a "percent decrease" statement.
5. Why is it correct to exponentiate each endpoint of the CI, but wrong in general to exponentiate the midpoint? Under what special condition does exponentiating the midpoint happen to land on the point estimate?