R Data Analysis Examples: Poisson Regression



Stat Computing > R > Data Analysis Examples

R Data Analysis Examples
Poisson Regression

Note: the code on this page works with R 2.4.1.

Examples of Poisson Regression

Example 1. School administrators study the attendance behavior of high school juniors at two schools. Predictors of the number of days of absence include gender of the student and standardized test scores in math and language arts.

Description of the Data

Let's pursue Example 1 from above. We have attendance data on 316 high school juniors from two urban high schools in the file poissonreg.csv The response variable of interest is days absent, daysabs. The variables math and langarts give the standardized test scores for math and language arts respectively. The variable male is a binary indicator of student gender.

Let's look at the data.

p<-read.table("http://www.ats.ucla.edu/stat/R/dae/poissonreg.csv", sep=",", header = TRUE) 
attach(p)
names(p)
[1] "id"       "school"   "male"     "math"     "langarts" "daysatt"  "daysabs"

The response variable of interest is days absent, daysabs. The variables math and langarts give the standardized test scores for math and language arts respectively. The variable male is a binary indicator of student gender.

Let's look at the data. For a nice layout for displaying the summary statistics, we made use of the package called fields.

library(fields)
c<-cbind(daysabs, math, langarts, male)
stats(c)
                  daysabs       math   langarts        male
N              316.000000 316.000000 316.000000 316.0000000
mean             5.810127  48.751012  50.063794   0.4873418
Std.Dev.         7.449003  17.880757  17.939211   0.5006325
min              0.000000   1.007114   1.007114   0.0000000
Q1               1.000000  37.725357  40.150265   0.0000000
median           3.000000  48.943764  50.000000   0.0000000
Q3               8.000000  61.043877  61.043877   1.0000000
max             45.000000  98.992889  98.992889   1.0000000
missing values   0.000000   0.000000   0.000000   0.0000000
var(daysabs)
[1] 55.48764

hist(daysabs, breaks=150, xlim=c(0, 50), col="gray")

table(male)

male
  0   1 
162 154

Some Strategies You Might Be Tempted To Try

Before we show how you can analyze this with a Poisson regression analysis, let's consider some other methods that you might use.

OLS Regression - You could try to analyze these data using OLS regression. However, count data are highly non-normal and are not well estimated by OLS regression.
Negative Binomial Regression - Negative binomial regression does better with over dispersed data, i.e. variance much larger than the mean.
Zero-inflated Regression Model - This might be a good if there are more zeros than would be expected by either a Poisson or negative binomial model.

R Poisson Regression Analysis

R function glm (generalized linear model) is used here for fitting the Poisson model.

m1<-glm(daysabs~math+langarts+male, family=poisson)
summary(m1)

Call:
glm(formula = daysabs ~ math + langarts + male, family = poisson)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-4.0117  -2.5455  -1.1014   0.8956  11.0374  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.687666   0.072651  36.994  < 2e-16 ***
math        -0.003523   0.001821  -1.934   0.0531 .  
langarts    -0.012152   0.001835  -6.623 3.52e-11 ***
male        -0.400921   0.048412  -8.281  < 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: 2409.8  on 315  degrees of freedom
Residual deviance: 2234.5  on 312  degrees of freedom
AIC: 3103.9

Number of Fisher Scoring iterations: 6

The output looks very much like the output from an OLS regression. The output begins with the information on deviance residuals. It gives us some idea on how well the model fits the data. You will then see the Poisson regression coefficients for each of the variables along with the standard errors, z-scores, and p-values. By default, we assume that there is no dispersion, meaning the mean and the variance are the same and this is indicated by the message stating that the dispersion parameter is fixed at 1. The last part of the output gives the model fit information, including AIC. The residual deviance and the null deviance are calculated based on the formula below.

The change in deviance between the null model and the full model measures how much better the full model fits the data than the null model (which is the intercept-only model.) We can also perform a Wald-test to get the overall model fit. This requires the lmtest package.

library(lmtest)
waldtest(m1, test="Chisq")

Wald test

Model 1: daysabs ~ math + langarts + male
Model 2: daysabs ~ 1
  Res.Df  Df  Chisq Pr(>Chisq)    
1    312                          
2    315  -3 176.24  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

This tells us that overall our model is doing significantly better than the null model.

Now, just to be on the safe side, let's rerun the poisson model with the robust standard errors. This requires one more package called sandwich, beside lmtest.

library(sandwich)
coeftest(m1, vcov=sandwich)

z test of coefficients:

              Estimate Std. Error z value  Pr(>|z|)    
(Intercept)  2.6876659  0.2177980 12.3402 < 2.2e-16 ***
math        -0.0035232  0.0076252 -0.4621  0.644044    
langarts    -0.0121521  0.0052867 -2.2986  0.021527 *  
male        -0.4009209  0.1393705 -2.8767  0.004019 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The robust standard errors attempt to adjust for heterogeneity in the model. Using the robust standard errors has resulted in a fairly large change in the standard error, which should be more appropriate. The z-tests still yield similar significant results, but give more realistic p-values.

Since math is not significant in the model with robust standard errors, we will rerun the model dropping that variable.

m2 <- glm(daysabs ~ langarts + male, family = poisson)
summary(m2)

Call:
glm(formula = daysabs ~ langarts + male, family = poisson)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-4.297  -2.510  -1.123   0.869  10.495  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.646976   0.069776  37.935   <2e-16 ***
langarts    -0.014670   0.001293 -11.342   <2e-16 ***
male        -0.409353   0.048219  -8.489   <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: 2409.8  on 315  degrees of freedom
Residual deviance: 2238.3  on 313  degrees of freedom
AIC: 3105.7

Number of Fisher Scoring iterations: 6

coeftest(m2, vcov=sandwich)

z test of coefficients:

              Estimate Std. Error z value  Pr(>|z|)    
(Intercept)  2.6469765  0.1823472 14.5161 < 2.2e-16 ***
langarts    -0.0146700  0.0034328 -4.2735 1.924e-05 ***
male        -0.4093528  0.1352244 -3.0272  0.002468 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

waldtest(m2, test="Chisq", vcov=vcovHC)

Wald test

Model 1: daysabs ~ langarts + male
Model 2: daysabs ~ 1
  Res.Df  Df  Chisq Pr(>Chisq)    
1    313                          
2    315  -2 26.462  1.794e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Finally, we will use the predict function to get the predicted values in days absent. For example, we fix the value of langarts at its mean and male at 0 and 1, the predict function returns the expected count for both male and female when the language arts score is held at its mean.

male<-c(0, 1)
ml<-mean(langarts)
langarts<-c(ml, ml)
pnew<-cbind(male, langarts)
fitted<-predict(m2, newdata=data.frame(pnew), type="response")
fitted

       1        2 
6.770267 4.495998

Sample Write-Up of the Analysis

Before we begin the sample write-up we need to get the output into a form more acceptable for publication, such as in Latex format. The Latex code can be generated easily using the function xtable from a package called xtable.

library(xtable)
x<-xtable(m2)
x
% latex table generated in R 2.4.1 by xtable 1.4-5 package
% Mon May 28 10:03:40 2007
\begin{table}[ht]
\begin{center}
\begin{tabular}{rrrrr}
  \hline
 & Estimate & Std. Error & z value & Pr($>$$|$z$|$) \\
  \hline
(Intercept) & 2.6470 & 0.0698 & 37.94 & 0.0000 \\
  male & $-$0.4094 & 0.0482 & $-$8.49 & 0.0000 \\
  langarts & $-$0.0147 & 0.0013 & $-$11.34 & 0.0000 \\
   \hline
\end{tabular}
\end{center}
\end{table}

The LaTex code above generates a table as shown below:

The Poisson regression model predicting days absent from school stay from language arts and gender was statistically significant (chi-squared = 26.462, df = 2, p<.0001). The predictors langarts and male were each statically significant. For these data, the expected log count for a one-unit increase in language arts was -0.0146. This translates to a decrease of about 1.46 days absent for a one standard deviation increase in language arts when gender is held constant. Male students had an expected log count -0.41 less than female students which amounts to about 2.27 fewer days absent than females while holding language arts constant.

Cautions, Flies in the Ointment

It is not recommended that Poisson models be applied to small samples. What constitutes a small sample does not seem to be clearly defined in the literature.

R Data Analysis Examples Poisson Regression