Poisson regression

Topic of Standard errors, and models for counts, part of Generalised linear models · Also in Standard errors, and models for counts

Poisson regression

Type	Topic
Taught in	Standard errors, and models for counts · Standard errors, and models for counts
Part	Generalised linear models
Questions using it	2

A generalised linear model for counts with a Poisson random component and, usually, the log link: \log\mu_i=\mathbf{X}_i^\top\boldsymbol{\beta}, so e^{\beta_j} multiplies the expected count per unit of x_j. It assumes the variance equals the mean, which overdispersion breaks.

In Standard errors, and models for counts

From Standard errors, and models for counts

A count is a non-negative integer: insurance claims in a year, particles in a Geiger counter, episodes of a behaviour. The mean has to stay non-negative for every covariate. The natural link, and the canonical one, is the log: \log\mu_i={\boldsymbol{X}}_i^\top\boldsymbol{\beta} \qquad\text{equivalently}\qquad \mu_i=\exp({\boldsymbol{X}}_i^\top\boldsymbol{\beta}). The identity link \mu_i={\boldsymbol{X}}_i^\top\boldsymbol{\beta} is the second thing people reach for. It can look acceptable on a short stretch of x, where an exponential and a straight line are hard to tell apart. It is still not a model we will use. Somewhere outside the window the line is negative, and a mean count cannot be. The log link is a legitimate model for any {\boldsymbol{X}} and any \boldsymbol{\beta}.

The outer distribution for a count, as a starting point, is Poisson: Y_i\mid {\boldsymbol{X}}_i\stackrel{\mathrm{ind}}{\sim}\mathrm{Poisson}\!\left(\mu_i=\exp({\boldsymbol{X}}_i^\top\boldsymbol{\beta})\right). Poisson is the French for fish.¹ There is a precise limit behind the “rare events” slogan. If Y\sim\mathrm{Binomial}(n,p) and n\to\infty with p=\lambda/n, then Y converges in distribution to \mathrm{Poisson}(\lambda). Five runs in a week is, in that caricature, a coin tossed once a second for 604{,}800 seconds with success probability 5/604800, provided every second is independent of every other second and the probability never changes. Those are the Poisson-process assumptions, and they are the reason a real count is so often not Poisson. A Geiger counter has a dead time after each click. Incidents in a behavioural study arrive in spurts. Even a careful laboratory has a rate that wobbles. The variance is then not exactly the mean.

Quasi-Poisson

Main article: Quasi-Poisson

Unless there is a strong reason to insist on a pure Poisson process, use \mathrm{Var}(Y_i\mid {\boldsymbol{X}}_i)=\phi\mu_i=\phi\exp({\boldsymbol{X}}_i^\top\boldsymbol{\beta}). Here \phi>1 is overdispersion and \phi<1 is underdispersion. The case \phi=1 is the Poisson. This is a quasi-Poisson model. The variance function is a legitimate description of spread. Inside the exponential family, a variance exactly proportional to the mean forces the Poisson, so there is no exponential-family distribution with a free \phi in Equation (4) in Standard errors, and models for counts. You cannot simulate a new count from a quasi-Poisson specification alone. We will use it as a description. A stricter view, which reserves the name model for something one can simulate from, would refuse the name. We will come back to that distinction when we treat quasi-likelihood properly. For fitting, glm only needs the mean and the variance, so family = quasipoisson runs.

Reading a log-linear coefficient

Main article: Reading a log-linear coefficient

Two covariates are enough: \log\mu=\beta_0+\beta_1 x_1+\beta_2 x_2. Increase x_1 by one unit and hold x_2 fixed. The mean is multiplied by e^{\beta_1}. The same sentence with the labels swapped is \beta_2. The intercept is the mean when every covariate is at its baseline, \exp(\beta_0), and it is meaningful only when that baseline is a real situation. Do not judge a slope by whether the coefficient looks small. A small coefficient may be a large effect written in a small unit, and the comparison that matters is with the standard error.

Mine injuries

Main article: Mine-injury data

Myers, Montgomery, and Vining report injuries in the upper seam of 44 West Virginia coal mines. Three covariates survived their selection: inner burden thickness in feet (x_1), the percent already extracted from the seam below (x_2), and years the mine has been open (x_3). The fitted mean used for the interpretation is \hat\mu=\exp(-3.2707-0.0015 x_1+0.0627 x_2-0.0317 x_3). The glm summary from the fit below prints an intercept of -3.7207, with slopes -0.00148, 0.06270, and -0.03165. The slopes match the rounded values in Equation (5) in Standard errors, and models for counts. The intercepts do not: -3.2707 and -3.7207 differ by a transposition. The reading below follows Equation (5) in Standard errors, and models for counts as written, and the standard errors and the dispersion are taken from the printed summary, where the slopes are the quantities being tested.²

On Equation (5) in Standard errors, and models for counts, the intercept is the mean number of injuries when the inner burden is 0 feet, nothing has been extracted below, and the mine is brand new: \exp(-3.2707)=0.038. A new mine and a mine with nothing extracted below are imaginable. Zero inner burden is a hole at the surface, which is not a mine. The intercept is algebra, not a sentence to publish.

A one-foot increase in inner burden multiplies the mean by \exp(-0.0015)=0.999. A foot is about thirty centimetres. A per-foot effect is not worth quoting. Per ten feet the coefficient is -0.015, and the mean is multiplied by \exp(-0.015)\approx 0.985. Still a small multiplicative change, and still something to compare with a standard error rather than with 1 by eye.

One extra percentage point of extraction below multiplies injuries by \exp(0.0627)=1.065. One extra year of age multiplies them by \exp(-0.0317)=0.969, which is about 1.03 times fewer injuries. The sign is not the story you might have expected. A dangerous mine may have been closed before it could appear in a data set of mines that are still open. Age in this table is entangled with having survived.

The Poisson log-likelihood is \sum_i\{y_i\log\mu_i-\mu_i\}, up to a constant. Differentiating and setting the score to zero produces \sum_{i=1}^n\frac{y_i-\mu_i}{\mu_i}\,\mu_i\,{\boldsymbol{X}}_i =\sum_{i=1}^n(y_i-\mu_i)\,{\boldsymbol{X}}_i={\boldsymbol{0}} when the link is log, because then \mu'=\mu and \mathrm{Var}=\mu, so the two cancel. It is still nonlinear in \boldsymbol{\beta}, because \mu_i=\exp({\boldsymbol{X}}_i^\top\boldsymbol{\beta}). The general score from the previous chapter is the same equation before that cancellation. Software does the iteration.

need <- function(pkg) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    install.packages(pkg, repos = "https://cloud.r-project.org")
  }
  suppressPackageStartupMessages(library(pkg, character.only = TRUE))
}
need("locfit")
data(mine, package = "locfit")
mines <- data.frame(
  y = mine$frac,
  x1 = mine$inb,
  x2 = mine$extrp,
  x3 = mine$time
)
fit1 <- glm(y ~ x1 + x2 + x3, data = mines, family = poisson)
fit2 <- glm(y ~ x1 + x2 + x3, data = mines, family = quasipoisson)
summary(fit1)
summary(fit2)

The Poisson summary is

              Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.72068    0.97888   -3.801  0.000144
x1          -0.001479   0.000824  -1.794  0.0728
x2           0.06270    0.01227    5.110  3.2e-07
x3          -0.03165    0.01631   -1.941  0.0523

with dispersion fixed at 1. Inner burden and mine age are borderline. Percent extraction is not. The quasi-Poisson summary has the same four estimates. Only the standard errors change, because the estimated dispersion is 0.905, a bit under 1: given these three covariates, the injury counts vary less than a Poisson process. The quasi-Poisson standard errors are smaller by about \sqrt{0.905}, the t statistics move accordingly, and mine age now sits below 0.05 (p=0.048) while inner burden remains borderline (p=0.067). All three terms stay in the model. A term that is hovering on the boundary, and that still carries information, is not dropped. The underdispersion fits the selection story: mines that would have been wild never stayed open long enough to be counted, so the ones we see are a little more deterministic than a pure count process.

In Standard errors, and models for counts, continued