2
$\begingroup$

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):

enter image description here

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:

Model:

m02 <- glmer(abundance ~ treatment * density + (1|block) +(1|ID), family = poisson)

Diagnostic plots: enter image description here

Summary:

> 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 

Overdispersion test:

> 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.
$\endgroup$

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.