# Deviance residuals in Poisson GLM

I am learning the concept of residuals in modelling. I performed a Poisson GLM for a 3x3 contingency table and I got the summary of the model. My question is: the deviance residual from the in-built function resid is differing from the residuals calculated manually using the formula 2 * (y * log(ifelse(y == 0, 1, y/mu)) - (y - mu)). Why are the values different? What am I doing wrong here?

Example:

counts=c(0,1,2,2,0,8,8,13,16)
outcome=gl(3,1,9)
treatment=gl(3,3)
mod=glm(counts~outcome+treatment,family=poisson())
resid(mod,type="deviance")
##         1          2          3          4          5          6          7
##-1.0954451  0.1694307  0.3374099  0.0000000 -2.3664319  1.1368934  0.2176802
##         8          9
## 0.7886223 -0.7609947
fitted.values=mod$fitted.values dev= 2*(counts*log(ifelse(counts==0,1,counts/fitted.values))-(counts-fitted.values)) ## 1 2 3 4 5 ## 1.200000e+00 2.870677e-02 1.138454e-01 -2.928413e-23 5.600000e+00 ## 6 7 8 9 ## 1.292527e+00 4.738466e-02 6.219251e-01 5.791129e-01  My guess is that I am giving a value 1 for all zero counts here. If I applied the formula without imputing a value 1 then the result shows infinite. How are the deviance residuals being calculated in glm? ## 2 Answers You're calculating the deviance per point. The deviance residual is$\textrm{sign}(y-\mu) \sqrt{d_i}$. You need to take the square root (after setting values less than zero, like the tiny residual value for observation 4, to zero), then apply the sign of the difference: d.res <- sqrt(pmax(dev,0)) sign(counts-fitted.values)*d.res  • Thank you very much.. But when I use pmax(dev,0) the result will not exactly same as I get in resid function because the negative value of dev is taken to be 0 here right?? How to sort this out? Sep 8, 2018 at 2:29 I had a similar question. If you crack open stats:::residuals.glm you can find the relevant code: d.res <- sqrt(pmax((object$family\$dev.resids)(y, mu, wts), 0))
ifelse(y > mu, d.res, -d.res)