# Pearson Residuals in the Poisson GLM

I am quite new to GLMs and have just fit a Poisson Regression in R to model a positive response $$y$$.

I now want to check how well it fits the data. In particular, I am unsure whether the variance really increases linearly with the mean. To check this, I plotted fitted values against the squared residuals ($$\hat \mu_i$$ vs $$r_i^2$$) hoping for a linear relationship.

Alternatively, I checked the plot of fitted values versus Pearson residuals. This should not show any pattern, right?

Does this make sense?

Yes, plotting the Pearson residuals vs fitted is one way to check for deviations from the assumed variance-mean relationship; the point cloud should be roughly equal in width (i.e. no "fanning" in either direction). However, this can be hard to judge if the sample is unbalanced (i.e. lots more points for some ranges of fitted values), because our eyes tend to read the total amplitude of the point cloud (which will be bigger if there are more samples, even if the variance is unchanged).

Your $$\mu_i$$ vs $$r_i^2$$ plot doesn't suffer as much from this problem.

However, the standard solution is the scale-location plot, which typically shows $$\sqrt{|r'_i|}$$, where $$r'$$ denotes the standardized Pearson residuals (i.e. $$(y_i-\mu_i)/\sqrt{\phi V(\mu_i) (1-h_i)}$$, where $$V()$$ is the variance function, $$\phi$$ is the estimated dispersion, and $$h$$ are the "hat values". If you have a fitted glm() model m, this is what plot(m, which=3) shows you (it's the third panel in the default 4-panel plot produced by plot(m)). (You can retrieve the standardized residuals for your model with rstandard(m).) This has some advantages: (1) standardizing takes care of the increased variance expected in values with greater leverage/farther from the center of the data; (2) according to ?plot.lm, plotting the square-root of the absolute values of the residuals reduces skew in the results (I'm not actually sure why that's important). The plot method also adds a smoothed line, which can help you assess the overall trend. (3) It is arguably easier to detect deviations from a constant/flat line than from linearity (e.g. see here for "worm plots", or detrended quantile-quantile plots, as an alternative to the standard Q-Q plot).

Below I show an example using the famous data set on deaths per year from horse kicks in the Prussian army. I use broom::augment() + ggplot2, which may be overkill for your case (it does allow a little more flexibility, in this case with showing replicated points) to generate the scale-location plot.

The scale-location plot only diagnoses a difference in the overall mean-variance relationship from linearity, not whether there is "too much" variance overall (e.g. $$V= \phi\mu$$ with $$\phi>1$$); don't forget to check for such overdispersion ...

You may also be interested in the DHARMa package.

library(ggplot2)
library(broom)
data("prussian", package="pscl")
m1 <- glm(y~year, data=prussian, family=poisson)
aa <- augment(m1)
ggplot(aa, aes(year, sqrt(abs(.std.resid)))) +
stat_sum(alpha=0.5) +   ## show repeated points
geom_smooth(method="loess", formula=y~x) +  ## smooth line
scale_size(range=c(3,10), breaks=c(1,10))   ## cosmetic


• Thank you for this helpful answer. One more question: How do you decide on the mean-variance relationship suitable for your GLM? Should we order the response, group it and check how the (group-wise) mean and variance relate? – Claudio Moneo May 17 at 13:57
• (1) I would generally try to decide a priori (e.g. assume variance equal/proportional to mean for count data); (2) you could try looking at e.g. AIC for models with different specifications (e.g. different parameterizations of the negative binomial, see here; (3) if you do want a general graphical diagnostic, do it on the residuals, i.e. plot abs(residuals(m)) vs fitted(m) (examining marginal variance is largely useless) – Ben Bolker May 17 at 16:05