I think there is a problem when computing Pearson residuals, in that they seem
to be computed at the raw residuals divided by the square root of the
corresponding diagonal element of the weight matrix W evaluated at the last step
of the iterative model fitting procedure (IWLS), instead of dividing by the
square root of the model variance. Here's an example of the weird effects it
produces:
me <- c(1:100)/2
y <- rpois(100,me)
yglm <- glm(y~me-1, family=poisson(link="identity"))
plot(fitted(yglm), residuals(yglm,type="response"))
plots the raw residuals... but
plot(fitted(yglm), residuals(yglm,type="pearson"))
plots the raw residuals MULTIPLIED by sqrt(fitted(yglm))... because it
computes the Pearson residuals dividing by the square root of the diagonal
element of the weight matrix W (at the end of the iterative fitting process)
rather than dividing by the sqrt of variance.