10
$\begingroup$

My question reveals my poor understanding of Poisson regression and GLMs in general. Here's some fake data to illustrate my question:

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

Some custom functions to return psuedo-R2:

### functions of pseudo-R2

psuR2 <- function(null.dev, model.dev) { 1 - (model.dev / null.dev)}

predR2 <- function(actuals, predicted) { 1 - (sum((actuals - predicted)^2)) / sum((actuals - mean(actuals))^2)}

Fit four models: OLS, Gaussian GLM with identity link, Poisson GLM with log link, Poisson GLM with identity link

#### OLS MODEL
mdl.ols=lm(y~x)
summary(mdl.ols)
pred.ols = predict(mdl.ols)

summary(mdl.ols)$r.squared
predR2(y, pred.ols)

#### GLM MODEL, family=gaussian(link="identity")
mdl.guass <- glm(y~x, family=gaussian(link="identity"), maxit=500)
summary(mdl.guass)
pred.guass = predict(mdl.guass)

psuR2(mdl.guass$null.deviance, mdl.guass$deviance)
predR2(y, pred.guass)

#### GLM MODEL, family=possion (canonical link)
mdl.poi_log <- glm(y~x, family=poisson(link="log"), maxit=500)
summary(mdl.poi_log)
pred.poi_log= exp(predict(mdl.poi_log))  #transform

psuR2(mdl.poi_log$null.deviance, mdl.poi_log$deviance)
predR2(y, pred.poi_log)

#### GLM MODEL, family=poisson((link="identity")
mdl.poi_id <- glm(y~x, family=poisson(link="identity"), start=c(0.5,0.5), maxit=500)
summary(mdl.poi_id)
pred.poi_id = predict(mdl.poi_id)

psuR2(mdl.poi_id$null.deviance, mdl.poi_id$deviance)
predR2(y, pred.poi_id)

Finally plot the predictions:

#### Plot the Fit
plot(x, y) 
lines(x, pred.ols)
lines(x, pred.guass, col="green")
lines(x,pred.poi_log, col="red")
lines(x,pred.poi_id, col="blue")

I have 2 questions:

  1. It appears that the coefficients and predictions coming out of OLS and Gaussian GLM with identity link are exactly the same. Is this always true?

  2. I'm very surprised that the OLS estimates and predictions are very different from the Poisson GLM with identity link. I thought both methods would try to estimate E(Y|X). What does the likelihood function look like when I use the identity link for Poisson?

$\endgroup$
  • $\begingroup$ Related: stats.stackexchange.com/questions/142338/… $\endgroup$ – kjetil b halvorsen Aug 18 '15 at 7:57
  • $\begingroup$ If you would like to do least squares to approximate the Poisson model with identity link you could also fit a weighted least squares model, mdl.wols=lm(y~x, weights=1/log(y+1.00000000001)) where the log(y+1.00000000001) is then taken as a first estimate of the variance (sqrt(y+1E-10)) also works - the estimates of such models would be very close to that of the Poisson GLM with identity link... $\endgroup$ – Tom Wenseleers Aug 15 '18 at 9:57
12
$\begingroup$
  1. Yes, they're the same thing. MLE for a Gaussian is least squares, so when you do a Gaussian GLM with identity link, you're doing OLS.

  2. a) "I thought both methods would try to estimate E(Y|X)"

    Indeed, they do, but the way that conditional expectation is estimated as a function of the data is not the same. Even if we ignore the distribution (and hence how the data enter the likelihood) and think about the GLM just in terms of mean and variance (as if it were just a weighted regression), the variance of a Poisson increases with the mean, so the relative weights on observations would be different.

    b) "What does the likelihood function look like when I use the identity link for Poisson?"

    $\mathcal{L}(\beta_0,\beta_1) = \prod_i e^{-\lambda_i}\lambda_i^{y_i}/y_i!$

    $\qquad\qquad\,=\exp(\sum_i -\lambda_i+{y_i}\log(\lambda_i)-\log{(y_i!)}\,)\quad$ where $\lambda_i=\beta_0+\beta_1 x_i$

    $\qquad\qquad\,=\exp(\sum_i -(\beta_0+\beta_1 x_i)+{y_i}\log(\beta_0+\beta_1 x_i)-\log{(y_i!)}\,)$

$\endgroup$
  • 4
    $\begingroup$ An elaboration on Glen_b's second point. One story I told myself, which I found quite clarifying, is that as the estimated conditional mean gets larger in the poisson model, the model get more tolerant of data values far away from the conditional mean. Contrast this with the straight linear model, which is uniformly tolerant no matter what the conditional mean is estimated to be. $\endgroup$ – Matthew Drury Aug 18 '15 at 20:48
  • $\begingroup$ @Glen_b, can I ask you to clarify what you said: "hence how the data enter the likelihood". Are you saying that the likelihood of the model fit is different between an OLS and POisson (link = identity), when fitted using MLE?. I.e., If fitting OLS using MLE, do you use the likelihood function for the normal distribution to calculate the likelihood of the fit, versus the likelihood function from the poisson distribution in the latter case? $\endgroup$ – Alex Feb 25 '16 at 5:14
  • 1
    $\begingroup$ @Alex Right; OLS is ML at the Gaussian and Gaussian likelihood is not Poisson likelihood $\endgroup$ – Glen_b Feb 25 '16 at 6:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.