As we saw in logistic regression, if we want to test and adjust for overdispersion we can add a scale parameter with the family=quasipoisson option. The estimated scale parameter will be labeled as "Overdispersion parameter" in the output. (As stated earlier we can also fit a negative binomial regression instead)

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

Below is the output when using the quasi-Poisson model. How does this compare to the output above from the earlier stage of the code? Do we have a better fit now?

> 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

To add the horseshoe crab color as a categorical predictor (in addition to width), we can use the following code. The change of baseline to the 5th color is arbitrary.

hcrabs$C = factor(hcrabs$C, levels=c('5','2','3','4'))
model = glm(Sa ~ W+C, family=quasipoisson(link=log), data=hcrabs)
summary(model)
anova(model, test='Chisq')

Note in the output that there are three separate parameters estimated for color, corresponding to the three indicators included for colors 2, 3, and 4 (5 as the baseline). We continue to adjust for overdispersion withfamily=quasipoisson, although we could relax this if adding additional predictor(s) produced an insignificant lack of fit.

> summary(model)

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

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.0415  -1.9581  -0.5575   0.9830   4.7523  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.09740    1.00260  -3.089  0.00235 ** 
W            0.14934    0.03748   3.985  0.00010 ***
C2           0.44736    0.37604   1.190  0.23586    
C3           0.24767    0.29339   0.844  0.39978    
C4           0.01100    0.32442   0.034  0.97300    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 3.233628)

    Null deviance: 632.79  on 172  degrees of freedom
Residual deviance: 559.34  on 168  degrees of freedom

The estimated model is: \(\log{\hat{\mu_i}}= -3.0974 + 0.1493W_i + 0.4474C_{2i}+ 0.2477C_{3i}+ 0.0110C_{4i}\), using indicator variables for the first three colors.

There does not seem to be a difference in the number of satellites between any color class and the reference level 5 according to the chi-squared statistics for each row in the table above. Furthermore, by the ANOVA output below we see that color overall is not statistically significant after we consider the width.

> 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.449e-06 ***
C     3    8.534       168     559.34    0.4507    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Is this model preferred to the one without color?

To add color as a quantitative predictor, we first define it as a numeric variable. Although the original values were 2, 3, 4, and 5, R will by default use 1 through 4 when converting from factor levels to numeric values. So, we add 1 after the conversion.

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')

Compare these partial parts of the output with the output above where we used color as a categorical predictor. We are doing this to keep in mind that different coding of the same variable will give us different fits and estimates.

> 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

In R, the lcases variable is specified with the OFFSET option, which takes the log of the number of cases within each grouping. This is our adjustment value \(t\) in the model that represents (abstractly) the measurement window, which in this case is the group of crabs with similar width. So, \(t\) is effectively the number of crabs in the group, and we are fitting a model for the rate of satellites per crab, given carapace width.

# 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)

Here is the output. Note the "offset = lcases" under the model expression.

 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

Based on the Pearson and deviance goodness of fit statistics, this model clearly fits better than the earlier ones before grouping width. For example, the Value/DF for the deviance statistic now is 1.0861. The estimated model is:

\(\log (\hat{\mu}_i/t)= -3.54 + 0.1729\mbox{width}_i\)

For each 1-cm increase in carapace width, the mean number of satellites per crab is multiplied by \(\exp(0.1729)=1.1887\).

> 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 

The residuals analysis indicates a good fit as well, and the predicted values correspond a bit better to the observed counts in the "SaTotal" cells. Let's compare the observed and fitted values in the plot below: