How can Poisson GLM work with non-count data (rate data)? [duplicate]

This question already has an answer here:

My question is related, but not the same as the following question: Fitting a Poisson GLM in R - issues with rates vs. counts

Here's some fake data:

### some fake data
x=c(1:14)
y=c(0,  1,  2,  3,  1,  4,  9, 18, 23, 31, 20, 25, 37, 45)
y_rate <- y / 1000

I'm going to use a Poisson GLM with a log link to predict y_rate:

### model
pois_mdl <- glm(y_rate ~ x, family=poisson(link="log"))
summary(pois_mdl)

Plot the fit:

### plot
plot(x, y_rate)
lines(x, pois_mdl$fitted.values) I am surprised that Poisson glm() allows for non-integer values in the dependent variable. Draws from a Poisson distribution are always integers (regardless of the value of the mean parameter). Why doesn't glm() blow up? marked as duplicate by amoeba, mdewey, Michael Chernick, kjetil b halvorsen, Matt KrauseNov 13 '17 at 21:55 This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question. 2 Answers I don't know why glm() doesn't blow up. To figure that out, you'll have to unpack all of the underlying code. (In addition, if your only question is how the R code works, this question is off topic here.) What I can say is that you are not modeling the rates correctly. If you want to model rates instead of counts, you need to include an offset in the model's formula. (There is a nice discussion on CV of what an offset is here: When to use an offset in a Poisson regression?) Using your example the code would be: pois_mdl2 <- glm(y~x+offset(log(rep(1000,14))), family=poisson(link="log")) Note that, although the coefficient estimates are the same, the standard errors are quite different: summary(pois_mdl2)$coefficients
#               Estimate Std. Error   z value      Pr(>|z|)
# (Intercept) -6.5681214 0.25118701 -26.14833 1.029521e-150
# x            0.2565236 0.02203911  11.63947  2.596237e-31
summary(pois_mdl)$coefficients # Estimate Std. Error z value Pr(>|z|) # (Intercept) -6.5681214 7.9431516 -0.8268911 0.4082988 # x 0.2565236 0.6969324 0.3680753 0.7128171 While I don't recommend looking at the source code for glm for those that wish to preserve their mental health, I looked at the source code to glm. The reason R doesn't blow up seems to be that it just doesn't bother to do the kind of defensive checks it probably should. The main iteratively re-weighted least squares loop works by using the methods attached to a family object of the appropriate type. In this case, that is poisson: > poi <- poisson() > class(poi)  "family" This family object has all that glm needs to fit the model, for example: > poi$linkfun(1)
 0
> poi$linkinv(1)  2.718282 Another, here's the derivative of the inverse link: > poi$mu.eta(1)
 2.718282

The y data comes in on line 258:

dev <- sum(dev.resids(y, mu, weights))

Unfortunately, dev.resids could not care at all whether y is positive integer valued:

> poi$dev.resid(1.5, 1, 1)  0.2163953 So I guess R doesn't blow up because it didn't think to blow up. • +1, lol. There was a reason I didn't actually look at the source code myself ;-). – gung Aug 20 '15 at 2:55 • The check, if there was one, would most likely occur in poisson()$initialize, which only checks that the ys are non-negative. Compare it to binomial()\$initialize, which, along with dealing with the possible one- or two-column specification of y, also checks that they are integer-valued. – Mark Aug 20 '15 at 7:20
• Nice @Mark, I missed that. – Matthew Drury Aug 20 '15 at 14:07