Model

In our dataset we have the ability to explore the relationship between death and smoking. We should make sure to think about this problem in the poisson frame with respect to \(\lambda\) (a rate), specifically the rate of death per individual. Note that our data just gives counts of death per group population.

So we have 2 things:

• the counts of deaths per group
• and \(\lambda\) the rate of death per person

They obviously relate with a rate = count of deaths / total individuals.

If we are looking at the regression

\[\begin{align*} \log(\lambda | X) &= X\beta \\ \log(count/total |X) &= X\beta \\ \Leftrightarrow \log(count |X) &= X\beta + \log(total) \end{align*}\]

In the final transformation we call \(log(total)\) the offset. There are cases when we don’t use offsets (e.g. number of papers a graduate student has published), but offsets can be like our example or more complicated, like area of region for growing mushrooms (vs number of mushrooms) etc. Also notice that the \(\log(total)\) does not have a \(\beta\) value associated with it.

Additionally, in the graduate student example (where the population is 1 individual), you can also think of the offset being \(\log(1) = 0\).

As such, in R we’d do the following:

poisson_glm <- glm(dead ~ age_int + smoke, offset = log(pop), family = poisson,
                   data = smoking)

The summary looks decently similar to our other glms.

summary(poisson_glm)

## 
## Call:
## glm(formula = dead ~ age_int + smoke, family = poisson, data = smoking, 
##     offset = log(pop))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.7322  -1.3737  -0.1642   0.7407   2.7613  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -6.249136   0.088890 -70.302  < 2e-16 ***
## age_int              0.068141   0.001154  59.073  < 2e-16 ***
## smokecigarretteOnly  0.390484   0.037330  10.460  < 2e-16 ***
## smokecigarrettePlus  0.189927   0.035892   5.292 1.21e-07 ***
## smokeno             -0.039428   0.046880  -0.841      0.4    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4055.984  on 35  degrees of freedom
## Residual deviance:   67.681  on 31  degrees of freedom
## AIC: 317.7
## 
## Number of Fisher Scoring iterations: 4

Hmmmmmm… if we look at those levels of smoking it isn’t really as we expected… we really want the first level to “no” so that all those \(\beta\) values compare against not smoking…

smoking <- smoking %>% mutate(smoke = fct_relevel(factor(smoke), "no"))

poisson_glm2 <- glm(dead ~ age_int + smoke, offset = log(pop), family = poisson,
                   data = smoking) 
poisson_glm2 %>% summary

## 
## Call:
## glm(formula = dead ~ age_int + smoke, family = poisson, data = smoking, 
##     offset = log(pop))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.7322  -1.3737  -0.1642   0.7407   2.7613  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -6.288564   0.086484 -72.713  < 2e-16 ***
## age_int              0.068141   0.001154  59.073  < 2e-16 ***
## smokecigarPipeOnly   0.039428   0.046880   0.841      0.4    
## smokecigarretteOnly  0.429912   0.039753  10.815  < 2e-16 ***
## smokecigarrettePlus  0.229356   0.038565   5.947 2.73e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4055.984  on 35  degrees of freedom
## Residual deviance:   67.681  on 31  degrees of freedom
## AIC: 317.7
## 
## Number of Fisher Scoring iterations: 4

Dispersion correction

Before we talk about inference of the \(\beta\) values and how to interpret them, let’s look at dispersion - which we talked about needing to do the get good confidence intervals.

Recall that since the variance and mean of the poisson are directly linked dispersion tells us if that “linking” is actually correct, we can caculate an estimate of \(\phi\) in the equation:

\[ Var(X) = \phi \lambda \]

when \(X\) is poisson.

We can test whether the dispersion is different than 1 with a \(\chi^2\) test:

dispersion <-  poisson_glm$deviance/poisson_glm$df.residual
pchisq(poisson_glm$deviance,df = poisson_glm$df.residual,
       lower.tail = FALSE)

## [1] 0.0001527248

Looks like \(\phi\) isn’t 1, that’s ok we can use a quasipoisson model.

qpoisson_glm  <- glm(dead ~ age_int + smoke, offset = log(pop),
                                   family = quasipoisson, data = smoking)

summary(qpoisson_glm)

## 
## Call:
## glm(formula = dead ~ age_int + smoke, family = quasipoisson, 
##     data = smoking, offset = log(pop))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.7322  -1.3737  -0.1642   0.7407   2.7613  
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -6.288564   0.125073 -50.279  < 2e-16 ***
## age_int              0.068141   0.001668  40.847  < 2e-16 ***
## smokecigarPipeOnly   0.039428   0.067797   0.582 0.565063    
## smokecigarretteOnly  0.429912   0.057490   7.478    2e-08 ***
## smokecigarrettePlus  0.229356   0.055772   4.112 0.000267 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 2.09146)
## 
##     Null deviance: 4055.984  on 35  degrees of freedom
## Residual deviance:   67.681  on 31  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4

Question: How does the raw \(\beta\) estimates change between the Quasi-Poisson model (that estimates \(\phi\)) and the standard Poisson model?