A simple Bayesian occupancy model for two interacting species with R and jags

We begin by simulating survey data for two commensal species for the following scenario. Species 1 has a 50% probability occupancy. Species 2 benefits from the presence of species 1, so its probability of occupancy is 83% when species 1 is also present, but only 50% in the absence of species 1. One hundred sites are surveyed 10 times each. The probability of detection for both species is 0.3.

data_generator_commensal_species = function(psi, p, nSites, nReps){ # occupancy prob psi, detection prob p
# true occupancy state
z1 <- rbinom(nSites, 1, psi)
z2 <- rbinom(nSites, 1, psi + z1/3)
# sampling of true occupancy state
y1 <- matrix(NA, nSites, nReps)
y2 <- matrix(NA, nSites, nReps)
for(i in 1:nSites) {
y1[i,] <- rbinom(nReps, 1, z1[i]*p) # detection history for species 1
y2[i,] <- rbinom(nReps, 1, z2[i]*p) # "" species 2
}
list(y1,y2)}
y = data_generator_commensal_species(psi = .5, p = .3, nSites = 100, nReps = 10)

We specify the model structure in jags.

model {
## Priors
a0 ~ dunif(-5, 5)
b0 ~ dunif(-5, 5)
b1 ~ dunif(-5, 5)
p ~ dunif(0, 1)

## Model
# State process
for(i in 1:nSites) {
logit(psi1[i]) <- a0
logit(psi2[i]) <- b0 + b1 * Z1[i]
Z1[i] ~ dbern(psi1[i])
Z2[i] ~ dbern(psi2[i])
# Detection process
for(j in 1:nOccs) {
y1[i, j] ~ dbin(p, Z1[i])
y2[i, j] ~ dbin(p, Z2[i])
}
}
}

Then we fit the model with R and jags.

library(rjags)
d = list( y1 = y[], y2 = y[], nSites = 100, nOccs = 10)
mod = jags.model("commensal_multispecies.txt", d)
update(mod)
post = coda.samples(mod, c("a0", "b0", "b1", "p"),1e4)
summary(post)
plot(post)

Here are the posterior distributions. A nice fit to the scenario for the simulated data. 