Ok now we're going to put some of these ideas together for a real estimation problem: estimating occupancy where detection in imperfect, but we have replicate samples of each site. The way this works conceptually that each site i is either occupied z[i]=1 or not occupied z[i]=0. If a site is occupied then it is detected or not in each of the j samples. In the simplest single-season occupancy model we have just 2 probabilities: psi the probability of site occupancy, and p the probability of detection. So here's how we proceed:
First we put some priors on the 2 parameters. Standard ones are uniform on (0,1):
psi~dunif(0,1)
p~dunif(0,1)
We then cycle over the sites. For each site z[i] is 1 or 1,
for(i in 1:n.sites)
{
z[i]~dbern(psi)
(the dbern is just a binomial with n.trials=1)
We also go ahead and compute a probability that a sampled site yield a detection, which in this case is the same across replicates so we compute it once outside that loop.
p.eff[i]<-z[i]*p
What this accomplishes is conditioning detection on occupancy. If the site is not occupied, there is no chance of detection (p.eff[i]<-0) . If it is occupied then we have p probability of detection.
We then loop over the replicate samples and model the y[i,j] detection outcomes (0 or 1) as Bernoulli outcomes with probability p.eff[i]. The complete likelihood is then:
for(i in 1:n.sites)
{
z[i]~dbern(psi)
p.eff[i]<-z[i]*p
for (j in 1:n.occas)
{
y[i,j]~dbern(p.eff[i])
}
}
We can also compute a derived quantity, an estimate of the number of sites occupied. By definition this is simply the sum of the z[i] over all sites:
site.occup<-sum(z[])
We then sample from the posterior distribution as before and get posterior statistics on psi, p, and site.occup.
For comparison I've also include a maximum likelihood analysis using the unmarked package in R. You should see that the results agree quite closely, as they should. Note that the Bayesian model is implemented in JAGSg JAGS via the jagsUI in R, which is much faster than OpenBUGs. There are some minor (in this example, practically nonexistent) differences in syntax between BUGS and JAGS. I now tend to use JAGS because it's faster, but it is easier (still not easy) to debug models in BUGS.
Here is the R script to simulate the data, and run the JAGS and unmarked analyses. Pay attention to some of the tricks, notably the initialization of the z vector, which should be 1 when the site is known to be occupied (1 or more detections).
What next
Given the above script and the previous binomial example, you should be well poised to generalize this model. Try adding a random effect over occasions to detection. Or over sites to occupancy or detection. With some care and trial and error (e.g., watch out for duplicate definitions of a node, e.g., do not include an indexed value inside a loop unless 1 or both indices are changing) you will be able to create some interesting and useful models. See also the covariate examples in the Kéry and Schaub chapter.