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 detection with an exponential detection function at rate $\lambda = 0.02$ per meter. Give each animal a latent detection range $X \sim \text{Exponential}(\lambda)$, the farthest perpendicular distance at which the observer would still have picked it up, so an animal at distance $x$ is seen exactly when $X \ge x$, and the detection function $g(x)$ is the survival function of $X$.

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. (Careful: $S(50)$ is the probability the animal is seen.)
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$.
4. The Exponential is memoryless. Explain what memorylessness implies about the detection hazard as perpendicular distance increases, and what that assumption commits you to biologically for an animal sitting almost on the trackline.

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.

What “no interaction” actually buys you

Problems 5–7 build one idea in three steps: an interaction means the effect of one predictor depends on the value of another. Everything you need is generated in the code block, so no data file is required. Work them in order.

Start with data built with no interaction, and deliberately fit a model that allows one anyway:

set.seed(123)
d5 <- data.frame(x1 = runif(200, -2, 2), x2 = runif(200, -2, 2))
d5$y <- 1 + 1.5*d5$x1 - 1*d5$x2 + rnorm(200, 0, 1)   # note: no x1*x2 term

m5 <- lm(y ~ x1 * x2, data = d5)    # x1 * x2 expands to x1 + x2 + x1:x2
summary(m5)

1. Report the four estimated coefficients against the values used to build the data (1, 1.5, −1, and no product term). What does the fitted x1:x2 coefficient come out as, and is it distinguishable from zero?
2. Write the full statistical model in mathematical notation, giving the distribution of $y_i$, its mean as a function of $x_{1i}$ and $x_{2i}$, and what each parameter is.
3. Differentiate the fitted mean with respect to $x_1$, and show your working. Why that is the right thing to differentiate is that the derivative is the change in $\hat y$ produced by a one-unit step in $x_1$, which is what anyone means by “the effect of $x_1$.” So asking whether $x_2$ turns up in it is asking whether that effect depends on where $x_2$ sits. Does $x_2$ appear anywhere in the result?
4. Using your fitted coefficients, compute the slope of $y$ on $x_1$ at $x_2 = -1$ and again at $x_2 = +1$. Report both. Why are they not exactly equal, given that the data were built with no interaction?
5. If you plotted the fitted surface over the $(x_1, x_2)$ plane, what shape would it be? Name it in one word, and say what the contour lines of equal $\hat y$ would look like.

When the effect of one predictor rides on another

Same design, same sample size, one change. The response is now built with a genuine product term.

set.seed(456)
d6 <- data.frame(x1 = runif(200, -2, 2), x2 = runif(200, -2, 2))
d6$y <- 2 + 1.5*d6$x1 - 1*d6$x2 + 2*d6$x1*d6$x2 + rnorm(200, 0, 1)

m6 <- lm(y ~ x1 * x2, data = d6)
summary(m6)

1. Report the four coefficients against the truth (2, 1.5, −1, 2).
2. Write out the full statistical model in mathematical notation, as in Problem 5(b).
3. Find the slope of the fitted mean in $x_1$ again. This time $x_2$ survives. Differentiate if you are comfortable doing that, but there is a route with no calculus in it. Gather the two terms containing $x_1$ and factor it out, so the fitted mean reads $\hat y = (\beta_0 + \beta_2 x_2) + (\beta_1 + \beta_3 x_2)\,x_1$. Whatever multiplies $x_1$ is the slope, and written that way you can read straight off it how the slope depends on $x_2$. Either way, write the slope as a formula in $x_2$, then evaluate it at $x_2 = -1$, $x_2 = 0$ and $x_2 = +1$. Something happens between the first and the last. Say what.
4. Solve for the value of $x_2$ at which the effect of $x_1$ vanishes entirely. Compare your fitted answer to the value implied by the parameters used to build the data.
5. Your collaborator writes: “the effect of $x_1$ on $y$ is 1.52.” Using your answer to (c), explain in two sentences why that sentence cannot be written for this model, and what it would have to say instead.
6. What shape is the fitted surface now, and what happened to the contour lines of equal $\hat y$ compared with Problem 5(e)?

The interaction you did not put in the model

Now the same structure with a binary response, which is where interactions stop being straightforward. Fit the model with the product term first:

set.seed(789)
d7 <- data.frame(x1 = runif(200, -2, 2), x2 = runif(200, -2, 2))
eta <- -1 + 1.2*d7$x1 - 0.8*d7$x2 + 1.5*d7$x1*d7$x2
d7$ybin <- rbinom(200, 1, plogis(eta))

m7 <- glm(ybin ~ x1 * x2, data = d7, family = binomial)
summary(m7)

1. Report the four coefficients. Which scale are they on, probability, odds, or log-odds?
2. Write the full model in mathematical notation, including the link function.
3. On the log-odds scale, the slope in $x_1$ is the same kind of formula you found in Problem 6(c). Write it down and evaluate it at $x_2 = -1$ and $x_2 = +1$.
4. Now refit the same data with no interaction term at all:

   m7add <- glm(ybin ~ x1 + x2, data = d7, family = binomial)   # additive on the log-odds scale
   coef(m7add)

   This model is additive by construction. On the log-odds scale, the effect of $x_1$ is one number that does not involve $x_2$. On the probability scale it is not, and here is the formula you need, taken as given: $$\frac{\partial p}{\partial x_1} = p(1-p)\,(\beta_1 + \beta_3 x_2).$$ The $p(1-p)$ factor is the logistic curve's own steepness, biggest at $p = 0.5$ and shrinking toward zero at either end, so the same step on the log-odds scale moves the probability by different amounts depending on where you started. Using that formula with $\beta_3 = 0$, compute $\partial p/\partial x_1$ at $x_1 = 0$ for $x_2 = -1$ and for $x_2 = +1$. Report both, and their ratio.
5. You just found that the effect of $x_1$ on the probability depends on $x_2$ in a model containing no interaction term. Explain in three or four sentences why this happens, and what it implies for a reader who wants to know whether two predictors “interact” in a GLM. Your answer should distinguish clearly between the two scales.

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, giving 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?