Poisson GLM in R – Issues with Rates vs. Counts

generalized linear modelpoisson distributionr

I'm currently working on a project involving GLMs (and eventually GAMs) of some count data over time. Normally I'd do this in SAS, but I'm trying to move to R, and having…issues.

When I fit a GLM to count data using the following:

cdi_model <- glm(counts ~ exposure + covariate + month, data=test, family = poisson)

I get:

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9825  -0.7903  -0.1187   0.5717   1.7649  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.97563    0.20117   9.821  < 2e-16 ***
exposure     0.94528    0.30808   3.068  0.00215 ** 
covariate   -0.01317    0.28044  -0.047  0.96254    
months      -0.03203    0.01303  -2.458  0.01398 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 40.219  on 29  degrees of freedom
Residual deviance: 29.297  on 26  degrees of freedom
AIC: 137.7

Number of Fisher Scoring iterations: 5

Ignore for a moment the performance, or lack thereof of the model itself – mostly playing with syntax and the like at this point.

However, when I try to fit rate data (counts/person-days) and use an offset like so:
cdi_model <- glm(count_rate ~ exposure + covariate + months + offset(log(pd)), data=test, family = poisson)

I get 50+ warnings, all "1: In dpois(y, mu, log = TRUE) : non-integer x = 0.002082" etc. That is more than one for each observation (there's only 30 in the data set).

Additionally, the model fit seems to go to pot. Output as follows:

 Deviance Residuals: 
       Min          1Q      Median          3Q         Max  
-0.0273656  -0.0122169   0.0002396   0.0072269   0.0258643  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept) -15.40110   15.12772  -1.018    0.309
exposure      0.84848   22.18012   0.038    0.969
covariate    -0.02751   21.31262  -0.001    0.999
months       -0.01889    0.95977  -0.020    0.984

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 0.0068690  on 29  degrees of freedom
Residual deviance: 0.0054338  on 26  degrees of freedom
AIC: Inf

Number of Fisher Scoring iterations: 9

Despite this, if I plot the predicted rate against the actual data, the fit doesn't look that much worse, and the actual effect estimate doesn't seem to change all that much.

Anyone have an idea what's going on – or if everything's going right and I'm missing something due to inexperience?

Best Answer

When you add the offset you don't need to (and shouldn't) also compute the rate and include the exposure.

I don't know whether this is the cause of the errors, but if the exposure per case is person days pd, then the dependent variable should be counts and the offset should be log(pd), like this:

cdi_model <- glm(counts ~ covariate + months + offset(log(pd)), 
                 data=test, family = poisson)

Related Solutions

Solved – Understanding over-dispersion as it relates to the Poisson and the Neg. Binomial

This answer comes in two parts. The first addresses the issue of the standard errors and why that implies the models are not identical, as @ndoogan observed in comments, and the second addresses, partially, the issue of when the coefficient estimates might be substantially different.

Consider, for example, a hypothesis test on the coefficient of log(Vmaj) where the null is that the coefficient equals 0.5. The two sets of model estimates would result in rejection of the null in the Poisson case, unless testing at a very high confidence level, and failure to reject the null in the NB case.

More generally, there is more to a collection of estimates than just the point estimates. In the presence of (Negative binomial-like) overdispersion, the standard errors produced by the Negative Binomial-based model will be more accurate estimates of the underlying standard deviation of the coefficient estimates. Thus, the NB model as a whole is more accurate.

For example, a simple model with $\mu_i = \exp\{1 + x_i\}$, and $y_i \sim NB$ with mean $\mu_i$ and variance $5\mu_i$:

N <- 1000
mu <- exp(1 + (x <- rnorm(N)))

p <- 0.2
r <- mu * p / (1-p)
y <- rnbinom(N, r, p)

poisson_summary = summary(glm(y~x, family="poisson"))
nb_summary = summary(glm.nb(y~x))

# Parametric bootstrap calculation of the s.e. of the coefficient of x

coef_x <- rep(0,1000)
for (i in 1:length(coef_x)) {
   mu <- exp(1 + (x <- rnorm(N)))
   r <- mu * p / (1-p)
   y <- rnbinom(N, r, p)
   coef_x[i] <- summary(glm(y~x, family="poisson"))$coefficients["x","Estimate"]
}
data.frame("Sim. SE" = sd(coef_x - 1), 
           "Poisson SE" = poisson_summary$coefficients["x", "Std. Error"],
       "N.B. SE" = nb_summary$coefficients["x", "Std. Error"])

     Sim..SE Poisson.SE    N.B..SE
1 0.03492467 0.01273354 0.03984743

where you can see that the simulated std. deviation of the coefficient estimate is roughly 3x the Poisson-based estimate of same, and much better estimated by the NB-based estimate of same.

The coefficient estimates are pretty much the same, as one might expect with this sample size and the std. errors above, so I won't take up space by displaying them.

Additionally, although this is often a minor effect when overdispersion is low, by weighting the observations with more accurate weights (i.e., better estimates of observation-specific variances), the parameter estimates themselves will be more accurate, well, asymptotically at any rate. The rule of thumb I learned long ago was that heteroskedasticity adjustments (for that is what they are, in essence) don't buy you much unless the differences between weights are on the order of 5x or more.

Note, however, that in small samples you may well get more accurate (in MSE terms) point estimates with the Poisson model if there isn't much overdispersion, because you are reducing the variability induced by estimating the dispersion parameter. Of course, this is almost certainly more than offset by the loss of accuracy in the standard errors of the coefficient estimates.

Solved – How to account for overdispersion in a glm with negative binomial distribution

I'm not sure how to correct the p-values. However you can typically examine the mean-variance assumption in a negative binomial regression by looking at the residuals versus fitted values plot.

If this plot of residuals versus fitted values is not (roughly) an amorphous, random cloud of data points, then you can try using quasi-Poisson regression. Another alternative is to construct your own mean-variance relationship using quasi-likelihood.

Hope this helps!

Best Answer

Related Solutions

Solved – Understanding over-dispersion as it relates to the Poisson and the Neg. Binomial

Solved – How to account for overdispersion in a glm with negative binomial distribution

Related Question