[R-sig-ME] glmmTMB hurdle model interpretation + post hoc with emmeans

[Guillaume Adeux]

Hi everyone,

I'm trying to analyze the influence of tillage regimes (CT:conventional,
RT:reduced, NT, no-till) on the density of blackgrass (Alopecurus
myosuroides). The thing with this species (as many other weed species), is
that it is often observed in patches (aggregated distribution). Therefore,
even in plots where the weed is present, a lot of quadrats will have 0
values (17 out of 48 observations). Based on these elements, hurdle models
seem like an appropriate way to distinguish presence/absence and density
when present (am I mistaking?).

My experiment was the following:
- 4 blocks
- 3 tillage strips within each block (CT, RT, NT)
- 4 pseudoreplications within each strip

#Reproducing the data frame#
block=c(rep(c(rep("B1",4),rep("B2",4)),3),rep(c(rep("B3",4),rep("B4",4)),3))

tillage=rep(c(rep("CT",8),rep("NT",8),rep("RT",8)),2)
surf_quadrat=as.numeric(as.character(rep(0.36,48)))
ALOMY=as.numeric(c(0,1,1,3,3,1,1,14,0,0,4,0,0,0,0,26,11,1,0,6,2,3,0,11,1,0,0,2,1,0,1,20,4,0,0,1,0,3,34,150,5,0,0,1,3,5,97,135))
weed_density=data.frame(block=block,tillage=tillage,surf_quadrat=surf_quadrat,ALOMY=ALOMY)

#Required packages#
require("glmmTMB")
require("lme4")
require("emmeans")

#Is there truly excess zeroes?#
sum(weed_density$ALOMY<1) #17
sum(weed_density$ALOMY>1) #21
mod_ALOMYpoisson=glmer(ALOMY~block+tillage+offset(log(surf_quadrat))+(1|block:tillage),data=weed_density,family="poisson")
mu <- predict(mod_ALOMYpoisson, type = "response")
exp <- sum(dpois(x = 0, lambda = mu))
round(exp)
    # 7 zeroes are predicted instead of 17

#Comparison between zero inflated and hurdle#
mod_ALOMY_zipoisson0=glmmTMB(ALOMY~block+tillage+offset(log(surf_quadrat))+(1|
block
:tillage),data=weed_density,family=list(family="poisson",link="log"),ziformula=~1)
mod_ALOMY_zipoisson=glmmTMB(ALOMY~ block
+tillage+offset(log(surf_quadrat))+(1| block
:tillage),data=weed_density,family=list(family="poisson",link="log"),ziformula=~
block +tillage)
mod_ALOMY_zineg1=glmmTMB(ALOMY~ block
+tillage+offset(log(surf_quadrat))+(1| block
:tillage),data=weed_density,family=list(family="nbinom1",link="log"),ziformula=~
block +tillage)
mod_ALOMY_zineg2=glmmTMB(ALOMY~ block
+tillage+offset(log(surf_quadrat))+(1| block
:tillage),data=weed_density,family=list(family="nbinom2",link="log"),ziformula=~
block +tillage)
mod_ALOMY_hurdlepois=glmmTMB(ALOMY~ block
+tillage+offset(log(surf_quadrat))+(1| block
:tillage),data=weed_density,family=list(family="truncated_poisson",link="log"),ziformula=~
block +tillage)
mod_ALOMY_hurdleneg2=glmmTMB(ALOMY~ block
+tillage+offset(log(surf_quadrat))+(1| block
:tillage),data=weed_density,family=list(family="truncated_nbinom2",link="log"),ziformula=~block
+tillage)            #truncated_nbinom1 doesn't converge

AIC(mod_ALOMYpoisson,mod_ALOMY_zipoisson0,mod_ALOMY_zipoisson,mod_ALOMY_zineg1,mod_ALOMY_zineg2,mod_ALOMY_hurdlepois,mod_ALOMY_hurdleneg2)
mod_ALOMY_hurdleneg2 # seems to be the best fit

#Significance of tillage effects#
Anova(mod_ALOMY_hurdleneg2,type="III")

#Post hoc analysis#
emmeans(mod_ALOMY_hurdleneg2,~tillage,offset=1,type="response")
pairs(emmeans(mod_ALOMY_hurdleneg2,~tillage)) # pairwise contrasts are on
the verge of significance
cld(emmeans(mod_ALOMY_hurdleneg2,~tillage))

So based on this script, emmeans states that blackgrass density is 7.5
plants/m² in CT, 48 in NT and 40 in RT. The ranking is coherent with the
litterature and field observations.
Is the conclusion this simple or am I overlooking the two components of the
model (conditional, zero inflation). Should any other aspects of the model
be presented for scientific publication (i.e. predictions taking into
account different components of the model)?

Any documentation and/or help is welcome :).

Thank you once again for your help.

Sincerely,

Guillaume ADEUX

	[[alternative HTML version deleted]]