3
$\begingroup$

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?

$\endgroup$
2
$\begingroup$

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
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Sri Priya Sep 8 '18 at 2:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.