2
$\begingroup$

Summary: I am trying to model some count data. I initially attempted to fit a poisson GLM, but diagnostics appear to indicate overdispersion. I have tried several different recommended remedies but these do not appear to improve the fit of the model. How can I model these overdispersed data effectively?

Data: My response variable y is count data varying from 0 to 11. I have two continuous predictor variables x1 and x2 that I have scaled to mean 0 and sd 1. Previous exploration has indicated x1 follows a non-linear relationship with y and is probably best fit with a quadratic term.

Model: I would like to fit the model y ~ x1 + I(x1^2) + x2 in order to test the associations of the two predictor variables with y.

Approach: Initially I fitted a poisson-family GLM. A residual qq plot indicated non-normal residuals and suggested right-skew. I therefore tested for overdispersion, which seemed to be present. I tried three alternative approaches to deal with the overdispersion:

  1. Negative binomial GLM
  2. Observation-level random effects
  3. Quasi-poisson GLM

In each case the residuals still appear to be right-skewed, if anything more extreme than the original model.

Question: have I diagnosed overdispersion correctly? Is there another problem with my model which is leading to the poor fit, and how can I remedy this?

Code:

# Load packages
library(MASS)
library(dplyr)
library(lme4)

# Load data
d <- tibble(y = c(0, 1, 4, 0, 0, 2, 2, 1, 1, 3, 0, 0, 1, 1, 5, 6, 4, 5, 1, 0, 2, 4, 3, 1, 1, 0, 0, 1, 5, 3, 3, 1, 0, 4, 3, 0, 3, 1, 0, 0, 3, 2, 0, 0, 1, 9, 3, 0, 4, 0, 1, 0, 1, 5, 2, 1, 3, 5, 4, 0, 0, 1, 0, 2, 0, 3, 0, 0, 1, 2, 1, 0, 0, 0, 0, 7, 0, 3, 2, 3, 1, 0, 6, 2, 0, 1, 4, 0, 5, 7, 3, 6, 0, 3, 11, 4, 0, 5, 2, 1), 
            x1 = c(-0.219, 1.278, -0.788, -0.788, 0.456, 0.372, 0.393, 0.92, -0.788, 0.983, -0.788, -0.556, -0.788, 0.034, -0.43, 2.164, 0.709, 0.751, -0.43, -0.198, 1.51, 5.032, 0.013, -0.43, 0.013, -0.788, -0.177, -0.43, 0.097, 0.688, 0.034, -0.577, -0.451, -0.198, -0.282, -0.43, 5.285, 0.92, -0.219, -0.788, 0.372, 0.667, 0.034, 1.342, 1.342, -0.156, 0.435, -0.198, -0.177, -0.788, 4.947, -0.556, 0.878, -0.198, 0.372, -0.409, 1.342, 1.278, 1.004, 0.604, -0.219, 0.097, 0.329, -0.577, -0.409, 1.953, 0.329, 0.414, 0.097, 0.161, 0.393, -0.788, -0.535, -0.767, -0.198, -0.43, -0.198, 0.646, 1.806, -0.788, 4.336, 0.962, 0.414, -0.788, -0.788, -0.577, 1.637, 0.983, -0.788, -0.198, -0.788, -0.198, -0.788, 0.393, 1.342, 1.806, -0.788, -0.788, 0.097, 0.161), 
            x2 = c(1.182, -1.088, -1.088, -1.088, 0.047, 1.182, -0.331, 0.533, -1.088, -0.52, 1.182, 1.182, 1.182, 1.182, 1.182, 0.955, 1.182, 1.182, 1.182, -0.331, 0.804, 1.182, -1.088, -1.088, -1.088, 1.182, -1.088, -1.088, 0.425, 1.182, 1.182, -1.088, 1.182, -0.331, 1.182, 1.182, 1.182, -1.088, 1.182, -1.088, -0.634, 0.274, -0.52, -1.088, 1.182, 1.182, 1.182, 1.182, 0.425, -1.088, 0.047, 1.182, -1.088, -1.088, -0.634, 1.182, 1.182, -1.088, 0.804, -0.331, 1.182, 1.182, -0.52, -1.088, 1.182, 1.182, 1.182, -1.088, -1.088, -0.52, -1.088, -1.088, 0.047, 1.182, -1.088, 1.182, -1.088, -0.18, 1.182, -1.088, 1.182, -0.709, 1.182, -1.088, -1.088, 1.182, -0.439, 0.425, 1.182, -0.331, 1.182, 1.182, -1.088, 1.182, 1.182, -0.331, -1.088, 1.182, 0.614, -0.52))

# Try fitting a poisson GLM
GLM1 <- glm(y ~ x1 + I(x1^2) + x2,
            family = "poisson",
            data = d)
# Check distribution of residuals
GLM1.res <- resid(GLM1)
qqnorm(GLM1.res)
qqline(GLM1.res)

# Assess overdispersion
mean(d$y)
var(d$y)
od.fac <- sum(GLM1.res^2) / df.residual(GLM1)
od.fac

# Try fitting negative binomial GLM
GLM2 <- glm.nb(y ~ x1 + I(x1^2) + x2,
            data = d)
GLM2.res <- resid(GLM2)
qqnorm(GLM2.res)
qqline(GLM2.res)

# Try fitting observation-level random effects
d <- mutate(d, ind = 1:n())
GLM3 <- glmer(y ~ x1 + I(x1^2) + x2 + (1|ind),
              family = "poisson",
              data = d2)
GLM3.res <- resid(GLM3)
qqnorm(GLM3.res)
qqline(GLM3.res)

# Try fitting quasi-poisson
GLM4 <- glm(y ~ x1 + I(x1^2) + x2 + (1|ind),
            family = "quasipoisson",
            data = d)
GLM4.res <- resid(GLM4)
qqnorm(GLM4.res)
qqline(GLM4.res)
Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ Either Anscombe or deviance residuals should look closer to normal than say Pearson or working residuals (but also don't expect them to look actually normal). You are getting deviance residuals, so that's fine. You have more skewness than you'd expect from a Poisson; you're right that it's overdispersed, but an overdispersed Poisson cannot change the skewness of the residuals.) .... I think the negative binomial model should be adequate but you can always look at a zero-inflated distribution. $\endgroup$
    – Glen_b
    Commented Apr 2, 2019 at 11:38
  • $\begingroup$ @Glen_b Thanks for your comment. Are you saying that the skewness of residuals I'm seeing here is within acceptable bounds, or that there is some other problem unrelated to overdispersion that I have not yet diagnosed? (Zero-inflation seems to be marginal, according to a vuong test comparing negbin against zinb model). $\endgroup$
    – user2390246
    Commented Apr 2, 2019 at 13:22
  • 2
    $\begingroup$ With a Poisson model, the skewness and the proportion of zeros are too high given the conditional mean (and the variance is considerably too high), but for a negative binomial there's less concern about them (it's not perfect by any means, but somewhat better). I'd caution against being overly focused on hypothesis testing of assumptions; hypothesis tests don't answer the right question, but instead give an inaccurate response to a question you already know the answer to. Such tests lead you to "fix" non problems and relax about things that matter. $\endgroup$
    – Glen_b
    Commented Apr 2, 2019 at 14:41

2 Answers 2

Reset to default
2
$\begingroup$

When there is overdispersion, you should not necessarily expect that including overdispersion in the model will improve fit, you include it mostly to get correct standard errors (and correct hypothesis tests.) So maybe your analysis are just fine, and you need not do anything extra.

But: residuals from glm's can be difficult to interpret (although Poisson regression is not the worst case), a better alternative might be to use simulated residuals as implemented in the R package DHARMa, an example can be found here: Why my residual-fitted plot looks like this? Specifically, there is no normal assumption in Poisson regression, so the usual qqplots against normality is strange. Simulated residuals give you residuals with a known null distribution, so more interpretable plots.

Quasi-poisson GLM

This yield the exact same estimates, fitted values and residuals as the Poisson GLM, the only difference is in the standard errors (as you observed.)

Negative binomial GLM

Still uses the same log link function, so estimates should not be very different from the Poisson case, but not identical. But no surprise that residuals are similar.

Observation-level random effects

Much the same can be said.

Improve this answer
$\endgroup$
1
$\begingroup$

Have you tried the Vuong's test for Poisson vs. Zero-Inflated Poisson? Although, I agree with @Glen_b about relying too heavily on hypothesis testing. You can also try Hurdle Models (both Poisson and Negative Binomial variants). It's a mixture distribution of binary and truncated count models.

However, conditional overdispersion can be observed in case of omitted variables. Do you have a reason to believe that X1 and X2 are the only predictors in the equation?

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks for the suggestions. I have tried Vuong test, although I have since found that this seems not be recommended: see here . Hurdle models I don't think make biological sense in my specific case. With regards to predictors, I have in fact included some others in my full model, but reduced it down for the sake of the example given here. But while I don't think these apply to my particular example, they are definitely useful considerations for similar situations to mine, so +1. $\endgroup$
    – user2390246
    Commented Apr 9, 2019 at 8:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.