I am an ecologist that has counted 11 different species of herbivores in about 110 blocks of two habitat types over 18 months, and I am interested in predicting the density across the study area from the model I hope to make. The approach so far is to build a spatially smoothed GLM, accounting for the over-dispersion with either a poisson distribution or the related negative binomial. Using the R software structure, the models looks like this:

herb.pois <-  glm(Spb~F1+X+Y+X2+Y2+Habitat*Month,
offset = lnArea, data = herbivore, family = poisson)

herb.nb <- glm(Spb~F1+X+Y+X2+Y2+Habitat*MonthF,
offset = lnArea, data = herbivore, family = neg.bin(0.23))

F1 is a remotely sensed vegetation variable I am hoping to correlate with animal density of some or all species, the offset lnArea is the log of the area surveyed in each block (all different sizes, produces an animals per square kilometer measure when prediction value of this is set to 0). The dispersion parameter for the negative binomial (NB) model is individually calculated for each species elsewhere.

Invariably the negative binomial model has greater traction as per a likelihood ratio test (LRT), or AIC comparison.

Model 1: [Poisson model]
Model 2: [Negative binomial model]
Resid. Df Resid. Dev Df Deviance            AIC
1      1298     1252.9                        2080.394
2      1300      763.3 -2   489.58            1942.656

However, both model predictions are wildly different, usually the poisson model is very tight (small computed confidence intervals, using the predict functions: se.fit*1.96)

Poisson model - please ignore Y label and blue line, it a representation of the F1 term

and the negative binomial has larger confidence intervals but importantly fewer important terms.

Negative Binomial model

The bottom line is I want to choose and use the model which best represents the changing density over time, retains important variables and has reasonable confidence intervals (CIs): The differences between these two models and respective likelihoods is leaving me unconvinced, and stretching my statistical knowledge. My gut is telling me that the CI's are too wide for the NB, and that more terms should be important. Should I trust the likelihood ratio test (LRT) to choose the best model? Should I be worried that the predicted values disagree so much? Am I missing something - would you recommend exploring other distributions, data transformations, model scenarios etc.? What other information might I try get to feel surer about the model selection?

I have explored zero inflated models for which predictions didn't seem reasonable (and NB seems to be sufficient?), stepwise dropping terms; for some herbivores spatial terms or the remotely sensed variable would drop out and I need to retain Habitat and Month regardless for prediction purposes.


I thought I had best include this as it may influence the above approach, for instance in multi variate is recommended - my next step is to take the same data, treated a little differently to analyse another aspect of the herbivore ecology. This question is a little more philosophical. On top of an interest in density of herbivores, I am interested in group density (equivalent to the rate at which a predator might encounter herbivores who are clumped) which is exactly equivalent to density in solitary individuals, but ranging to un-correlated in dynamic fission-fusion herding herbivores whose group formation doesn't depend on density alone. Imagine the extremes of a single herd of 500 in the wet season, breaking down into 25 groups of 20 animals in the dry season - density is the same, the first model wouldn't do it justice) The best idea I have had so far is run the models the same way with the dependent term equal to the count of groups of that herbivore in that block (replacing what was the count of animals) i.e. groups per square kilometer. This would allow me to predict another grid for each month of group densities.

Does this seem reasonable? Would I better consider a multi-variate model or something else.

  • 1
    $\begingroup$ As a small note, the confidence intervals in se.fit are not on the response scale, but on the linear scale... you have to transform the upper and lower bounds through the inverse link function. See for example stat.columbia.edu/~martin/W2024/R11.pdf, where it's done for logistic regression. $\endgroup$ – Joe May 1 '14 at 17:01

Your Answer

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

Browse other questions tagged or ask your own question.