Zero-inflated Poisson regression: how can I calculate contrasts for BACI (Before-after-control-impact) experiment?

Question

I am analyzing count data with a lot of zeros and found that although the data do not fit a poisson glm, they fit the zero-inflated poisson (ZIP) regression significantly better than the standard poisson glm.

This analysis is for a BACI study, in which I have data before, during, and after the treatment, in three zones: control, on-trail (treatment1) and off-trail (treatment2).

I am interested in the difference in change of detection rate of a species between the on-trail (treatment) site vs. control site, before and after the treatment. I have performed contrast on this data to determine this difference with a simple linear regression model (using lm), but I'm unsure how to find this difference using the zero-inflated poisson model.

The results of my ZIP model are here (ZP = Zone+Phase combined into one variable; the "After" phase is called "Open" in the dataset)

Call:
zeroinfl(formula = Deer ~ ZP | ZP, data = zinb)

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-0.6756 -0.5180 -0.4137 -0.1243 14.8998 

Count model coefficients (poisson with log link):
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.076730   0.031034  34.695  < 2e-16 ***
ZPCDuring    0.080611   0.055391   1.455    0.146    
ZPCOpen     -0.696793   0.092432  -7.538 4.76e-14 ***
ZPFBefore    0.062467   0.042638   1.465    0.143    
ZPFDuring    0.112727   0.067928   1.659    0.097 .  
ZPFOpen     -0.765391   0.080475  -9.511  < 2e-16 ***
ZPTBefore   -0.008428   0.045729  -0.184    0.854    
ZPTDuring   -0.063361   0.063193  -1.003    0.316    
ZPTOpen     -0.717266   0.078422  -9.146  < 2e-16 ***

Zero-inflation model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.40428    0.06617   6.110 9.97e-10 ***
ZPCDuring    0.85610    0.11076   7.729 1.08e-14 ***
ZPCOpen      0.79893    0.13448   5.941 2.84e-09 ***
ZPFBefore   -0.04894    0.09287  -0.527    0.598    
ZPFDuring    1.51339    0.13035  11.610  < 2e-16 ***
ZPFOpen      0.20254    0.12932   1.566    0.117    
ZPTBefore   -0.08598    0.09830  -0.875    0.382    
ZPTDuring    0.98729    0.11720   8.424  < 2e-16 ***
ZPTOpen      0.17416    0.12823   1.358    0.174    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Number of iterations in BFGS optimization: 25 
Log-likelihood: -8572 on 18 Df

and the code i have been using to get the contrasts using the lm model is here (based on this tutorial):

result.lm <- lm(Deer ~ Zone + PhaseBDO + Zone:PhaseBDO, data=counts)
options(contrasts=c(unordered="contr.sum", ordered="contr.poly")) 
result.lsmo.SP <- lsmeans::lsmeans(result.lm, ~Zone:PhaseBDO)

contrast(result.lsmo.SP, list(bact=c(1, 0, -1, 0, 0, 0, -1, 0, 1)))
confint(contrast(result.lsmo.SP, list(bact=c(1, 0, -1, 0, 0, 0, -1, 0, 1))))

Can anyone suggest how I can calculate the contrast from my ZIP regression model to test for significant differences between the difference in detection rates between the two time periods (before/after) in the two zones (control/treatment)?

For example, I want to be able to say if there was a significantly larger increase (or decrease) in the detection rate of deer after the treatment was implemented, in the treatment zone compared with the control zone.

Thanks in advance for your suggestions.

Answer 1

I think you can get what you need by inputting these methods for zeroinfl objects:

recover.data.zeroinfl = lsmeans:::recover.data.lm

lsm.basis.zeroinfl = function(object, trms, xlev, grid, ...) {
    m = model.frame(trms, grid, na.action = na.pass, xlev = xlev)
    contr = object$contrasts
    contr = contr[setdiff(names(contr), c("zero", "count"))]
    X = model.matrix(trms, m, contrasts.arg = contr)
    use = seq_len(ncol(X)) # use only the elts related to response
    bhat = coef(object)[use] 
    V = vcov(object, ...)[use, use]
    nbasis = estimability::all.estble
    misc = list(tran = "log", inv.lbl = "rate")
    dffun = function(k, dfargs) NA
    dfargs = list()
    list(X = X, bhat = bhat, nbasis = nbasis, V = V, 
        dffun = dffun, dfargs = dfargs, misc = misc)
}

Then run the lsmeans and contrasts stuff you had before. This will work for the count portion of the model. Another setup would be needed to use the zero part.

I'll work on supporting zeroinfl objects more fully in lsmeans.