Diagnostic plots for count regression

Question

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:

• Interpreting residual diagnostic plots for glm models?
• Assumptions of generalized linear models
• GLMs - Diagnostics and Which Family

Answer 1

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

Answer 2

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.

Answer 3

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.

Answer 4

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.

Answer 5

I would definitely recommend the {performance} package. It has a check_model(mod1) function that shows the relevant diagnostic plots.

library(MASS)
data(quine) 
  
mod1 <- glm(Days~Age+Sex, data=quine, family="poisson")
summary(mod1)
#> 
#> Call:
#> glm(formula = Days ~ Age + Sex, family = "poisson", data = quine)
#> 
#> Deviance Residuals: 
#>    Min      1Q  Median      3Q     Max  
#> -6.647  -2.964  -1.299   1.465  10.258  
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  2.63090    0.05693  46.215  < 2e-16 ***
#> AgeF1       -0.25232    0.06804  -3.708 0.000209 ***
#> AgeF2        0.35964    0.06083   5.913 3.37e-09 ***
#> AgeF3        0.29915    0.06412   4.665 3.08e-06 ***
#> SexM         0.10476    0.04181   2.506 0.012221 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for poisson family taken to be 1)
#> 
#>     Null deviance: 2073.5  on 145  degrees of freedom
#> Residual deviance: 1908.3  on 141  degrees of freedom
#> AIC: 2506.8
#> 
#> Number of Fisher Scoring iterations: 5

performance::check_model(mod1)
#> Loading required namespace: qqplotr
#> `geom_smooth()` using formula 'y ~ x'
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: ggrepel: 138 unlabeled data points (too many overlaps). Consider
#> increasing max.overlaps

^{Created on 2021-02-05 by the reprex package (v1.0.0)}