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[[1]], y2 = y[[2]], 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.