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

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.

• um .. what's a BACI? Nov 11 '15 at 22:20
• before-after-control-impact. It's a common term in ecology/environmental impact assessment. Nov 12 '15 at 0:11
• Thanks for the question - I've edited the title to make this more evident. Nov 12 '15 at 13:41

I think you can get what you need by inputting these methods for zeroinfl objects:
recover.data.zeroinfl = lsmeans:::recover.data.lm
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. • THANK YOU, rvl! This does solve the issue and I'm able to calculate lsmeans and contrasts for the counts portion of the model. Do you have any code that would make the zero part of the model work as well? Nov 13 '15 at 22:13 • All in due time. It is definitely possible. Also, the most useful thing would probably be to combine them into an estimate of the mean with the zeros included. But this will take several days at least. Nov 13 '15 at 23:36 • Thanks, that is great, no hurry. I'm a bit confused what you mean by "combining them into an estimate of the mean with zeros included". Do you mean to take a raw mean of all lines to get the mean detection rate (number of observations/effort)? Or do you mean that it would be possible to combine the "counts" and "zero" portions of the zero-inflaed poisson output, and run lsmeans and contrasts on that combined result? The second (combining counts and zeros from zeroinfl output) would be very interesting and I'd be interested in how to interpret this as well. Thanks again! Nov 14 '15 at 15:44 • I'm talking about the mean of the distribution that is zero with probability$p$, and Poisson with probability$1-p\$. If you send an email to the lsmeans developer, I'll try to keep you apprised of new developments in this type of model. Nov 16 '15 at 15:52