counts <- c(18, 17, 15, 20, 10, 20, 25, 13, 12)
outcome <- gl(3, 1, 9)
treatment <- gl(3, 3)
glmmod <- glm(counts ~ outcome + treatment, family = poisson)
glmmod$resid


#working
resid(glmmod, type="working")
(counts - glmmod$fitted.values)/exp(glmmod$linear)

#deviance
resid(glmmod, type = "dev")
fit <- exp(glmmod$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(glmmod, type = "resp")
counts - fit

#pearson
resid(glmmod, type = "pear")
(counts - fit)/sqrt(fit)

# Q-Q plot
x11()
plot(glmmod, which = 2)

x11()
qqnorm(rstandard(glmmod, type = "pear"))
qqline(rstandard(glmmod, type = "pear"))