# Why are Pearsons residuals from a poisson regression so large?

As I understand it, Pearsons residuals are ordinary residuals expressed in standard deviations.

I've ran this poisson regression:

library(ggplot2)

glm_diamonds <- glm(price ~ carat, family = "poisson", data=diamonds)


I then saved the Pearsons residuals and fitted values from the model:

resid <- resid(glm_diamonds, type = "pearson")
fitted <- fitted(glm_diamonds)
df <- data.frame(resid, fitted)


I then plotted the Pearsons residuals against fitted values:

ggplot(df, aes(fitted, resid)) + geom_point() + ylab("Pearsons residuals") + xlab("Fitted values")


It can be seen in the plot that many of residuals are hundreds of units away from zero. If Pearsons residuals are standard deviations, why are some residuals hundreds of units away from zero? Or in other words, why don't the residuals range from about -3 to 3 if they are standard deviations?

-
As far as I can see, the variable price is continuous in this dataset (price in US dollars). This makes the use of Poisson regression questionable. What was your reasoning when choosing Poisson regression? –  COOLSerdash May 17 at 12:01
I randomly picked a built-in R dataset, without any sort of thinking as to an appropriate model. Question is purely regarding size of Pearsons residuals –  luciano May 17 at 12:06
Well, diamonds are usually valued by "the 4 C's" (colour, clarity, cut, carat) - so you have an overdispersion problem without even looking at the data. Also, Poisson has no dispersion parameter, which is like a normal distribution with a fixed variance.. –  probabilityislogic May 17 at 13:00
I know that you didn't use absolute residuals. What I meant is that absolute values of Pearson residuals indicate that the model doesn't fit the data very well. Further, the residuals should only fall between -3 to +3 under the assumption that your model is correct. –  COOLSerdash May 17 at 13:15
@luciano: Do you know the definition of a Pearson residual? See here if not. –  Scortchi May 17 at 13:49

The key point is that the standardization method for Pearson residuals is to divide the difference between observed values $y_i$ and the fitted Poisson mean $\hat\mu_i$ by the theoretical standard deviation implied by that fitted mean:
$$r_i=\frac{y_i - \hat\mu_i}{\sqrt{\hat\mu_i}}$$
So if the model is badly mis-specified the assumed relation $\operatorname{Var} \mu_i=\mu_i$ can be wildly inaccurate: you have over-dispersion as @probabilityislogic says; moreover the fitted means are much too large for high-carat stones, indicating the assumed linear relation between the log mean and carat is too simple.