BACI contrasts - calculating lsmeans from Mark results

I am working on analysis for a before-after-control-impact study and have been using the occupancy model in Rmark to calculate occupancy of ten species in an area. As this is a BACI study, I would like to calculate the BACI contrast to answer the question, "Is the amount of change between before and after phases significantly different for the impact site as compared with the control site?".
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:
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;
Thanks in advance for any suggestions you can offer!
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!