Continuous distributions, linear models, and GLMs

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

You record body weight (g) for every rainbow trout kept in a 30-day creel survey on a mountain lake. A histogram of the weights has almost no fish near zero. It rises to a single peak at roughly 80 g and then trails off in a long right tail of a few much larger fish approaching a kilogram. The distribution is strictly positive and strongly right-skewed, but its most common value sits well above zero. (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? Which candidate does the peak above zero rule out, and why?
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. Read the note below first. Does that back-transformed value estimate the mean or the median weight at that length? Analyst B's residual standard deviation on the log scale is $\hat\sigma = 0.12$. By what percentage does the mean weight exceed the plain back-transformed value? Compare that with the trout in Problem 1, where $\sigma = 0.8$, and say what makes the correction matter so much more there.

Bat passes at an acoustic detector

An ultrasonic detector on a forest edge logs every bat that flies past. Through a summer night, passes arrive at random at an average rate of $\lambda = 4$ per hour, each independent of the others. That one assumption describes the same process in two ways.

• The number of passes in a window of $t$ hours is Poisson with mean $\lambda t$.
• The waiting time $T$ until the next pass is Exponential with rate $\lambda$, with density $f(t) = \lambda e^{-\lambda t}$ and cumulative distribution function $F(t) = P(T \le t) = 1 - e^{-\lambda t}$ for $t \ge 0$.

Work in hours, converting minutes as you go, so 30 minutes is $t = 0.5$.

1. What are the mean and the median waiting time between passes, in minutes? (The median is the $t$ at which $F(t) = 0.5$.) Which is larger, and what about the shape of the Exponential makes it so?
2. Find the probability that 30 minutes go by with no bat at all, in two ways. First use the Poisson, as the probability of zero passes in half an hour. Then use the Exponential, as the probability that the wait for the next pass is longer than half an hour. Explain in one sentence why the two answers have to be identical.
3. The detector has already been silent for 20 minutes. What is the probability that it stays silent for at least another 30 minutes, that is, $P(T > 50 \text{ min} \mid T > 20 \text{ min})$? Compare your answer with (b). This is memorylessness: under a constant rate, having already waited 20 minutes tells you nothing about how much longer you will wait. In one or two sentences, say what memorylessness requires of the bats, and give one biological reason real bat passes might not be memoryless.
4. Check the connection by simulation. Generate waiting times, turn them into the clock time of each pass, and count the passes in each hour:

   set.seed(11)
   gaps   <- rexp(5000, rate = 4)                       # waiting times between passes, in hours
   times  <- cumsum(gaps)                               # clock time of each pass
   times  <- times[times < 1000]                        # keep the first 1000 hours
   counts <- tabulate(floor(times) + 1, nbins = 1000)   # passes in each hour
   mean(counts); var(counts)

   Report the mean and variance of the hourly counts. Which distribution do those two numbers point to, and why?

Are the sockeye counts Poisson? Diagnose overdispersion

Load the adult sockeye counts and look at them before fitting any predictors. Parts (a)–(d) are about the marginal count distribution: every day’s count pooled together, with nothing in the model to explain why one day differs from another. Part (e) asks what changes once a predictor enters the picture.

d <- read.csv(here::here("Day2_GLM", "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.
5. Overdispersed counts do not always mean a Poisson regression is wrong. A salmon run builds, peaks, and fades over the season, so different days have very different expected counts. Simulate a season in which every day’s count is exactly Poisson given that day’s expected count, with the expected count following a run-timing curve that peaks at about 400 fish:

   set.seed(2026)
   day <- 1:90
   mu  <- exp(6 - ((day - 45) / 20)^2)   # expected count on each day
   y   <- rpois(90, mu)                  # each day's count is exactly Poisson
   var(y) / mean(y)                      # dispersion ratio of the pooled counts
   
   m <- glm(y ~ day + I(day^2), family = poisson)
   sum(residuals(m, type = "pearson")^2) / df.residual(m)   # dispersion after accounting for day

   Report both dispersion numbers. The first compares the variance of the pooled counts with their mean, exactly as in (a). The second makes the same comparison after the model has accounted for day, using how far each count falls from its own day’s fitted mean; a value near 1 means the counts scatter around their fitted means just as much as a Poisson says they should, and no more. Why do the pooled counts look heavily overdispersed when every single count is Poisson? (Hint: where does the variation in y come from when different days have different expected counts?) What does this mean for the conclusion you drew about the sockeye counts in (a)?

Simulating with and without an interaction

Problem 5 builds one idea: an interaction means the effect of one predictor depends on the value of another. Everything you need is generated in the code blocks, so no data file is required.

Simulate two datasets with the same design and sample size. The first is built with no interaction and the second with a genuine product term. Then fit a model that allows an interaction to each one.

set.seed(123)
d_none <- data.frame(x1 = runif(200, -2, 2), x2 = runif(200, -2, 2))
d_none$y <- 1 + 1.5*d_none$x1 - 1*d_none$x2 + rnorm(200, 0, 1)   # no x1*x2 term
m_none <- lm(y ~ x1 * x2, data = d_none)    # x1 * x2 expands to x1 + x2 + x1:x2
summary(m_none)

set.seed(456)
d_int <- data.frame(x1 = runif(200, -2, 2), x2 = runif(200, -2, 2))
d_int$y <- 2 + 1.5*d_int$x1 - 1*d_int$x2 + 2*d_int$x1*d_int$x2 + rnorm(200, 0, 1)   # a genuine x1*x2 term
m_int <- lm(y ~ x1 * x2, data = d_int)
summary(m_int)

1. Report the four estimated coefficients for each model against the values used to build its data, which are 1, 1.5, −1 and no product term for m_none, and 2, 1.5, −1 and 2 for m_int. In each model, is the fitted x1:x2 coefficient 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. Explain why this one model statement describes both fits.
3. Find the slope of the fitted mean in $x_1$. That slope is the change in $\hat y$ produced by a one-unit step in $x_1$, which is exactly what anyone means by “the effect of $x_1$.” You can differentiate the fitted mean with respect to $x_1$. Or, with no calculus, gather the two terms that contain $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. Write the slope as a formula in $x_2$, then evaluate it with each model’s fitted coefficients, for m_none at $x_2 = -1$ and $x_2 = +1$, and for m_int at $x_2 = -1$, $0$ and $+1$. Why are the two m_none slopes not exactly equal, given that its data were built with no interaction? And what happens to the m_int slope between $x_2 = -1$ and $x_2 = +1$?
4. For m_int, 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 about m_int: “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. Picture each fitted mean as a surface over the $(x_1, x_2)$ plane. First suppose the interaction coefficient were exactly zero, as it is in the equation that generated d_none, so that $\hat y = \hat\beta_0 + \hat\beta_1 x_1 + \hat\beta_2 x_2$. What shape is that surface, in one word, and what do its contour lines of equal $\hat y$ look like? Then describe the surface for m_int, and say what happened to its contour lines.

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? Work on the link scale. A probability of 0.5 means the odds are $0.5/(1-0.5) = 1$, and $\log(1) = 0$, so $p = 0.5$ exactly when the log-odds equal 0. The model says the log-odds equal $-3.0 + 0.08 \times \text{basal area} + 0.04 \times \text{shrub cover}$, so set that expression equal to 0 and solve for basal area.
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?