rmodelspscl

Any difference between modelobject$residuals and residuals(modelobject) in R package pscl?


I'm trying to build some models using Zero-Inflated Poisson regression using pscl package and after having manipulated the output object which turns to be zeroinfl, I find that doing residuals(fm_zip) is not equal to fm_zip$residuals.

The following is an example of what I'm talking about:

library("pscl")
data("bioChemists", package = "pscl")
fm_zip <- zeroinfl(art ~ . | 1, data = bioChemists)
names(fm_zip)
fm_zip$residuals
residuals(fm_zip)
all.equal(fm_zip$residuals,residuals(fm_zip))
qplot(fm_zip$residuals,residuals(fm_zip))

As you might realize, the results are not equal. I would say that both ways are equivalent but it seems like they're not. Could you explain me what is wrong with this? According to residuals R help, those two alternatives are supposed to return the difference (observed - fitted). By contrast, I did the same with a plain vanilla linear regression and they are equal.

My R version is:

sessionInfo()
R version 3.0.1 (2013-05-16)
Platform: x86_64-w64-mingw32/x64 (64-bit)...

and the package verion is pscl_1.04.4

Any help is appreciated.


Solution

  • To get equal result you should set type to response ( pearson by default)

     all.equal(fm_zip$residuals,residuals(fm_zip,'response'))
    [1] TRUE
    

    From the ?residuals.zeroinfl:

    The residuals method can compute raw residuals (observed - fitted) and Pearson residuals (raw residuals scaled by square root of variance function).

    The perason variance is defined as:

    mu <- predict(fm_zip, type = "count")
    phi <- predict(fm_zip, type = "zero")
    theta1 <- switch(fm_zip$dist, poisson   = 0, 
                                  geometric = 1, 
                                  negbin    = 1/object$theta)
    variance <- fm_zip$fitted.values * (1 + (phi + theta1) * mu)
    

    EDIT Don't hesitate to read the code behind, it is generally a source of learning and you can also avoid many confusions. To get the code behind the S3 method residuals.zeroinfl, you can use something like this :

    getS3method('residuals','zeroinfl')