8.7 Revisisting Overdispersion
Let’s run our model fit test based on the deviance for model2.
1 - pchisq(summary(model2)$deviance, summary(model2)$df.residual )## [1] 0Again, we reject the hypothesis of a close fit between model and data. To gain a little more insight we can plot estimates of the variance against the expected value, alongside a line with an intercept of zero and a slope of 1.
We expect the data points to fall somewhat evenly along that line. Here, it appears our variance is consistently larger than our mean, indicating the possibility of overdispersion.
plot(
log(fitted(model2)),
log((ferraro2016$morbidityw1-fitted(model2))^2),
xlab=expression(hat(mu)),
ylab=expression((y-hat(mu))^2),
pch=20,col="blue"
)
abline(0,1) ## 'varianc = mean' line
Two ways we might handle this overdispersion are (1) estimate the overdispersion parameter directly within the model, or (2) use a negative binomial model.
8.7.1 Quassi-Poisson Family
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.
model3 <- glm(
formula = morbidityw1 ~ 1 + abuse_rare + income_star + abuse_rare:income_star,
family = quasipoisson(link=log),
data = ferraro2016,
na.action = na.exclude
)
summary(model3)##
## Call:
## glm(formula = morbidityw1 ~ 1 + abuse_rare + income_star + abuse_rare:income_star,
## family = quasipoisson(link = log), data = ferraro2016, na.action = na.exclude)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4949 -1.3511 -0.4745 0.6512 8.6487
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.0766511 0.0217425 49.518 < 2e-16 ***
## abuse_rare1 -0.2382614 0.0421387 -5.654 1.71e-08 ***
## income_star -0.0311089 0.0153759 -2.023 0.0431 *
## abuse_rare1:income_star -0.0000822 0.0300781 -0.003 0.9978
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 2.812535)
##
## Null deviance: 8068.4 on 2966 degrees of freedom
## Residual deviance: 7956.7 on 2963 degrees of freedom
## (55 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5The new standard errors (in comparison to the model without the overdispersion parameter), are larger, Thus, the Wald statistics will be smaller and less likely to be significant.
8.7.2 Negative Binomial Regression
We can also fit a negative binomial model, in which allows dispersion or variance in the outcome. The negative binomial distribution has one parameter more than the Poisson regression. This parameters adjusts the variance independently from the mean. In fact, the Poisson distribution is a special case of the negative binomial distribution.
Importantly, the Poisson and Negative Binomial have the same mean structure, so we can interpret coefficients in the same way. However, we will need to use the MASS package (Venables and Ripley 2002) to fit the Negative Binomial regression.
model4 <- MASS::glm.nb(
formula = morbidityw1 ~ 1 + abuse_rare + income_star + abuse_rare:income_star,
data = ferraro2016,
na.action = na.exclude
)
summary(model4)##
## Call:
## MASS::glm.nb(formula = morbidityw1 ~ 1 + abuse_rare + income_star +
## abuse_rare:income_star, data = ferraro2016, na.action = na.exclude,
## init.theta = 1.529056769, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8428 -0.8704 -0.2930 0.3709 3.7762
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0766474 0.0221378 48.634 < 2e-16 ***
## abuse_rare1 -0.2382366 0.0406825 -5.856 4.74e-09 ***
## income_star -0.0311705 0.0156452 -1.992 0.0463 *
## abuse_rare1:income_star -0.0004478 0.0289957 -0.015 0.9877
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.5291) family taken to be 1)
##
## Null deviance: 3375.4 on 2966 degrees of freedom
## Residual deviance: 3334.7 on 2963 degrees of freedom
## (55 observations deleted due to missingness)
## AIC: 12770
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.5291
## Std. Err.: 0.0682
##
## 2 x log-likelihood: -12760.4610We notice the smaller dispersion parameter, approximately \(1.53\), providing some support the negative binomial is a better fit for our data.
We can also test the hypothesis of overdispersion formally using a likelihood ratio test. The difference between the two models is captured by estimating a dispersion parameter that is held constant in a Poisson model. Thus, the Poisson model is actually nested in the negative binomial model. We can then use a likelihood ratio test to compare the two and test this model assumption.
pchisq(2 * (logLik(model4) - logLik(model2)), df = 1, lower.tail = FALSE)## 'log Lik.' 0 (df=5)In this example the associated chi-squared value estimated from \(2*(logLik(\mathrm{negative \:binomial}) – logLik(\mathrm{poisson}))\) is \(1968.679\) with one degree of freedom. This strongly suggests the negative binomial model, estimating the dispersion parameter, is more appropriate than the Poisson model.