poisson-regressionresiduals
Zuur 2013 Beginners Guide to GLM & GLMM suggests validating a Poisson regression by plotting Pearsons residuals against fitted values. Zuur states we shouldn't see the residuals fanning out as fitted values increase, like attached (hand drawn) plot.
But I thought a key characteristic of the Poisson distribution is that variance increases as mean increases. So shouldn't we expect to see increasing variation in the residuals as fitted values increase?
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.
A) "If a poisson regression is used, then Zuur states the standard errors will be wrong." This is similar to the what happens in a linear model if the errors are heteroskedastic: the variance of $\hat{\boldsymbol{\beta}}_{OLS}$ is $\sigma^2(\mathbf{X'X})^{-1}$, where $\sigma^2$ is the supposed homogeneous variance, but if the variance is not homogeneous an adjustment has to be made to the standard errors. The Poisson model restricts the variance to equal the mean (equidispersion). If the variance exceeds the mean (overdispersion), the Poisson standard errors are wrong because their values depend on the estimated variance: if variance is wrong, standard errors are wrong. As you know, if data are heteroskedastic and you apply OLS, you can (should) throw away the wrong standard errors and prefer heteroskedasticity-consistent standard errors, also called White or sandwich standard errors. Youl could compute sandwich standard errors in a Poisson regression too, but a better specified model would be... better (see Zeileis, Kleiber, and Jackman.)
B) "But why does the model need to make assumptions about the distribution of observations around the fitted values?" Every model depends on the assumption that the distribution of observations around the fitted values is optimal! Think about prediction: if the distribution of observations around the fitted values is far from optimal, predictions are unreliable.
C) "Why can't it use information contained in the distribution of the observed residuals to calculate the standard errors?" Standardized residuals are: $$z_i=\frac{y_i-\hat{y}_i}{\hat\sigma_{y_i}}$$ in a Poisson model $\hat\sigma_{y_i}=\sqrt{E(y_i)}$. If the model is true, then the $z_i$'s sould be approximately independent, each with standard deviation 1. If data are overdispersed, then the $z_i$'s will be larger in absolute value because of the extra variation beyond what is predicted under the Poisson model. In such cases standard errors are actually calculated using "information contained in the distribution of the observed residuals", but that distribution is wrong because your model is misspecified: you assume equidispersion but your assumption does not hold.
D) "My best guess is that for every fitted values, there might only be a few observed values, and there is not enough information contained in these few observed values to calculate the standard errors." What is wrong in your figure is that there are many residuals greater than 1 in absolute value, so they cannot be standardized residuals. More observed values may help, but a better specified model would be more effective. Look at Zeileis, Kleiber, and Jackman. They use a dataset with 4406 observations. Look at the fitted vs Person residuals plots for a Poisson and a negative binomial model:
Which one do you prefer?
However, while the negative binomial already improves the fit dramatically, it can in turn be improved by the hurdle and zero-inflated models. See Zeileis, Kleiber, and Jackman.
Best Answer
The distinction is clear as soon as you understand what a Pearson residual is.
You are correct that for a Poisson model, variance increases as mean increases.
As a result, ordinary raw residuals ($r_i=y_i-\hat\mu_i$) should have a spread that increases with fitted values (though not in proportion).
However, Pearson residuals are residuals divided by the square root of the variance according to the model ($r^P_i=\frac{y_i-\hat\mu_i}{\sqrt{\hat\mu_i}}$ for a Poisson model). This means that if the model is correct, the Pearson residuals should have constant spread.