14. Poisson regression and the GLM idea

In 1973 two biologists, Michael Johnson and Peter Raven, published a count. For each of 30 islands in the Galapagos archipelago they recorded how many plant species grew there. Alongside each count they noted a handful of geographic facts: the island’s area, its highest elevation, its distance to the nearest neighbor, its distance to the big island of Santa Cruz, and the area of the closest adjacent island. The question was old and simple. Why do some islands hold hundreds of species and others only a dozen?

You have met this dataset twice already. In Chapter 9 it was a diagnostics cautionary tale. Fit an ordinary least-squares line, and one island, Isabela, the largest in the chain, sits so far out in predictor space that it bends the whole fit toward itself. In Chapter 10 you tried the standard remedy, transforming the response, and it helped the residual plots but never quite felt right. Something about the data resisted the linear model, and we promised you an honest answer later. Later is now.

The trouble is visible in one picture. Figure 1 plots species against island area with the least-squares line drawn through it. The fit does something impossible. Its prediction band has a single fixed width, so pushed toward small islands the band dips below zero and promises negative species counts. And one fixed width cannot fit both ends: small islands cluster tightly near zero while large islands scatter across hundreds. A species count cannot be negative, cannot be fractional, and its spread grows with its size. The linear model assumes none of that. This chapter builds the model that does.

By the end of this chapter the species counts will have a model that respects what they are, the same model will handle ship-damage rates with wildly different exposures, and all of it, together with the straight lines of Chapter 2 and the logistic curves of Chapter 13, will collapse into a single family with one estimating idea. That family is the generalized linear model, and it is the payoff of this whole part of the course.

14.1 Why counts break the linear model¶

Intuition¶

A count is a whole number of things that happened: species on an island, damage incidents on a ship, customers in an hour, seizures in a week. Counts have three stubborn features. They are never negative. They are whole numbers, not fractions. And their variability grows with their average: an island that averages 5 species varies by a species or two, while an island that averages 200 varies by dozens. The linear model of Chapters 2 through 8 assumes a response that can wander anywhere on the number line, comes in any real value, and scatters by the same $\sigma^2$ at every level of the predictors. Every one of those assumptions is false for a count.

The natural probability model for a count is the Poisson distribution (Definition 14.1). It describes the number of events that occur when each event is rare, events happen independently, and they arrive at a steady average rate. A Poisson distribution has a single parameter, its mean $\mu$, and a notable property: its variance equals its mean. Figure 2 shows three Poisson distributions. As the mean climbs from 1 to 10, the distribution slides right, spreads out, and grows more symmetric, and at every step its spread is locked to its center.

That mean-equals-variance link is exactly the pattern Figure 1 showed and exactly what constant-variance least squares cannot represent. Figure 3 makes the contrast explicit: the linear model draws a flat variance line, while count data ride the diagonal where variance tracks the mean.

Formula¶

A random count $Y$ follows a Poisson distribution with mean $\mu > 0$ when

• $Y$ is the count, a nonnegative whole number.

• $\mu$ is the mean of the distribution, the average count.

• $y!$ is the factorial of $y$, which normalizes the probabilities to sum to one.

In words: the probability of seeing exactly $y$ events is set entirely by the average rate $\mu$. The distribution carries two facts we will use constantly,

In words: for a Poisson count the mean and the variance are the same number. This single equation is what separates count regression from everything before it, and it is the thread we follow to the end of the chapter.

Before we build a model for counts, spend a minute with the distribution itself, because the slider below is the only parameter a Poisson has.

14.2 The Poisson log-linear model and rate ratios¶

Intuition¶

We want the mean count $\mu_i$ to depend on predictors, the same ambition as ordinary regression, but with two repairs. The mean must stay positive no matter what the predictors do, and the model should describe how predictors multiply the count rather than add to it, because count processes are usually multiplicative: doubling an island’s area does not add a fixed number of species, it scales the expected number.

Both repairs come from one move: model the logarithm of the mean as a linear function of the predictors. Because the logarithm can be any real number while its inverse, the exponential, is always positive, a linear model on the log scale can range freely while the mean it produces stays positive. And because a sum on the log scale is a product on the original scale, coefficients become multipliers. Figure 5 shows the link at work: a straight line in the linear predictor $\eta$ becomes a curve in the mean $\mu$ that hugs zero for negative $\eta$ and climbs steeply for positive $\eta$, never once dipping below zero.

Formula¶

The Poisson log-linear regression model says that for each case $i$,

• $Y_i$ is the observed count for case $i$, assumed Poisson given the predictors.

• $\mu_i = E\{Y_i\}$ is its mean.

• $\eta_i$ is the linear predictor, the familiar straight-line combination of predictors.

• $\log$ is the link function here: it maps the positive mean onto the whole real line where the linear predictor lives.

• $\beta_0, \dots, \beta_{p-1}$ are the regression coefficients, estimated by $b_0, \dots, b_{p-1}$.

Solving the link for the mean gives the multiplicative form

In words: the mean count is a product of factors, one per predictor. This is why coefficients read as rate ratios. Increase $X_{ij}$ by one unit, holding the others fixed, and $\eta_i$ rises by $\beta_j$, so $\mu_i$ is multiplied by $e^{\beta_j}$:

In words: $e^{\beta_j}$ is the factor by which the expected count changes for a one-unit rise in $X_j$. A coefficient of zero means a ratio of one, no effect; a positive coefficient means a ratio above one, the count grows; a negative coefficient means a ratio below one, the count shrinks. This is the count-data cousin of the odds ratio $e^{\beta_j}$ you met for logistic regression in 13.3 Reading coefficients as odds ratios, and of the percent-change reading of a log-scale coefficient from 10.2 The log transformation and reading its coefficients.

Here is the picture to keep in your head. On the log scale the model is an ordinary straight line, so equal steps in a predictor add the same amount $\beta_j$ each time. Exponentiate, and that even adding turns into even multiplying: each step multiplies the count by the same factor $e^{\beta_j}$. Figure 6 shows both sides at once, using the elevation effect from the fit below, where one extra hundred meters multiplies the expected species count by about 1.42.

Derivation (the Poisson log-likelihood and its score equations)¶

Least squares will not estimate this model: the response is not normal and the mean is not linear. As in logistic regression, we estimate by maximum likelihood, the principle from 13.2 Fitting by maximum likelihood that picks the coefficients making the observed counts most probable.

Proof. The counts $Y_1, \dots, Y_n$ are independent given the predictors, so the likelihood is the product of their Poisson probabilities, and the log-likelihood is the sum of the logs:

Substitute the model $\log \mu_i = \mathbf{x}_i'\boldsymbol\beta$, so that $\mu_i = e^{\mathbf{x}_i'\boldsymbol\beta}$, where $\mathbf{x}_i = (1, X_{i1}, \dots, X_{i,p-1})'$ is the predictor vector for case $i$:

To maximize, differentiate with respect to the coefficient vector and set the result to zero. The last term does not involve $\boldsymbol\beta$. For the first two, using $\partial (\mathbf{x}_i'\boldsymbol\beta)/\partial \boldsymbol\beta = \mathbf{x}_i$ and $\partial e^{\mathbf{x}_i'\boldsymbol\beta}/\partial \boldsymbol\beta = e^{\mathbf{x}_i'\boldsymbol\beta}\mathbf{x}_i = \mu_i \mathbf{x}_i$,

These are the score equations for Poisson regression. Read the first component (the one for the intercept, where $\mathbf{x}_i$ has a 1): it says $\sum_i (y_i - \hat\mu_i) = 0$, the fitted means reproduce the total count. Each other component says $\sum_i (y_i - \hat\mu_i) X_{ij} = 0$, the residuals $y_i - \hat\mu_i$ are orthogonal to every predictor. These are the analogues of the normal equations of 2.2 Least squares from first principles, with the fitted mean $\hat\mu_i$ playing the role of the fitted value. The difference is that $\hat\mu_i = e^{\mathbf{x}_i'\mathbf{b}}$ is nonlinear in the coefficients, so unlike least squares there is no closed-form solution. The equations are solved by the same iteratively reweighted least squares that fits logistic regression, sketched in 13.2 Fitting by maximum likelihood; software runs it in a handful of steps. $\blacksquare$

Notice how little changed from logistic regression. There the score equations were $\sum_i (y_i - \hat\pi_i)\mathbf{x}_i = \mathbf{0}$ with $\hat\pi_i$ a probability; here they are $\sum_i (y_i - \hat\mu_i)\mathbf{x}_i = \mathbf{0}$ with $\hat\mu_i$ a mean count. Same shape, different mean model. Section 14.6 shows this is no accident.

R¶

The fitting function is glm, the same one used for logistic regression, with family = poisson. The default link for the Poisson family is the log, so you rarely name it.

gala <- read.csv("data/gala.csv")
pois_fit <- glm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent,
                family = poisson, data = gala)
summary(pois_fit)

Call:
glm(formula = Species ~ Area + Elevation + Nearest + Scruz +
    Adjacent, family = poisson, data = gala)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)  3.155e+00  5.175e-02  60.963  < 2e-16 ***
Area        -5.799e-04  2.627e-05 -22.074  < 2e-16 ***
Elevation    3.541e-03  8.741e-05  40.507  < 2e-16 ***
Nearest      8.826e-03  1.821e-03   4.846 1.26e-06 ***
Scruz       -5.709e-03  6.256e-04  -9.126  < 2e-16 ***
Adjacent    -6.630e-04  2.933e-05 -22.608  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 3510.73  on 29  degrees of freedom
Residual deviance:  716.85  on 24  degrees of freedom
AIC: 889.68

Number of Fisher Scoring iterations: 5

Steering the log link by hand is the fastest way to feel what maximum likelihood is doing, so take the two coefficients yourself and let the log-likelihood tell you when you are close.

14.3 Offsets: modeling a rate instead of a count¶

Intuition¶

Raw counts are only comparable when they come from equal opportunity. A cargo ship at sea for 40,000 months has far more chance to suffer wave damage than one at sea for 100 months, so comparing their raw incident counts is meaningless. What we want is the rate: incidents per month of service. The exposure, months of service, differs from ship to ship and must be built into the model, not ignored.

The fix is elegant. Model the rate as log-linear, then rewrite it as a model for the count with one fixed piece added to the linear predictor. That fixed piece, the log of the exposure, is called an offset (Definition 14.4): a predictor whose coefficient is not estimated but pinned to one, because we know exactly how the count should scale with exposure. Double the months at sea and, all else equal, you expect double the incidents. Figure 8 shows why raw counts mislead: ship type B has by far the most incidents, but also by far the most service; per month of service its damage rate is ordinary.

Formula¶

Let $t_i$ be the exposure for case $i$ (months of service, person-years, area surveyed) and let $\lambda_i$ be the event rate per unit exposure. The count’s mean is rate times exposure, $\mu_i = \lambda_i t_i$. Model the log rate as linear:

Derivation (the offset). Since $\mu_i = \lambda_i t_i$, take logs:

In words: modeling the log rate is the same as modeling the log count with $\log t_i$ added to the linear predictor as a term with known coefficient exactly one. That extra term is the offset. It shifts each case’s linear predictor by its log-exposure without costing a parameter, so the estimated coefficients describe the rate $\lambda_i$, and $e^{\beta_j}$ is now a rate-ratio per unit exposure. If you wrongly dropped the offset, the model would explain incidents partly through whatever predictor happens to correlate with service time, confounding rate with exposure. $\blacksquare$

R¶

Ships with zero recorded service contribute no exposure and no opportunity for incidents, so $\log t_i$ is undefined for them; we drop those rows first. The offset goes in through the offset argument.

ships <- read.csv("data/ships.csv")
ships <- ships[ships$service > 0, ]
ships$year   <- factor(ships$year)
ships$period <- factor(ships$period)
ship_fit <- glm(incidents ~ type + year + period,
                family = poisson, offset = log(service), data = ships)
round(summary(ship_fit)$coefficients, 4)

            Estimate Std. Error  z value Pr(>|z|)
(Intercept)  -6.4059     0.2174 -29.4600   0.0000
typeB        -0.5433     0.1776  -3.0595   0.0022
typeC        -0.6874     0.3290  -2.0891   0.0367
typeD        -0.0760     0.2906  -0.2614   0.7938
typeE         0.3256     0.2359   1.3803   0.1675
year65        0.6971     0.1496   4.6587   0.0000
year70        0.8184     0.1698   4.8207   0.0000
year75        0.4534     0.2332   1.9446   0.0518
period75      0.3845     0.1183   3.2507   0.0012

An offset is easier to believe once you have tried the alternatives, so put the exposure exponent under your finger and see what each setting is claiming about the ships.

14.4 Overdispersion and the quasi-Poisson fix¶

Intuition¶

The Poisson model makes a strong promise: the variance equals the mean. Real count data often break that promise by scattering more than a Poisson allows, a condition called overdispersion (Definition 14.5). Its causes are ordinary: predictors you did not measure, clustering of events, individual-to-individual differences in the underlying rate. Whatever the cause, the symptom is the same, the counts spread wider than their fitted means.

Overdispersion rarely changes the coefficient estimates much, but it wrecks their standard errors. A Poisson fit that assumes variance equals mean, when the true variance is several times larger, reports standard errors several times too small, turning noise into false significance. So overdispersion is not a nuisance to shrug off; it is the difference between an honest analysis and a misleading one. Figure 10 shows the diagnosis for the Galapagos fit: the squared residuals, which should hover around the fitted mean under a true Poisson, sit far above the variance-equals-mean line.

Formula¶

Measure dispersion with the Pearson statistic (Definition 14.6) divided by its degrees of freedom. Define the Pearson residual for case $i$ and the estimated dispersion $\hat\phi$:

• $r_i^{P}$ standardizes each residual by the Poisson standard deviation $\sqrt{\hat\mu_i}$.

• $X^2 = \sum_i (r_i^{P})^2$ is the Pearson chi-square statistic.

• $n - p$ is the residual degrees of freedom, cases minus estimated parameters.

In words: $\hat\phi$ is the average squared standardized residual. Under a correct Poisson model it should be near 1. A value well above 1 signals overdispersion; a value of $\hat\phi$ around 2 means the true variance is roughly twice the mean, and so on.

The quasi-Poisson fix keeps the Poisson mean model and its coefficient estimates untouched, because the score equations $\sum_i (y_i - \mu_i)\mathbf{x}_i = \mathbf{0}$ do not involve the variance. It only rescales the variance, assuming

for an unknown constant $\phi > 1$, estimated by $\hat\phi$. Every standard error is then multiplied by $\sqrt{\hat\phi}$:

In words: inflate each Poisson standard error by the square root of the estimated dispersion. This is the same idea you saw in 10.4 Weighted least squares, where modeling nonconstant variance changed the weights and the standard errors while leaving the story of the mean intact. Here the variance is allowed to be a constant multiple of the mean rather than exactly the mean.

R¶

glm reports the raw dispersion assumption, so compute $\hat\phi$ directly from the Pearson residuals, then refit with family = quasipoisson, which does the rescaling for you.

phi_hat <- sum(residuals(pois_fit, type = "pearson")^2) / df.residual(pois_fit)
phi_hat

[1] 31.74914

Figure 11 drives the point home across all five predictors: dividing each test statistic by $\sqrt{\hat\phi}$ pulls two of them below the rough significance threshold that the naive Poisson had them clearing with room to spare.

A nod to the negative binomial¶

Quasi-Poisson is a quick, assumption-light patch: it does not name a probability distribution, it just rescales variances. When you want a genuine likelihood, and therefore AIC, likelihood-ratio tests, and honest prediction intervals, the standard choice is the negative binomial model. It adds one parameter $\theta$ and lets the variance grow quadratically, $\operatorname{Var}\{Y\} = \mu + \mu^2/\theta$, which fits heavy overdispersion better than the strictly linear quasi-Poisson variance. In R it is MASS::glm.nb; the Galapagos data give $\hat\theta = 1.67$, small enough to confirm severe overdispersion. For this course, diagnose with $\hat\phi$, reach for quasi-Poisson as the everyday fix, and know the negative binomial is the next step when you need a full likelihood.

Dispersion is just a number until you watch it move a p-value, so drag it yourself and see which Galapagos predictors survive the correction.

14.5 The Galapagos thread, resolved¶

We can now close the story that has run since Chapter 9. Recall the two earlier visits. Chapter 9 made the Galapagos data its lesson in influence: Isabela, the largest island, is the high-leverage point that 9.1 Leverage values and the hat matrix points to, one an ordinary least-squares fit cannot help but chase, because it sits alone at the far edge of the predictor space. In 10.2 The log transformation and reading its coefficients and 10.3 Box-Cox: letting the data choose the power you tried the classic remedy, transforming the response, taking logs or a Box-Cox power to tame the skew and the widening scatter. It improved the residual plots, but it never removed the discomfort, because the real problem was never the scale of a continuous response. The response was a count, and no transformation of a count turns it into a normal, constant-variance measurement.

The Poisson model fixes the actual defect. It builds in the nonnegativity, the multiplicative structure, and the mean-variance link that species counts genuinely have, and once we account for overdispersion its standard errors are honest. Figure 13 puts the two philosophies side by side: the log-transform-then-least-squares curve of Chapter 10, and the quasi-Poisson mean curve of this chapter, both against the raw counts on a log-area axis. They agree in the middle and part company at the extremes, where the count model keeps predictions positive and lets the scatter widen as it should.

The lesson is bigger than one archipelago. Chapter 10’s transformations are the right tool when a continuous response is on an awkward scale. When the response is a fundamentally different kind of thing, a count, a yes or no, a rate, you do not bend the response to fit the linear model; you change the model. That change is the subject of the last section, and it is the reason logistic and Poisson regression belong in the same course as least squares.

So the first decision in any analysis is not which predictors to use; it is what kind of thing the response is. Figure 14 turns that decision into a short flowchart: name the response type, follow the arrow, and you land on the matching model. The next section shows why all four destinations are really one machine.

14.6 One family: the generalized linear model¶

Intuition¶

Step back and look at what the last three chapters share. Ordinary regression models a continuous mean as a straight line. Logistic regression models a probability by passing a straight line through the logit link. Poisson regression models a count by passing a straight line through the log link. In every case there is a linear predictor $\eta = \mathbf{x}'\boldsymbol\beta$ doing the same work; what changes is only two choices: which probability distribution the response follows, and which link function connects its mean to the linear predictor. Fix those two choices and everything else, the estimation by maximum likelihood, the deviance, the workflow, is shared machinery. That is the generalized linear model, or GLM (Definition 14.8).

Formula¶

A generalized linear model has three parts:

• The random component names the response distribution (normal, binomial, Poisson, and others), which fixes how the variance depends on the mean.

• The linear predictor $\eta_i$ is the same linear combination of predictors used all book long.

• The link function $g$ maps the mean $\mu_i$ onto the scale of the linear predictor; its inverse maps back, $\mu_i = g^{-1}(\eta_i)$.

In words: choose a distribution and a link, and you have specified a model; the same estimation engine fits them all. The three models of this part are three rows of one table.

Figure 15 draws the same idea as a diagram: one linear predictor feeding three links into three response types.

R¶

The unification is not just conceptual; it is literally how the software is built. The function glm fits every row of the table, and choosing family = gaussian(link = "identity") reproduces ordinary least squares exactly.

14.7 Chapter summary¶

You can now model counts. This chapter began with why ordinary least squares is the wrong tool for a count, whose nonnegative, whole-number, mean-equals-variance behavior violates the unbounded constant-variance response the linear model assumes. The repair was the Poisson log-linear model, which models the log of the mean as a linear predictor so the fitted mean stays positive and every coefficient reads as a rate ratio $e^{\beta_j}$. Maximizing the Poisson log-likelihood gave score equations of exactly the same form as the normal equations of Chapter 2. An offset, the log of the exposure entered with fixed coefficient one, turned the count model into a rate model for the ship data. The Pearson dispersion $\hat\phi$ diagnosed the severe overdispersion in the Galapagos counts, and the quasi-Poisson correction inflated the standard errors to honest sizes while leaving the estimates alone. That finally closed the Galapagos thread from Chapters 9 and 10: the trouble was never the scale of the response, it was that the response is a count. And all of it, together with the linear and logistic models, is three rows of one generalized linear model table.

Key results at a glance

Key terms. Poisson distribution, Poisson log-linear regression model, linear predictor, link function, rate ratio, score equations, offset, exposure, overdispersion, Pearson dispersion, quasi-Poisson, negative binomial, generalized linear model.

You should now be able to

• Explain why ordinary least squares is the wrong tool for count responses, naming the two model assumptions it violates.

• State the Poisson log-linear model and interpret its coefficients as rate ratios, in the units of the problem.

• Derive the Poisson log-likelihood and its score equations, and connect them to the maximum-likelihood machinery of Chapter 13.

• Use an offset to model a rate, and explain why the log of exposure enters the linear predictor with a fixed coefficient of one.

• Diagnose overdispersion from the Pearson dispersion and apply the quasi-Poisson correction, and say when a negative binomial model is the better fix.

• Resolve the Galapagos species question with a count model, and explain why it succeeds where the transformations of Chapter 10 fell short.

• Assemble the one-table GLM framework (family plus link) that contains linear, logistic, and Poisson regression as instances of a single idea.

Where this fits. In the workflow spine of The modeling workflow this chapter lives mostly in FIT and CHECK, with a strong USE payoff. The Poisson log-linear model is a FIT step for a response that the earlier chapters’ FIT tools could not handle; the overdispersion diagnosis is a CHECK step, the count-data analogue of the residual-variance checks of Chapter 9 and the variance modeling of 10.4 Weighted least squares; and reading coefficients as rate ratios is the USE step, the count cousin of the odds ratios of 13.3 Reading coefficients as odds ratios. Looking back, this chapter completes the arc that began when Chapter 13 first broke free of the normal response: logistic regression handled yes-or-no, Poisson regression handles how-many, and the GLM table shows both as extensions of the least squares you have known since Chapter 2. Looking forward, the same “choose a family and a link” habit carries into gamma models for skewed positive data, multinomial models for several categories, and the mixed and hierarchical models that a second course in regression takes up.

14.8 Frequently asked questions¶

Q1. Can’t I just log-transform the counts and run ordinary least squares? You can, and people did for decades, but it has real costs. Counts of zero make $\log 0$ undefined, forcing awkward “add 0.5” patches. The back-transformed predictions describe a median, not a mean, so they do not add up correctly. And the constant-variance assumption is still wrong even after logging. Poisson regression sidesteps all three by modeling the count directly on its own terms.

Q2. Why the log link specifically? Could I use a different one? The log link is the natural choice because it makes coefficients multiplicative (rate ratios), keeps the mean positive, and is the “canonical” link that makes the Poisson math cleanest. You can use others (the identity link gives an additive count model, occasionally useful), but the log link is the default for good reasons and is what “Poisson regression” means unless stated otherwise.

Q3. My dispersion is $\hat\phi = 0.8$, below 1. Is that underdispersion, and do I fix it? Values modestly below 1 are common and usually harmless; they make Poisson standard errors slightly conservative (too large), so your tests are if anything too cautious. Genuine underdispersion (from, say, counts that are more regular than random, like eggs in a full carton) exists but is rare. Worry about $\hat\phi$ well above 1, not modestly below it.

Q4. Does quasi-Poisson change my coefficient estimates? No. The coefficients solve the score equations $\sum_i (y_i - \mu_i)\mathbf{x}_i = \mathbf{0}$, which never mention the variance, so quasi-Poisson leaves them identical to the Poisson fit. Only the standard errors, and therefore the p-values and confidence intervals, change. That is the whole point: the estimates were fine, the uncertainty was understated.

Q5. When should I use the negative binomial instead of quasi-Poisson? Use quasi-Poisson for a fast, distribution-free correction to standard errors. Reach for the negative binomial when you need a genuine likelihood: to compare models by AIC, to run likelihood-ratio tests, or to form prediction intervals for new counts. The negative binomial also fits very heavy overdispersion better because its variance grows quadratically in the mean.

Q6. Why did dropping the zero-service ships not bias the analysis? A ship with zero months of service had no opportunity to suffer damage, so it carries no information about the damage rate, and $\log 0$ is undefined anyway. Removing rows with zero exposure is standard and correct; you are not discarding informative data, you are discarding rows that had no exposure to model.

Q7. The Poisson and quasi-Poisson gave the same residual deviance. Shouldn’t the fix change the fit? The deviance measures how far the fitted means sit from the data, and quasi-Poisson does not move the fitted means, so the deviance is unchanged. What changes is how we judge the deviance and the coefficients against chance: quasi-Poisson divides by the dispersion, so the same deviance now tells a more cautious story.

14.9 Practice problems¶

1. (A) List the three features of count data that the ordinary linear model fails to respect, and say which model assumption each one violates.

2. (A) State the mean and variance of a Poisson random variable. Why does this single fact make constant-variance least squares the wrong model for counts?

3. (A) In a Poisson log-linear model a predictor has coefficient $b = 0.5$. State its rate ratio and interpret it in words.

4. (A) Explain in one or two sentences what an offset is and why its coefficient is fixed at one rather than estimated.

5. (A) A Poisson fit and its quasi-Poisson refit report identical coefficients but different standard errors. Explain why.

6. (A) Define overdispersion and name two ordinary causes of it in real count data.

7. (A) Give the family and link function for each GLM row: linear regression, logistic regression, Poisson regression.

8. (A) A colleague computes crime rates as crimes / population and fits ordinary least squares. Give two reasons an offset-based Poisson model is better.

9. (A) Why can a Poisson model never predict a negative count, whatever the predictor values? Point to the piece of the model that guarantees it.

10. (B) Starting from the Poisson pmf $P(Y=y) = e^{-\mu}\mu^y/y!$, write the log-likelihood for independent counts $y_1, \dots, y_n$ with $\log\mu_i = \mathbf{x}_i'\boldsymbol\beta$, and derive the score equations $\sum_i (y_i - \mu_i)\mathbf{x}_i = \mathbf{0}$ (Theorem 14.3).

11. (B) From the score equations in problem 10, show that the intercept component forces $\sum_i \hat\mu_i = \sum_i y_i$: the fitted means reproduce the total observed count.

12. (B) Show that the Poisson score equations have the same form $\sum_i (y_i - \hat m_i)\mathbf{x}_i = \mathbf{0}$ as the logistic score equations of Chapter 13 and the normal equations of Chapter 2, identifying what $\hat m_i$ is in each case.

13. (B) Derive the offset result: starting from $\mu_i = \lambda_i t_i$ and $\log\lambda_i = \mathbf{x}_i'\boldsymbol\beta$, show that $\log\mu_i = \mathbf{x}_i'\boldsymbol\beta + \log t_i$, and explain why $\log t_i$ enters with coefficient one.

14. (B) Show that in a Poisson log-linear model with a single binary predictor $X \in \{0, 1\}$, the coefficient $\beta_1$ equals $\log(\mu_1/\mu_0)$, the log of the ratio of mean counts between the two groups. Hence $e^{\beta_1}$ is that ratio.

15. (B) The quasi-Poisson model assumes $\operatorname{Var}\{Y_i\} = \phi\mu_i$. Explain why this rescaling multiplies every coefficient standard error by $\sqrt\phi$ but leaves the coefficient estimates unchanged, referring to the score equations.

16. (B) Write the Pearson statistic $X^2 = \sum_i (y_i - \hat\mu_i)^2/\hat\mu_i$ and explain why $E\{X^2\} \approx n - p$ under a correct Poisson model, so that $\hat\phi = X^2/(n-p)$ should be near 1.

17. (B) For two Poisson means $\mu_A$ and $\mu_B$ with exposures $t_A, t_B$ entered as offsets, show that the fitted rate ratio $\lambda_B/\lambda_A = e^{\beta}$ does not depend on the exposures, so the comparison is a fair per-exposure comparison.

18. (B) Explain why maximum likelihood, not least squares, is used to fit the Poisson model, and describe in words how iteratively reweighted least squares reaches the estimates (you may cite the Chapter 13 sketch).

19. (C) Fit the Poisson model Species ~ Area + Elevation + Nearest + Scruz + Adjacent to gala.csv in R or Python. Report the coefficient on Area and its rate ratio for a 100 km$^2$ increase.

20. (C) For the Galapagos Poisson fit, compute the Pearson dispersion $\hat\phi$ from the Pearson residuals and confirm it matches the value 31.7 in the chapter. State what it says about the model.

21. (C) Refit the Galapagos model with family = quasipoisson (R) or the scale="X2" option (Python). Report which predictors are significant at the 5% level under quasi-Poisson but were significant under naive Poisson, and explain the change.

22. (C) Fit the ships damage model incidents ~ type + year + period with offset = log(service) on the positive-service rows. Report the rate ratio for period75 and interpret it.

23. (C) Refit the ships model without the offset (drop the offset argument). Compare the type coefficients to the offset model and explain why they change, connecting the change to the confounding of rate and exposure.

24. (C) Compute the ships Pearson dispersion $\hat\phi$, refit as quasi-Poisson, and report how the standard error of the typeB coefficient changes. Is the ship data as overdispersed as the Galapagos data?

25. (C) Fit the Toluca model hours ~ lotsize two ways: lm and glm with family = gaussian(link = "identity"). Show the coefficients agree, and explain what this says about the relationship between least squares and the GLM.

26. (C) Using the Galapagos Poisson fit, predict the expected species count for an island with Area = 40, Elevation = 200, Nearest = 1, Scruz = 30, Adjacent = 5. Do the prediction on the linear-predictor scale and then exponentiate, and confirm software gives the same number with type = "response".

27. (C) Fit a negative binomial model to the Galapagos data with MASS::glm.nb (R) or statsmodels NegativeBinomial (Python) and report the estimated $\theta$. Compare its Elevation coefficient to the Poisson and quasi-Poisson fits.

28. (C) Make a diagnostic plot of squared Pearson residuals against fitted means for the Galapagos Poisson fit, with the variance-equals-mean reference line, and describe what it shows about overdispersion.

29. (A) Chapter 10 tried log and Box-Cox transformations on the Galapagos species counts, and they helped the residual plots but never fully satisfied. Explain why the Poisson model succeeds where those transformations fell short. Name at least two specific defects of transforming a count that the count model avoids.

14.10 Exam practice¶

EP 14.1. Poisson regression makes two departures from ordinary least squares: it models the log of the mean rather than the mean itself, and it estimates the coefficients by maximum likelihood rather than by minimizing a sum of squares. Explain why each departure is forced on us by the nature of count data. Tie the log link to one specific property of counts, and tie the choice of maximum likelihood to another.

EP 14.2. A biologist fits a two-predictor Poisson model to the Galapagos data and hands you this output. Interpret the Elevation coefficient as a rate ratio for a rise of 100 meters, say what the exponentiated intercept means and whether it describes any real island, and state whether you would trust the tiny standard errors as they stand.

gala <- read.csv("data/gala.csv")
f <- glm(Species ~ Elevation + Nearest, family = poisson, data = gala)
summary(f)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.574e+00  3.819e-02  93.603  < 2e-16 ***
Elevation   1.475e-03  3.341e-05  44.151  < 2e-16 ***
Nearest     6.559e-03  1.444e-03   4.541 5.59e-06 ***

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 3510.7  on 29  degrees of freedom
Residual deviance: 1806.4  on 27  degrees of freedom

EP 14.3. In the ships damage model, fit with an offset of $\log(\text{service})$, a student sees the coefficient period75 $= 0.385$ and writes: “Ships operated in the later period suffer 0.385 more damage incidents than those in the earlier period.” Two separate things are wrong with that sentence. Identify and correct both, then write the statement the student should have made.

EP 14.4. The two fits below are the same ships model, once with the exposure offset and once without it. The type coefficients change sharply. Explain what happens to the type B coefficient when the offset is dropped and why, using the exposure figures, and state which of the two fits answers the question “which ship type is the most dangerous to operate.”

ships <- read.csv("data/ships.csv"); ships <- ships[ships$service > 0, ]
with_offset <- glm(incidents ~ type + year + period, family = poisson,
                   offset = log(service), data = ships)
no_offset   <- glm(incidents ~ type + year + period, family = poisson, data = ships)

type coefficients WITH offset:    typeB -0.5433  typeC -0.6874  typeD -0.0760  typeE  0.3256
type coefficients NO offset:      typeB  1.7957  typeC -1.2528  typeD -0.9045  typeE -0.1463

mean months of service by type:   A 1356   B 19760   C 885   D 635   E 855
total incidents by type:          A   42   B   253   C  12   D  17   E  32

EP 14.5. For the full five-predictor Galapagos Poisson fit the Pearson dispersion is $\hat\phi = 31.7$. (a) Explain in context what this single number says about the species counts and about the Poisson standard errors. (b) A colleague refits the same data with a negative binomial model and reports $\hat\theta = 1.675$ and an Elevation coefficient of 0.003855, against the Poisson value 0.003541. Interpret $\hat\theta$, and explain why the two elevation coefficients are nearly identical even though their standard errors are not.

Chapter game¶