6
$\begingroup$

I have a set of count data that I have fitted both a linear model to and a Poisson generalized linear model. The mean of the raw data is 233.375 and the standard deviation 279.983. I have been surprised when I fit a Poisson glm in R that the intercept = 5.53865. A negative bionomial glm in R gives an intercept of 5.4526.

Is this because the data are sparse and/or over-dispersed?

The data and R code are:

data <-c(2, 25,1121,361,251,123,123,81,25,215,4,196,0,353,968,336,179,229,92,204,35,299,8,371)
lm(data ~ 1)
glm(data ~ 1, family = poisson)
glm.nb(data ~1, link = log)
Improve this question
$\endgroup$
3
  • 8
    $\begingroup$ Note that $e^{5.53} = 232.75$ That's why the Poisson and NB give different intercepts: They aren't on the same scale. The Poisson model is $\log(E(Y|X)) = \theta'X$. The intercepts are very close. $\endgroup$
    – Peter Flom
    Commented Aug 30, 2012 at 17:35
  • $\begingroup$ For the linear model are you doing ordinary linear regression? If so that ignores the fact that the response is a count variable. The mean and standard deviation are close to equal the variance is a lot larger than the mean and so the negative binomial is probably a better model due to overdispersion. $\endgroup$
    – Michael R. Chernick
    Commented Aug 30, 2012 at 17:37
  • $\begingroup$ The linear model is regular linear regression. My plan is to fit a negative binomial but I'd like to interpret the parameters in terms of the biology of the study. The counts are of tree saplings from 6 10m^2 plots that received 2 treatments. The plots were surveyed in 2004 and 2011. (6 plots * 2 treatments * 2 years = 24observations). I had thought that the Poisson/Neg-bin intercept could be interpreted as the mean number of saplings per plot, but ~5 is way too small, so I think I don't understand what the model is doing. $\endgroup$
    – N Brouwer
    Commented Aug 30, 2012 at 19:04

1 Answer 1

Reset to default
1
$\begingroup$

The difference in estimated intercepts is not because of the overdispersion in the data. Peter Flom's comment is the correct answer. To see this, change the lm() model into a glm() model with a gaussian family:

glm(data ~ 1, family = gaussian)
glm(data ~ 1, family = gaussian(link="log"),start=c(20))

The canonical link for the gaussian family is the identity link, so you get exactly the same estimate as for lm(). Changing the link to the log link function gives you the same estimate of the intercept that you're getting from the Poisson and NB models. The gaussian model with log link is $log(E(Y|X))=θ^′X$, while the glm with identity link is $E(Y|X)=θ^′X$. That's why exponentiating the estimated intercept for the log link models $e^{5.453} = 233$ gives you the estimated intercept for the identity link models - you are using the inverse link function. Getting the expected # of saplings per plot for this simple model is easy with just the value of the coefficient, but once you add treatment effects and other covariates it will be more difficult. You should use the predict() function like this:

data = data.frame(saplings=data,
                  treat=gl(2,6,24,label=c("control","treat")),
                  year=gl(2,12,label=c("2004","2011")))

test.glm = glm.nb(saplings~treat*year,link=log,data=data)
nd = data.frame(treat=gl(2,1,4,label=c("control","treat")),
            year=gl(2,2,label=c("2004","2011")))
predict(test.glm,newdata=nd,type="response")

Also see this question, and read Chapter 6 of Zuur et al (2007) "Analysing Ecological Data"

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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