I have done a poisson regression on my data set and am now looking to investigate the model fit. I notice that the fitted values from predict() in r give me the pre exp transformed values. For diagnostic purposes, should I use the predict() values or the estimated response(a la type=response) to pair up with deviance, pearson or response residuals?
As explained here, neither deviance nor Pearson residuals are ideal for diagnosing Poisson models, as they will appear visually inhomogeneous for low count rates, even if the model is entirely correct.
Instead, you can use the DHARMa package, which implements the idea of randomized quantile residuals by Dunn and Smyth (1996). Here an example with a missing quadratic predictor in the glm. Note the comments about testing for overdispersion in the vignette.
library(DHARMa) dat = createData(replicates = 1, sampleSize = 300, intercept = 0, fixedEffects = 1, quadraticFixedEffects = 2, randomEffectVariance = 0, family = poisson()) fit = glm(observedResponse ~ Environment1 , data = dat, family = poisson) res = simulateResiduals(fit) plot(res)