I have been able to do this using count data run with a glm poisson regression model, as well as with the zero-inflated poisson regression model. For the zero-inflated poisson regression model, the developer of the lsmeans package helped me by supplying a script to calculate lsmeans from zero-inflated poisson regression results, which was not already a capability of the lsmeans package.
I imagine others have used occupancy results to answer BACI questions in the past. If anyone has other suggestions for a method I could use to answer this type of question I would be happy to hear that as well. I've tried splitting my data into groups based on site and determining if there is a significant difference between phases (before/after) (see: http://www.phidot.org/forum/viewtopic.php?f=1&t=1906&p=5670&hilit=baci#p5670), but this still leaves the question unanswered of whether there was a difference in the difference between sites. The data were collected over one year, so seasonality is an issue in the analysis. Camera traps were used to collect data.
The model I'm using in Rmark is:
- Code: Select all
species.models=mark(species, model="Occupancy", group = c("Zone", "Phase"), model.parameters = list(p=list(formula=~WeekDay+DetectionDist), Psi=list(formula=~Zone+Phase)), invisible=FALSE)
And the code I'm using in R to calculate the contrasts from zero-inflated poisson regression is below. Perhaps this helps explain my question further. {ZP = Zone+Phase factor;
- Code: Select all
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)
}
bact <- read.csv("C:/....CountCameraDay_AllSp.csv",header=TRUE)
summary(m1 <- zeroinfl(Skunk ~ ZP + ProjectDay + DetectionDist + offset(SurveyLength))|ZP + ProjectDay + DetectionDist, data = bact))
result.lsmo.SP1 <- lsmeans::lsmeans(m1, ~ZP+ProjectDay+DetectionDist)
summary(result.lsmo.SP1)
contrast(result.lsmo.SP1, list(bact=c(1, -1, -1, 1)))
confint(contrast(result.lsmo.SP1, list(bact=c(1, -1, -1, 1))))
Thanks in advance for any suggestions you can offer!