1

I want to know how R finds standardized residuals or rstandard when fitting a poisson GLM. I have already found this but it is not talking about standardized residuals. So here is a simple code:

counts <- c(18,17,15,20,10,20,25,13,12)
x=rpois(9,1)
E=counts*10+rpois(9,2)
glm.D93 <- glm(counts ~ x+offset(log(E)), family=poisson)

Now I need to know how R calculates:

> rstandard(glm.D93)
           1            2            3            4            5            6 
 0.018364902 -0.009725560  0.011933387 -0.026455659 -0.036635623 -0.002118836 
           7            8            9 
 0.036552236 -0.022033897  0.034771135

I used help on rstandard to find that there are actually two types of rstandard. One is based on "deviance" (default) and the other on based on "pearson" residuals. Both of these can be easily obtained by the following functions:

> resid(glm.D93,type="dev")
> resid(glm.D93,type="pear")

I am guessing that to find rstandard, I should divide above two residuals by the standard deviation of the ith residual. Thanks.

Community
  • 1
  • 1
Stat
  • 671
  • 7
  • 19
  • Well, you could start with looking at the code of `rstandard.glm` and then go from there. – Roland Nov 04 '13 at 16:48
  • Thanks. I got it. `resid(glm.D93,type="dev")/sqrt(summary(glm.D93)$dispersion * (1 - influence(glm.D93)$hat))` – Stat Nov 04 '13 at 16:59

1 Answers1

3

It's pretty easy to see the code which makes clear the features you deduced. Since rstandard is clearly generic (and most stats package functions are S3, you just add the class to the "stem" function name:

> rstandard.glm
function (model, infl = influence(model, do.coef = FALSE), 
                 type = c("deviance", "pearson"), ...) 
{
    type <- match.arg(type)
    res <- switch(type, pearson = infl$pear.res, infl$dev.res)
    res <- res/sqrt(summary(model)$dispersion * (1 - infl$hat))
    res[is.infinite(res)] <- NaN
    res
}
IRTFM
  • 258,963
  • 21
  • 364
  • 487