# Diagnostic plots for count regression

What diagnostic plots (and perhaps formal tests) do you find most informative for regressions where the outcome is a count variable?

I'm especially interested in Poisson and negative binomial models, as well as zero-inflated and hurdle counterparts of each. Most of the sources I've found simply plot the residuals vs. fitted values without discussion of what these plots "should" look like.

Wisdom and references greatly appreciated. The back story on why I'm asking this, if it's relevant, is my other question.

Related discussions:

Here is what I usually like doing (for illustration I use the overdispersed and not very easily modelled quine data of pupil's days absent from school from MASS):

1. Test and graph the original count data by plotting observed frequencies and fitted frequencies (see chapter 2 in Friendly) which is supported by the vcd package in R in large parts. For example, with goodfit and a rootogram:

library(MASS)
library(vcd)
data(quine)
fit <- goodfit(quine$Days) summary(fit) rootogram(fit)  or with Ord plots which help in identifying which count data model is underlying (e.g., here the slope is positive and the intercept is positive which speaks for a negative binomial distribution): Ord_plot(quine$Days)


or with the "XXXXXXness" plots where XXXXX is the distribution of choice, say Poissoness plot (which speaks against Poisson, try also type="nbinom"):

distplot(quine$Days, type="poisson")  2. Inspect usual goodness-of-fit measures (such as likelihood ratio statistics vs. a null model or similar): mod1 <- glm(Days~Age+Sex, data=quine, family="poisson") summary(mod1) anova(mod1, test="Chisq")  3. Check for over / underdispersion by looking at residual deviance/df or at a formal test statistic (e.g., see this answer). Here we have clearly overdispersion: library(AER) deviance(mod1)/mod1$df.residual
dispersiontest(mod1)

4. Check for influential and leverage points, e.g., with the influencePlot in the car package. Of course here many points are highly influential because Poisson is a bad model:

library(car)
influencePlot(mod1)

5. Check for zero inflation by fitting a count data model and its zeroinflated / hurdle counterpart and compare them (usually with AIC). Here a zero inflated model would fit better than the simple Poisson (again probably due to overdispersion):

library(pscl)
mod2 <- zeroinfl(Days~Age+Sex, data=quine, dist="poisson")
AIC(mod1, mod2)

6. Plot the residuals (raw, deviance or scaled) on the y-axis vs. the (log) predicted values (or the linear predictor) on the x-axis. Here we see some very large residuals and a substantial deviance of the deviance residuals from the normal (speaking against the Poisson; Edit: @FlorianHartig's answer suggests that normality of these residuals is not to be expected so this is not a conclusive clue):

res <- residuals(mod1, type="deviance")
plot(log(predict(mod1)), res)
abline(h=0, lty=2)
qqnorm(res)
qqline(res)

7. If interested, plot a half normal probability plot of residuals by plotting ordered absolute residuals vs. expected normal values Atkinson (1981). A special feature would be to simulate a reference ‘line’ and envelope with simulated / bootstrapped confidence intervals (not shown though):

library(faraway)
halfnorm(residuals(mod1))

8. Diagnostic plots for log linear models for count data (see chapters 7.2 and 7.7 in Friendly's book). Plot predicted vs. observed values perhaps with some interval estimate (I did just for the age groups--here we see again that we are pretty far off with our estimates due to the overdispersion apart, perhaps, in group F3. The pink points are the point prediction $\pm$ one standard error):

plot(Days~Age, data=quine)
prs  <- predict(mod1, type="response", se.fit=TRUE)
pris <- data.frame("pest"=prs[[1]], "lwr"=prs[[1]]-prs[[2]], "upr"=prs[[1]]+prs[[2]])
points(pris$pest ~ quine$Age, col="red")
points(pris$lwr ~ quine$Age, col="pink", pch=19)
points(pris$upr ~ quine$Age, col="pink", pch=19)


This should give you much of the useful information about your analysis and most steps work for all standard count data distributions (e.g., Poisson, Negative Binomial, COM Poisson, Power Laws).

• Great thorough answer! It was helpful to also run through these diagnostics with Poisson-simulated data to train my eye with what the plots should look like. – half-pass Sep 20 '13 at 16:31
• No, this was perfect. You gave lots of references, which was helpful :) – half-pass Sep 20 '13 at 16:35
• Interesting side note: I'm finding that the NB distribution rarely appears to fit simulated NB data based on the GOF test, rootogram, Ord plot, and NB-ness plot. The exception seems to be very "tame" NB data that is nearly symmetric -- high mu, high theta. – half-pass Sep 20 '13 at 17:03
• Hm, are your sure you use type="nbinomial" as argument? E.g. fm <- glm.nb(Days ~ ., data = quine); dummy <- rnegbin(fitted(fm), theta = 4.5) works fine. – Momo Sep 20 '13 at 17:06
• For the life of me, I can't fathom why this didn't get more upvotes. – gung Sep 23 '13 at 1:24

For the approach of using standard diagnostic plots but wanting to know what they should look like, I like the paper:

 Buja, A., Cook, D. Hofmann, H., Lawrence, M. Lee, E.-K., Swayne,
D.F and Wickham, H. (2009) Statistical Inference for exploratory
data analysis and model diagnostics Phil. Trans. R. Soc. A 2009
367, 4361-4383 doi: 10.1098/rsta.2009.0120


One of the approaches mentioned there is to create several simulated datasets where the assumptions of interest are true and create the diagnostic plots for these simulated datasets and also create the diagnostic plot for the real data. put all these plots on the screen at the same time (randomly placing the one based on real data). Now you have a visual reference of what the plots should look like and if the assumptions hold for the real data then that plot should look just like the others (if you cannot tell which is the real data, then the assumptions being tested are likely close enough to true), but if the real data plot looks clearly different from the other, then that means that at least one of the assumptions don't hold. The vis.test function in the TeachingDemos package for R helps implement this as a test.

• An example with the above data, for the record: mod1 <- glm(Days~Age+Sex, data=quine, family="poisson"); if(interactive()) {vis.test(residuals(mod1, type="response"), vt.qqnorm, nrow=5, ncol=5, npage=3)} – half-pass Sep 25 '13 at 1:06

This is an old question, but I thought it would be useful to add that my DHARMa R package (available from CRAN, see here) now provides standardized residuals for GLMs and GLMMs, based on a simulation approach similar to what is suggested by @GregSnow.

From the package description:

The DHARMa package uses a simulation-based approach to create readily interpretable scaled residuals from fitted generalized linear mixed models. Currently supported are all 'merMod' classes from 'lme4' ('lmerMod', 'glmerMod'), 'glm' (including 'negbin' from 'MASS', but excluding quasi-distributions) and 'lm' model classes. Alternatively, externally created simulations, e.g. posterior predictive simulations from Bayesian software such as 'JAGS', 'STAN', or 'BUGS' can be processed as well. The resulting residuals are standardized to values between 0 and 1 and can be interpreted as intuitively as residuals from a linear regression. The package also provides a number of plot and test functions for typical model mispecification problem, such as over/underdispersion, zero-inflation, and spatial / temporal autocorrelation.

@Momo - you may want to update your recommendation 6, it is misleading. Normality of deviance residuals is in general not expected under a Poisson, as explained in the DHARMa vignette or here; and seing deviance residuals (or any other standard residuals) that differ from a straight line in a qqnorm plot is therefore in general no concern at all. The DHARMa package provides a qq plot that is reliable for diagnosing deviations from Poisson or other GLM families. I have created an example that demonstrates the problem with the deviance residuals here.

There is a function called glm.diag.plots in package boot, to generate diagnostic plots for GLMs. What it does:

Makes plot of jackknife deviance residuals against linear predictor, normal scores plots of standardized deviance residuals, plot of approximate Cook statistics against leverage/(1-leverage), and case plot of Cook statistic.