SECR in BUGS/JAGS with patchy habitat

Faced with patches of suitable habitat surrounded by inhospitable terrain, or a large extent of habitat punctuated with patches of non-habitat, we had the choice of ML methods or one of the packages designed specifically for Bayesian SECR analysis, such as SPACECAP or SCRbayes. But then we are limited to the range of models provided by package authors: we don't have the flexibility to specify our own models that comes with WinBUGS, OpenBUGS or JAGS.

Here I present a way to incorporate patchy habitat into a BUGS/JAGS model specification.

With uniform habitat

Let's review how a JAGS analysis would work with a large area of uniformly good habitat. (The analysis in WinBUGS or OpenBUGS would be very similar, but the code may need to be modified slightly.) If you aren't familiar with the principles of SECR, take a look here before reading on.

The data consist of capture histories for a number of animals, recording when and where each was captured. Our example will use camera-trap data, so an animal could be caught more than once in a single occasion (eg, one night), and the captures could be in different locations, so we will use a Poisson encounter model (Royle & Gardner 2011). We also know the locations of all the traps, including those that caught nothing.

We designate a state space : in uniform habitat this is simply a rectangle including the traps, and large enough so that any animals with an activity centre (AC) outside the state space have a negligible probability of being captured.

We assume there are more animals in the state space than we caught, and we want to know how many. For this we use data augmentation : to the capture histories we have, we add a large number of all-zero capture histories. Some of these represent real animals we didn't catch, while some are 'imaginary' animals. If we know how many real animals are in the state space, we can calculate the density. For more on data augmentation, see Royle & Dorazio (2012).

A typical JAGS model specification for the analysis is shown below (adapted from WinBUGS code in Royle & Gardner, 2011 p181).

 1 model {
 2   sigma2 ~ dunif(0, 20)      #set up the priors for the parameters
 3   sigma <- sqrt(sigma2 / 2)
 4   lam0 ~ dgamma(0.1, 0.1)
 5   psi ~ dunif(0, 1)
 6
 7   for (i in 1:M){         #loop through the augmented population
 8     z[i] ~ dbern(psi)     #state of individual i (real or imaginary)
 9     SX[i] ~ dunif(xl, xu) #priors for the activity centre for each individual
10     SY[i] ~ dunif(yl, yu) #  xl, yl = lower coordinate; xu, yu = upper value
11
12
13     for(j in 1:J) {       #loop through the J camera trap locations
14       Dsq[i,j] <- pow(SX[i]-trapmat[j,1], 2) + pow(SY[i]-trapmat[j,2],2)
15       g[i,j] <- exp(-Dsq[i,j]/sigma2)
16       lambda[i,j] <- K * g[i,j] * lam0 * z[i]
17       y[i,j] ~ dpois(lambda[i,j])
18     }
19   }
20   N <- sum(z[1:M]) #derive number (check against M)
21   D <- N / A       #derive density
22 }

Let's start with the data, and see how the model says they were generated. (Skip ahead to the next section if you are familiar with this kind of model.)

In line 17, the data are y[i, j] , the number of times that animal i was detected in trap j ; y is a matrix with J columns for the traps and M rows, one for each animal in the augmented population: it has all-zero rows for animals that were not caught.

The y[i, j] are draws from a Poisson distribution with rate parameter lambda[i,j] , which is calculated (line 16) from the number of trapping occasions ( K ), the detection function value ( g[i,j] ), the number of times the animal would be caught if the distance between the trap and the activity centre were zero ( lam0 ) and z[i] . z[i] is an indicator variable, 1 if the animal is real, 0 if it is imaginary: imaginary animals are never caught. Each animal is assigned a value of z in line 8, being 1 with probability psi .

We use a half-normal detection function, and this calculation is done in line 15. Probability of detection depends on the parameter sigma2 and the distance between the trap and the animal's activity centre: the squared distance is Dsq[i, j] . Here sigma2 corresponds to 2σ ² in the usual Gaussian kernel.

The x and y coordinates of the animal's AC are SX[i] and SY[i] , and the differences between the x and y coordinates of the AC and the trap are squared and added to give Dsq[i, j] . trapmat is a matrix with a row for each trap and columns for the x and y coordinates.

Lines 9 and 10 encode the priors for the x and y coordinates of the AC. The variables xl , xu , yl and yu describe the edges of the rectangular state space, so this prior corresponds to a uniform distribution over the whole area.

In lines 2 to 5 we provide minimally-informative priors for the remaining parameters. In lines 20 and 21 we calculate the total number of real animals (ie, those with z[i] = 1) and derive a density by dividing this by the total area of the state space ( A ). After running the model we check that N << M; if not, the data augmentation was inadequate: we need to increase M (add more all-zero rows to the capture history matrix) and rerun.

With patchy habitat

The model above allows ACs to be located anywhere in the entire state space. If part of the space is  unsuitable as habitat for the target species, or if the area of interest is a patch of suitable habitat in a larger area of inhospitable country (see diagram right), this will not work. We need to provide information on which parts of the state space are suitable habitat, and insert two lines of code to ensure that ACs are only located in suitable habitat.

Specifying suitable habitat

In the diagram right we have an isolated patch of habitat delineated by the black line. This is overlaid with a grid of squares (pixels), in this case 50 pixels east-west and 46 north-south.

We create a matrix with 50 rows and 46 columns, where cell [x, y] = 1 if the corresponding pixel is good habitat, 0 otherwise. This matrix is habmat in the code below. In the plot, cells with 1 are red, those with 0, grey. Such a matrix can be generated as a raster layer in a GIS and imported into R for further manipulation. This is the way other SECR programs deal with habitat data.

To streamline the JAGS estimation, we work in units of pixels instead of metres (or kilometres or whatever were the original units), with the SW corner of the state space having coordinates (0, 0). In particular, we convert the trap coordinates to the new units, so that x coordinates are between 0 and 50, and y coordinates between 0 and 46. For example, if the original trap location was (10750, 20750) metres, this is converted to (8.287463, 15.99673) pixels.

Modifying the JAGS code

We allow JAGS to draw AC locations from the uniform distribution as before: lines 9 and 10 in the code block above don't change, except that now xl = yl = 0 and - in our example - xu = 50 and yu = 46. Then we insert two additional lines of code:

 9    SX[i] ~ dunif(0, xu) #priors for the activity centre for each individual
10    SY[i] ~ dunif(0, yu) #  xu, yu are the upper limits of x, y
11    pOK[i] <- habmat[trunc(SX[i]+1), trunc(SY[i]+1)] # habitat check
12    OK[i] ~ dbern(pOK[i]) # OK[i] = 1, the ones trick

In line 11, we see which pixel the proposed AC location falls into, and get the habitat value, 1 or 0, from the habitat matrix described in the previous section. To get the indices into the matrix, we round up the x and y coordinates of the AC; JAGS has no round-up function, so we add one and round down with trunc . For example, if the proposed AC was at (14.069414, 27.99673), we'd look in habmat[15, 28] .

Now we indulge in a bit of trickery: it's referred to as the "ones trick" in the WinBUGS documentation .

We treat the value pulled out of the habitat matrix as the "probability" that the AC is in good habitat (I called it pOK in the code), though in practice the probability is 1 or 0. And we provide "data" showing whether or not the AC is in good habitat; of course they all are, so OK is just a vector of ones.

In line 12, we compare them: if AC i is in bad habitat, pOK[i] = 0 and the probability that OK = 1 is also zero, ie, it's impossible. So this location will be rejected. If it's in good habitat, pOK[i] = 1 and the probability that OK = 1 is now 1; whether the location is accepted or rejected depends on the rest of the model and the AC to trap distances.

Starting values

You will need to provide initial values for the AC locations. If you leave this to JAGS, it will just draw random values from the priors in lines 9 and 10, many ACs will be in impossible locations, and it will throw an error.

A clever way to do this would be to use the average of the capture locations for the animals actually captured, and give random locations in the good-habitat area for the rest. The spatstat package has functions that will do this.

A lazy but effective way is just to choose a location somewhere in the good habitat and use this as the starting value for all the ACs.

Using WinBUGS or OpenBUGS

Line 11 in the code above will throw an error at the syntax-checking stage in WinBUGS or OpenBUGS: "logical function not allowed in integer expression". They don't like having trunc inside the [ ] of habmat . Replace line 11 with the following:

11.1  xdex[i] <- trunc(SX[i]+1)
11.2  ydex[i] <- trunc(SY[i]+1)
11.3  pOK[i] <- habmat[xdex[i], ydex[i]] # habitat check

I tried this code with WinBUGS via the R2WinBUGS package and it worked just as well as JAGS via rjags but somewhat slower. OpenBUGS via R2OpenBUGS did not converge with the short run I tried (which was adequate for JAGS and WinBUGS).

Does it work?

UPDATED: The plot below left shows the simulated ACs for the four animals which were captured. The spread of the points reflects uncertainty about the location of the AC ( not the size of the home range), so we are more certain about the locations for the black and green animals than for the red animal and especially the blue guy, which was only caught in the trap on the corner of the array.

The right hand plot shows simulated locations for real animals that were not caught. These probably aren't inside the array of traps, but could be anywhere else in the habitat patch. (The plot I first posted included phantom animals too, and they can't be caught even if they have ACs inside the trap array!)

I have now run thousands of simulations with known true densities to investigate the effectiveness of different study designs to detect trends in density, and code similar to that shown here gives good results.

Example code

Click here to download a ZIP file with example R code.

References

Royle, J A; R M Dorazio. 2012. Parameter-expanded data augmentation for Bayesian analysis of capture�recapture models. Journal of Ornithology 152:521-537.

Royle, J A; B Gardner. 2011. Hierarchical spatial capture�recapture models for estimating density from trapping arrays. 163�190 in O'Connell, A F, J D Nichols, and K U Karanth, editors. Camera traps in animal ecology: methods and analyses . Springer, New York.