In Poisson regression, the response variable \(Y\) is an occurrence count recorded for a particular measurement window. Usually, this window is a length of time, but it can also be a distance, area, etc. For example, \(Y\) could count the number of flaws in a manufactured tabletop of a certain area. If the observations recorded correspond to different measurement windows, a scale adjustment has to be made to put them on equal terms, and we model the rate or count per measurement unit \(t\).

For the random component, we assume that the response \(Y\) has a Poisson distribution. That is, \(Y_i\sim Poisson(\mu_i)\), for \(i=1, \ldots, N\) where the expected count of \(Y_i\) is \(E(Y_i)=\mu_i\). The link function is usually the (natural) log, but sometimes the identity function may be used. The systematic component consists of a linear combination of explanatory variables \((\alpha+\beta_1x_1+\cdots+ \beta_kx_k\)); this is identical to that for logistic regression. Thus, in the case of a single explanatory, the model is written

\(\log(\mu)=\alpha+\beta x\)

This is equivalent to

\(\mu=\exp(\alpha+\beta x)=\exp(\alpha)\exp(\beta x)\).

Interpretations of these parameters are similar to those for logistic regression. \(\exp(\alpha)\) is the effect on the mean of \(Y\) when \(x = 0\), and \(\exp(\beta)\) is the multiplicative effect on the mean of \(Y\) for each 1-unit increase in \(x\).

• If \(\beta = 0\), then \(\exp(\beta) = 1\), and the expected count, \( \mu = E(Y) = \exp(\beta)\), and \(Y\) and \(x\) are not related.
• If \(\beta > 0\), then \(\exp(\beta) > 1\), and the expected count \( \mu = E(Y)\) is \(\exp(\beta)\) times larger than when \(x = 0\).
• If \(\beta< 0\), then \(\exp(\beta) < 1\), and the expected count \( \mu = E(Y)\) is \(\exp(\beta)\) times smaller than when \(x = 0\).

Compared with the model for count data above, we can alternatively model the expected rate of observations per unit of length, time, etc. to adjust for data collected over differently-sized measurement windows. For example, if \(Y\) is the count of flaws over a length of \(t\) units, then the expected value of the rate of flaws per unit is \(E(Y/t)=\mu/t\). For a single explanatory variable, the model would be written as

\(\log(\mu/t)=\log\mu-\log t=\alpha+\beta x\)

The term \(\log t\) is referred to as an offset. It is an adjustment term and a group of observations may have the same offset, or each individual may have a different value of \(t\). The term \(\log(t)\) is an observation, and it will change the value of the estimated counts:

Similar to the case of logistic regression, the maximum likelihood estimators (MLEs) for \(\beta_0, \beta_1 \dots \), etc.) are obtained by finding the values that maximize the log-likelihood. In general, there are no closed-form solutions, so the ML estimates are obtained by using iterative algorithms such as Newton-Raphson (NR), Iteratively re-weighted least squares (IRWLS), etc.

The usual tools from the basic statistical inference of GLMs are valid:

• Confidence Intervals and Hypothesis tests for parameters
  • Wald statistics and asymptotic standard error (ASE)
  • Likelihood ratio tests
• Distribution of probability estimates

In the next, we will take a look at an example using the Poisson regression model for count data with SAS and R. In SAS we can use PROC GENMOD which is a general procedure for fitting any GLM. Many parts of the input and output will be similar to what we saw with PROC LOGISTIC. In R we can still use glm(). The response counts are recorded for the same measurement windows (horseshoe crabs), so no scale adjustment for modeling rates is necessary.

In the above model, we detect a potential problem with overdispersion since the scale factor, e.g., Value/DF, is greater than 1.

What does overdispersion mean for Poisson Regression? What does it tell us about the relationship between the mean and the variance of the Poisson distribution for the number of satellites? Recall that one of the reasons for overdispersion is heterogeneity, where subjects within each predictor combination differ greatly (i.e., even crabs with similar width have a different number of satellites). If that's the case, which assumption of the Poisson model is violated?

proc genmod data=crab;
model Sa=w / dist=poi link=log scale=pearson;
output out=data2 p=pred;
run;
proc print;
run;

model.disp = glm(Sa ~ W, family=quasipoisson(link=log), data=hcrabs)
summary.glm(model.disp)
summary.glm(model.disp)$dispersion

> summary.glm(model.disp)

Call:
glm(formula = Sa ~ W, family = quasipoisson(link = log), data = hcrabs)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.8526  -1.9884  -0.4933   1.0970   4.9221  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.30476    0.96729  -3.417 0.000793 ***
W            0.16405    0.03562   4.606 7.99e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 3.182205)

    Null deviance: 632.79  on 172  degrees of freedom
Residual deviance: 567.88  on 171  degrees of freedom

With this model, the random component does not technically have a Poisson distribution any more (hence the term "quasi" Poisson) because that would require that the response has the same mean and variance. From the estimate given (Pearson \(X^2/171 = 3.1822\)), the variance of the number of satellites is roughly three times the size of the mean.

The new standard errors (in comparison to the model without the overdispersion parameter), are larger, (e.g., \(0.0356 = 1.7839( 0.02)\) which comes from the scaled SE (\(\sqrt{3.1822}=1.7839\)); the adjusted standard errors are multiplied by the square root of the estimated scale parameter. Thus, the Wald statistics will be smaller and less significant.

 What could be another reason for poor fit besides overdispersion? Can we improve the fit by adding other variables?

As it turns out, the color variable was actually recorded as ordinal with values 2 through 5 representing increasing darkness and may be quantified as such. So, we next consider treating color as a quantitative variable, which has the advantage of allowing a single slope parameter (instead of multiple indicator slopes) to represent the relationship with the number of satellites.

proc genmod data=crab;
model Sa=w c / dist=poi link=log  scale=pearson;
run;

data(hcrabs)
colnames(hcrabs) = c("C","S","W","Sa","Wt")
hcrabs$C = as.numeric(hcrabs$C)+1
model = glm(Sa ~ W+C, family=quasipoisson(link=log), data=hcrabs)
summary(model)
anova(model, test='Chisq')

> summary(model)

Call:
glm(formula = Sa ~ W + C, family = quasipoisson(link = log), 
    data = hcrabs)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.9919  -1.9539  -0.5526   0.9836   4.7596  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -2.35058    1.15156  -2.041   0.0428 *  
W            0.14957    0.03716   4.025 8.55e-05 ***
C           -0.16940    0.11112  -1.524   0.1292    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 3.228764)

    Null deviance: 632.79  on 172  degrees of freedom
Residual deviance: 560.20  on 170  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 6

> anova(model, test='Chisq')
Analysis of Deviance Table

Model: quasipoisson, link: log

Response: Sa

Terms added sequentially (first to last)


     Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                   172     632.79              
W     1   64.913       171     567.88 7.332e-06 ***
C     1    7.677       170     560.20    0.1231    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The estimated model is now

\(\log{\hat{\mu_i}}= -2.3506 + 0.1496W_i - 0.1694C_i\)

And the interpretation of the single slope parameter for color is as follows: for each 1-unit increase in the color (darkness level), the expected number of satellites is multiplied by \(\exp(-.1694)=.8442\).

In terms of the fit, adding the numerical color predictor doesn't seem to help; the overdispersion seems to be due to heterogeneity. Is there something else we can do with this data? We can either (1) consider additional variables (if available), (2) collapse over levels of explanatory variables, or (3) transform the variables.

Let's consider grouping the data by the widths and then fitting a Poisson regression model that models the rate of satellites per crab.

In this approach, we create 8 width groups and use the average width for the crabs in that group as the single representative value. While width is still treated as quantitative, this approach simplifies the model and allows all crabs with widths in a given group to be combined. The disadvantage is that differences in widths within a group are ignored, which provides less information overall. However, another advantage of using the grouped widths is that the saturated model would have 8 parameters, and the goodness of fit tests, based on \(8-2\) degrees of freedom, are more reliable.

To account for the fact that width groups will include different numbers of crabs, we will model the mean rate \(\mu/t\) of satellites per crab, where \(t\) is the number of crabs for a particular width group. This variable is treated much like another predictor in the data set. The difference is that this value is part of the response being modeled and not assigned a slope parameter of its own. We will see more details on the Poisson rate regression model in the next section.

The data, after being grouped into 8 intervals, is shown in the table below. Two columns to note in particular are "Cases", the number of crabs with carapace widths in that interval, and "Width", which now represents the average width for the crabs in that interval.

    Interval     Cases   Width   AvgSa    SDSa    VarSa
      =23.25        14  22.693  1.0000  1.6641   2.7692
 23.25-24.25        14  23.843  1.4286  2.9798   8.8792
 24.25-25.25        28  24.775  2.3929  2.5581   6.5439
 25.25-26.25        39  25.838  2.6923  3.3729  11.3765
 26.26-27.25        22  26.791  2.8636  2.6240   6.8854
 27.25-28.25        24  27.738  3.8750  2.9681   8.8096
 28.25-29.25        18  28.667  3.9444  4.1084  16.8790
      >29.25        14  30.407  5.1429  2.8785   8.2858

In this approach, each observation within a group is treated as if it has the same width.

data crabgrp;
input width cases SaTotal AvgSa;
lcases=log(cases);
cards;
22.69 14 14 1.00
23.84 14 20 1.428571
24.77 28 67 2.392857
25.84 39 105 2.692308
26.79 22 63 2.863636
27.74 24 93 3.875
28.67 18 71 3.944444
30.41 14 72 5.142857
;
proc genmod data=crabgrp order=data;
model SaTotal = width / dist=poisson link=log 
 offset=lcases obstats residuals;
 output out=data3 p=pred;
run;

# creating the 8 width groups
width.long = cut(hcrabs$W, breaks=c(0,seq(23.25, 29.25),Inf))

# summarizing per group
Cases = aggregate(hcrabs$W, by=list(W=width.long),length)$x
lcases = log(Cases)
Width = aggregate(hcrabs$W, by=list(W=width.long),mean)$x
AvgSa = aggregate(hcrabs$Sa, by=list(W=width.long),mean)$x
SdSa = aggregate(hcrabs$Sa, by=list(W=width.long),sd)$x
VarSa = aggregate(hcrabs$Sa, by=list(W=width.long),var)$x
SaTotal = aggregate(hcrabs$Sa, by=list(W=width.long),sum)$x
cbind(table(width.long), Width, AvgSa, SdSa, VarSa)

hcrabs.grp = data.frame(Width=Width, lcases=lcases, SaTotal=SaTotal)
model.grp = glm(SaTotal ~ Width, offset=lcases, family=poisson, data=hcrabs.grp)
summary(model.grp)

 summary(model.grp)

Call:
glm(formula = SaTotal ~ Width, family = poisson, data = hcrabs.grp, 
    offset = lcases)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5322  -0.7750  -0.3454   0.7151   1.0347  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.54018    0.57658  -6.140 8.26e-10 ***
Width        0.17290    0.02125   8.135 4.11e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 72.3772  on 7  degrees of freedom
Residual deviance:  6.5164  on 6  degrees of freedom

> rstandard(model.grp)
         1          2          3          4          5          6          7          8 
-1.7312682 -1.1474695  1.1980218  0.7447990 -0.3414285  1.0400246 -0.4260710 -1.0213629 
> predict.glm(model.grp, type="response", newdata=hcrabs.grp)
       1        2        3        4        5        6        7        8 
20.54097 25.05952 58.88382 98.57261 65.55897 84.23620 74.18733 77.96057 

In this lesson, we showed how the generalized linear model can be applied to count data, using the Poisson distribution with the log link. Compared with the logistic regression model, two differences we noted are the option to use the negative binomial distribution as an alternate random component when correcting for overdispersion and the use of an offset to adjust for observations collected over different windows of time, space, etc. Much of the properties otherwise are the same (parameter estimation, deviance tests for model comparisons, etc.).

One other common characteristic between logistic and Poisson regression that we change for the log-linear model coming up is the distinction between explanatory and response variables. The log-linear model makes no such distinction and instead treats all variables of interest together jointly. This allows greater flexibility in what types of associations can be fit and estimated, but one restriction in this model is that it applies only to categorical variables.