roffsetglmpscl

Pseudo R² for a Poisson GLM with offset


My question may be of technical nature: I am trying to model disease counts (d) by using population (p) as offset to control for exposure. In R, I found two possible ways to go:

m1 -> glm(d ~ 1 + offset(log(n)), family=poisson, data=dat)
m2 -> glm(d ~ 1, family=poisson, data=dat, offset=log(n))

The summary of m1 and m2 shows that summary(m1) = summary(m2) but if I try to calculate the McFadden through the pR2 (pscl package): McFadden(m1)McFadden(m2).

Does someone have an explanation for that?


Solution

  • Here is the source code of pscl:::pR2.glm:

    function (object, ...) 
    {
        llh <- logLik(object)
        objectNull <- update(object, ~1)
        llhNull <- logLik(objectNull)
        n <- dim(object$model)[1]
        pR2Work(llh, llhNull, n)
    }
    <environment: namespace:pscl>
    

    If the offset is specified in the formula, it gets lost in the second line (update to compute the intercept-only model).

    See this example:

    library("foreign")
    ceb <- read.dta("http://data.princeton.edu/wws509/datasets/ceb.dta")
    ceb$y <- round(ceb$mean*ceb$n, 0)
    ceb$os <- log(ceb$n)  
    
    m0 <- glm(y ~ res + offset(os), data=ceb, family=poisson)
    m1 <- glm(y ~ res, offset=os, data=ceb, family=poisson)
    
    all.equal(coef(m0), coef(m1))
    # [1] TRUE
    
    ### compute null models
    coef(update(m0, ~1))  # wrong, offset not considered
    # (Intercept) 
    #        5.02 
    coef(update(m1, ~1))
    # (Intercept) 
    #       1.376 
    coef(update(m0, ~1, offset=os))
    # (Intercept) 
    #       1.376