I am fitting a random effects model to some ecological count data (abundance). The aim is to analyse the effect of influence of chemical pest control on insect population. I run the analysis on cumulative abundance data as treatment could not be applied on the same day in all blocks, but it was applied block-wise. I used a balanced randomized block design with the following variables:
- abundance = response variable, continuous (count data)
- treatment = categorical predictor, levels = A, B, Control
- density = pre-treatment pest density, categorical predictor, levels = high, low
- block = the random term, 10 levels
I have fitted the following model:
m01 <- glmer(abundance ~ treatment * density + (1|block), family = poisson)
I get the following diagnostic plots for m01 (I checked residual plots, normal qqplot and actual vs fitted):
Although this looks good, I get rather suspicious p-values out of the summary
> summary(m01) Random effects: Groups Name Variance Std.Dev. block (Intercept) 0.2323 0.482 Number of obs: 60, groups: block, 10 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) 5.18330 0.15412 33.63 < 2e-16 *** treatmentA 0.38389 0.02911 13.19 < 2e-16 *** treatmentB 0.84773 0.02683 31.59 < 2e-16 *** densityLow -0.57396 0.03740 -15.34 < 2e-16 *** treatmentA:densityLow 0.16299 0.04754 3.43 0.000608 *** treatmentB:densityLow 0.03709 0.04455 0.83 0.405027
The response seems actually overdispersed:
> aods3::gof(m01) D = 3106.982, df = 53, P(>D) = 0 X2 = 3093.423, df = 53, P(>X2) = 0
I tried to correct for overdispersion by using observation-level random effect (OLRE) as often recommended:
m02 <- glmer(abundance ~ treatment * density + (1|block) +(1|ID), family = poisson)
> summary(m02) Random effects: Groups Name Variance Std.Dev. ID (Intercept) 0.2563 0.5063 block (Intercept) 0.1582 0.3978 Number of obs: 60, groups: ID, 60; block, 10 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) 5.06166 0.20559 24.620 < 2e-16 *** treatmentA 0.36948 0.22938 1.611 0.107226 treatmentB 0.86452 0.22896 3.776 0.000159 *** densityLow -0.41699 0.23042 -1.810 0.070350 . treatmentA:densityLow 0.08149 0.32496 0.251 0.801991 treatmentB:densityLow -0.16620 0.32451 -0.512 0.608540
> gof(m02) D = 1.614, df = 52, P(>D) = 1 X2 = 1.5875, df = 52, P(>X2) = 1
It seems that now the response is underdispersed! Although p-values seem much more realistic and consistent with my data when I plot means and se, the diagnostic plots seems far from satisfying (the actual versus fitted in particular looks very weird..).
So my questions are:
- Is this a known issue from using OLRE to correct for overdispersion? I am rather new to GLMMs and I would like to understand what is going on!
- Would there be a way to deal with this problem? I used negative binomial GLMMs and LMM (by log-transforming the response) to address the problem but I would like to avoid using those if possible ... GLMM nb seems appropriate but it does not give me much options for statistical inference (reviewers often ask for p-values) and I am not a fan of transforming and back-transforming count data. So I would like to stick to poisson GLMM provided that I can word around this issue.