4
$\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)
$\endgroup$
  • 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 Aug 30 '12 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 Chernick Aug 30 '12 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 Aug 30 '12 at 19:04
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"

$\endgroup$

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.