Can you tell me what is returned by glm$residuals and resid(glm) where glm is a quasipoisson object. e.g. How would I create them using glm$y and glm$linear.predictors.
glm$residuals
n missing unique Mean .05 .10 .25 .50 .75 .90 .95
37715 10042 2174 -0.2574 -2.7538 -2.2661 -1.4480 -0.4381 0.7542 1.9845 2.7749
lowest : -4.243 -3.552 -3.509 -3.481 -3.464
highest: 8.195 8.319 8.592 9.089 9.416
resid(glm)
n missing unique Mean .05 .10 .25
37715 0 2048 -2.727e-10 -1.0000 -1.0000 -0.6276
.50 .75 .90 .95
-0.2080 0.4106 1.1766 1.7333
lowest : -1.0000 -0.8415 -0.8350 -0.8333 -0.8288
highest: 7.2491 7.6110 7.6486 7.9574 10.1932
Calling resid(model) will default to the deviance residuals, whereas model$resid will give you the working residuals. Because of the link function, there is no single definition of what a model residual is. There are the deviance, working, partial, Pearson, and response residuals. Because these only rely on the mean structure (not the variance), the residuals for the quasipoisson and poisson have the same form. You can take a look at the residuals.glm
function for details, but here is an example:
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
glm.D93 <- glm(counts ~ outcome + treatment, family=quasipoisson())
glm.D93$resid
#working
resid(glm.D93,type="working")
(counts - glm.D93$fitted.values)/exp(glm.D93$linear)
#deviance
resid(glm.D93,type="dev")
fit <- exp(glm.D93$linear)
poisson.dev <- function (y, mu)
sqrt(2 * (y * log(ifelse(y == 0, 1, y/mu)) - (y - mu)))
poisson.dev(counts,fit) * ifelse(counts > fit,1,-1)
#response
resid(glm.D93,type="resp")
counts - fit
#pearson
resid(glm.D93,type="pear")
(counts - fit)/sqrt(fit)