0
$\begingroup$

I have been trying to fit count data that is zero-inflated and overdispersed using generalized linear mixed models. My research led me to the glmmadmb function in the glmmADMB package. I am fitting these models for several plant species and am modeling the effect of my treatment on seed abundances for a given species. I have Site as a random factor modeled as (1|Site). I first tested for overdispersion. I then fit models with poisson, nbinom and nbinom1 distributions and used AICtab to compare AICs and pick the best model. However, when I look at the residual plots, I see that the models do not fit very well (and AICs are very high, ~400). I am really struggling to figure out if I can improve the fit of these models, or if I should just report my results as they are. If anyone has advice on how I can improve these models, I would really appreciate the help.

Here is a link to my data: SR.15.EQ CSV file

And here is the code I am trying for one of the species (note, they all have weird residual plots, this is just on of them):

glmmNB<- glmmadmb(CON_XAL~Treatment+(1|Site), data = SR.15.EQ, zeroInflation = TRUE, family = "nbinom")

glmmpois<- glmmadmb(CON_XAL~Treatment+(1|Site), data = SR.15.EQ, zeroInflation = TRUE, family = "poisson")

AICtab(glmmNB, glmmpois)

summary(glmmNB)

And for residuals (not sure if the smooth spline line is appropriate here, but residuals are weird even without it)

op <- par(mfrow = c(2, 2))  
plot(fitted(glmmNB), residuals(glmmNB, type="pearson"), main= "Residuals vs Fitted", xlab = "Fitted Values", ylab = "Residuals")  
abline(h = 0, lty = 2)  
lines(smooth.spline(fitted(glmmNB), residuals(glmmNB, type="pearson")))  

plot(SR.15.EQ$Treatment, residuals(glmmNB, type="pearson"), main= "Residuals vs Fitted", xlab = "Fitted Values", ylab = "Residuals")  

qqnorm(residuals(glmmNB, type="pearson"))  
par(op)

Here is my model summary and my residual plots: enter image description here

Let me know if you need any more information! Really looking forward to your feedback.

$\endgroup$
1
$\begingroup$

It may just be that your model doesn't explain a lot of the variance in the data, but here are some points to consider. I couldn't access you data btw, but here goes.

Having an AIC of ~400 doesn't mean that your model fits badly. "AIC as a value by itself is meaningless. It derives meaning from comparison with the AIC values of other models" (Symonds & Moussalli, 2011). What's important is the difference between the AIC of your models. There is some debate over whether a minimum difference in AIC between models is useful. I recommend the Symonds & Moussalli paper as an introduction to understanding AIC, particularly the part of AIC weights.

Unless you have some other variables you haven't included I would report the models as they are. Either choosing the model that you think best applies to you data, i.e. zero-inflation and a chosen distribution, OR as you have relatively few models do an all subset AIC approach. Fit a model with and without treatment for all three distributions. Then calculate the Akaike weights for them to indicate the probability of each model being the best model. If treatment is very important then some combination of the models with it in should have much larger weights than the other models.

Looking at your residuals plots, the first one shows some points with high fitted values. First I would colour the points by treatment to see if the high ones are all from the same treatment. Also I would see if predict(glmmNB, type="response") is clearer. I'm not 100% on this but I think fitted defaults to using the random effect whereas predict() uses the population level.

The second plot looks alright although possibly low variance in the third group. It's important to check if unequal variance is affecting your results: either transform the data or drop the treatment just to see what changes.

The third plot is a qq plot for investigating normality which is not necessary for a poisson/negative binomial model.

I would also point out that just because your data is overdispersed doesn't mean it is necessarily zero-inflated. You could check it by comparing models where zeroInflated = T and F.

Assuming that your data are zero-inflated you may want to consider a hurdle model instead (turn off zeroInflation first), which can be fitted in glmmamdb with family="truncpoiss", "truncbinom" or "truncbinom1". In essence use zero-inflated models if you believe some of your zeros are "false" i.e. due to measurement error (Chapter 11 of Zuur et al 2009).

Zuur, A.F., Ieno, E.N., Walker, N.J., Saveliev, A.A. and Smith, G.M., 2009. Mixed effects models and extensions in ecology with R. New York: Springer. 574 p. Symonds, M.R. and Moussalli, A., 2011. A brief guide to model selection, multimodel inference and model averaging in behavioural ecology using Akaike’s information criterion. Behavioral Ecology and Sociobiology, 65(1), pp.13-21.

$\endgroup$
  • $\begingroup$ Could you give the citation details for the paper and book you mention? $\endgroup$ – mdewey Oct 25 '16 at 13:51

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.