Help with zero-inflated generalized linear mixed models with random factor in R

Question

My study has a complicated design and I am not sure if I am modeling my zero-inflated data correctly. I have seed abundances and seedling abundances for 11 species. I have one main "treatment" with four levels (C, P, I, R). For the sampling design: we sampled at 5 sites. At each site, each of the 4 treatments were replicated 4 times (e.g. PA, PB, PC, PD). So there are 4 replicates of 4 treatments nested in 5 sites (so Site is my random factor). The response variable (one species' seed abundances) are zero-inflated (count data). I checked for overdispersion by running a glm and dividing the residual deviance by the degrees of freedom (I believe this is the right approach?). All variables were overdispersed (ratio greater than 1); so I went with negative binomial, instead of Poisson, distributions. I've tried two ways of modeling this: zeroinfl and glmmadmb. But after learning zeroinfl is not useful for mixed models, I am trying to model with glmmADMB:

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

My outcome looks like this:

Summary of my response variable:

 summary(SR.year.raw$CON_XAL)
 Min.   1st Qu.    Median      Mean   3rd Qu.      Max.   
(0.0)     (0.0)     (0.0)   (232.1)     (7.5)  (6245.0)    

My questions are:

1. Am I specifying the random effect correctly in the glmmadmb model?
2. What post-hoc tests can I do with glmmADMB to look at how treatments differ from one another? For example, I have a treatment called "control". But in the model output, there is no significance value for the control. So I suppose the output is saying that the "island" treatment is significantly different from the control. But how do I know if the "Island" treatment is significantly different from the "plantation" treatment?
3. Can you recommend some ways of graphing the results, i.e., something I ultimately could include in my paper? I ask this because from everything I have read, all of the example graphs are for comparing different models (Poisson vs ngbinom), but I can't seem to find code for a good final conclusions graph.
4. Is the warning message a problem? Can I ignore it?

Answer 1

1. No, zeroinfl() currently does not support random effects. So the formula you specified actually means something different: You use a fixed treatment effect in the count part and a fixed site effect in the zero-inflation part. See vignette("countreg", package = "pscl") for more details.

2. If you want random effects, then no. If you use fixed interaction effects instead, you could still try to find a suitable model with zeroinfl(). But with your number of observations this is probably not the best solution.

3. As the model is not the one you would want to fit, this is not relevant here.

4. For zeroinfl() there would be and I suppose that for glmmADMB there are as well. But I'm not an expert on that.

5. You could employ effect plots for the covariate effects or rootograms for the goodness of fit. It depends on what you really want to show.