I have some overdispersed count data (mean=8.6, var=263.5) that I am hoping to model using either a negative binomial (NB) or quasipoisson GLM. A rootogram suggests that the Poisson distribution is a poor fit.

I compared three different models: Poisson, NB, and quasipoisson.

fake = structure(list(factor = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                       1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                                       2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
                                       3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
                                       4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("a", "b", "c", 
                                                                                       "d"), class = "factor"), freq = c(0L, 0L, 2L, 2L, 0L, 16L, 0L, 
                                                                                                                         29L, 0L, 0L, 0L, 0L, 0L, 4L, 0L, 4L, 0L, 0L, 1L, 8L, 0L, 11L, 
                                                                                                                         9L, 37L, 0L, 0L, 1L, 0L, 5L, 10L, 0L, 1L, 2L, 4L, 14L, 33L, 0L, 
                                                                                                                         64L, 19L, 56L, 2L, 2L, 0L, 1L, 0L, 12L, 10L, 15L, 7L, 2L, 19L, 
                                                                                                                         7L, 1L, 78L, 1L, 49L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 10L)), .Names = c("factor", 
                                                                                                                                                                                             "freq"), row.names = c(NA, -64L), class = "data.frame")
rs1 = glm(freq ~ factor, data=fake, family="poisson")

rs2 = glm.nb(freq ~ factor, data=fake)

rs3 = glm(freq ~ factor, data=fake, family="quasipoisson")

#compare Poisson with NB
lrtest(rs1, rs2)

The model coefficients turn out exactly the same for all 3 models, but with different coefficient SEs. Nevertheless, the LR test suggests that the NB fits much better than the Poisson.

  1. What is going on here?
  2. When I figure out this behavior, would it be reasonable to compare response residual magnitude and distribution to choose between quasipoisson and NB, since the former doesn't have likelihood-based statistics?
  • $\begingroup$ What makes you think the results are weird, exactly? $\endgroup$ – Glen_b Jun 25 '13 at 3:14
  • $\begingroup$ I don't understand why the three models should have exactly the same coefficients. I know NB reduces to Poisson if there is no overdispersion (which is definitely not the case here), but quasipoisson should be different, right? Also why is the LR test so strongly favoring the NB? $\endgroup$ – half-pass Jun 25 '13 at 4:54
  • 2
    $\begingroup$ No, quasipoisson should not be different; the only difference between quasi-Poisson and Poisson is the estimate of dispersion. $\endgroup$ – Glen_b Jun 25 '13 at 5:13
  • 1
    $\begingroup$ What? How do you conclude that "the NB and quasi-Poisson are estimating no overdispersion"? That's not what the output says -- the NB has a $\theta$ parameter of 0.29 and the quasi-poisson has an estimated dispersion of 26. Actually the way you worded that suggests you're focused on the ratio of unconditional variance to overall mean, which doesn't tell you about conditional overdispersion. $\endgroup$ – Glen_b Jun 25 '13 at 22:42
  • 1
    $\begingroup$ Thank you. I will read this and come back if I still have questions. $\endgroup$ – half-pass Jun 26 '13 at 16:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.