3
$\begingroup$

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.

$\endgroup$
3
  • 1
    $\begingroup$ um .. what's a BACI? $\endgroup$
    – Glen_b
    Nov 11 '15 at 22:20
  • $\begingroup$ before-after-control-impact. It's a common term in ecology/environmental impact assessment. $\endgroup$
    – Ben Bolker
    Nov 12 '15 at 0:11
  • 1
    $\begingroup$ Thanks for the question - I've edited the title to make this more evident. $\endgroup$
    – abmiller8
    Nov 12 '15 at 13:41
2
$\begingroup$

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.

$\endgroup$
5
  • 1
    $\begingroup$ 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? $\endgroup$
    – abmiller8
    Nov 13 '15 at 22:13
  • $\begingroup$ 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. $\endgroup$
    – Russ Lenth
    Nov 13 '15 at 23:36
  • $\begingroup$ 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! $\endgroup$
    – abmiller8
    Nov 14 '15 at 15:44
  • $\begingroup$ 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. $\endgroup$
    – Russ Lenth
    Nov 16 '15 at 15:52
  • $\begingroup$ Thanks, I'll email the developer. I am also trying to run the baci contrasts on results from an occupancy analysis performed using Rmark. Do you know if there's a way to calculate lsmeans from occupancy results to complete the contrasts analysis on this type of model? $\endgroup$
    – abmiller8
    Nov 16 '15 at 16:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.