Making a habitat mask for SECR in JAGS

Back to home page I have just [August 2016] put together an R package to automate the conversion of habitat masks produced with Murray Efford's secr package to the format needed to run in JAGS or WinBUGS/OpenBUGS. [Rewritten 28 June 2019 to work with version 0.1.0 of the package. See here for information on the new version.]

In 2014 I blogged about a way of including information about suitable habitat in a JAGS model for spatially explicit capture-recapture (SECR), where proposed activity centre (AC) locations are checked against a habitat matrix and rejected if they fall in unsuitable habitat. This check is performed for every animal in the augmented data set for every iteration of the MCMC chain, which soon runs into millions of checks. So we use bare-bones code in JAGS (and slightly more complicated code in BUGS). The check involves truncating the x and y coordinates of the activity centre and using those as indices into the habitat matrix. But that means fiddly formatting of the habitat matrix and conversion of the trap coordinates to match.

Having seen colleagues struggling to get the habitat matrix right, I decided to do a function to automate this, starting from a mask created with secr::make.mask.

You can install it from Github using the devtools package by opening R and running
packageVersion("makeJAGSmask") # Should be 0.1.0 or later.

You will need several other packages which you may need to update or install:
install.packages(c("abind", "plotrix", "raster", "secr", "rgeos"))

Using secr::make.mask

The makeJAGSmask package has some example data to play with, which we can display with:

# Get the example data and plot them:
str(simSCR, 1)
attach(simSCR) # remember to detach later
MASS::eqscplot(patchDF, type='l')
points(traps, pch=3, col='red')

Note that we have a trap layout with systematic random placement of clusters. Some of the traps are very close to the edge of the polygon defining suitable habitat.

The object defining the polygon and the trap locations must be in metres, not degrees, and must both use the same coordinate reference system, eg, UTM.

Next we use Murray's code to create a habitat mask; we'll plot it and look at the structure:

secrmask <- secr::make.mask(traps, spacing=1000, type='polygon', poly=patchDF)
plot(traps, add=TRUE)

secr::make.mask has a range of useful options. For example, if your traps are close to the edge of a large patch, you should use type = "trapbuffer" to limit the mask to the area close to the traps:

t2 <- traps[1:7, ]
secrmask2 <- secr::make.mask(t2, buffer=1e4, spacing=1000,
  type='trapbuffer', poly=patchDF)
MASS::eqscplot(patchDF, type='l')
plot(secrmask2, add=TRUE)
plot(t2, add=TRUE) # (plot not shown)

The spacing parameter will become the pixel width. This affects the size of the habitat matrix, which doesn't matter much, though you'll probably want to keep it to a few thousand cells. The matrix will be a "pixilated" version of the habitat polygon, so the boundaries will not match exactly, and the match will be better with smaller pixels. This will only affect traps or activity centres near the boundary. Use the cell.overlap argument to control how pixels are allocated near the edge. Here we set spacing = 1000, so the units for sigma will be kilometres.

Converting the mask to a habitat matrix

This is easily done with the convertMask function:

mymask <- convertMask(secrmask, traps)
detach(simSCR) # no longer needed

By default, the mask is plotted together with the traps, and details of traps which appear to be in unsuitable habitat will be displayed. This is not a problem if they are at the edge of the habitat.

mymask is an object of class JAGSmask, which is list containing:

  • habMat : the 0/1 matrix indicating whether the location at {x, y} is in suitable habitat or not, provided {x, y} is given in pixel-widths from a false origin chosen so that the lower-left (southwest) corner of the mask is at {1, 1}.
  • trapMat : a 2-column matrix with the x and y coordinates of the traps, using the same coordinate system as habMat.
  • upperLimit : a vector of length 2 giving the upper limits for x and y in habMat (the lower limits are both 1); this can be used to specify the priors for activity centres in the JAGS code.
  • pixelWidth : scalar, the width of the pixels in the original units (usually metres).
  • area : the area of the cells of suitable habitat in habMat in square metres; use this to calculate density.

The list has attributes with the location of the false origin, the pixel width, and the bounding box of the original mask, which can be used to convert the JAGS output back to metres.

Class JAGSmask has plotting function, invoked with plot(mymask), which allows you to set several of the graphical parameters.

Implementation in the JAGS model

Take a look at the previous post for code for a typical SECR model in JAGS. The full script to run an example is given at the end of this post. Here we'll just deal with the priors for activity centres and the habitat check. Here are the relevant lines, with the line numbers used in the previous post:

 9 S[i, 1] ~ dunif(1, upperLimit[1]) # priors for the activity centre for each individual
10 S[i, 2] ~ dunif(1, upperLimit[2])
11 pOK[i] <- habMat[trunc(S[i, 1]), trunc(S[i, 2])] # habitat check
12 OK[i] ~ dbern(pOK[i]) # OK[i] = 1, the ones trick

Note that we are using S for the activity centre coordinates in pixel units and AC for the coordinates in metres.

upperLimit and habMat are elements in the JAGSmask object and are included in the data passed to JAGS. OK is a vector of ones, which is also included in the data.

When you set priors for sigma (the scale parameter for the detection function), remember that this will be in units of pixel widths.

You will need to provide JAGS with starting values for the S matrix to ensure that all animals begin in good habitat. See the example script below for ways to do this.

Transforming the JAGS output to the original units

The output for the detection function scale parameter - usually denoted sigma - will be in units of pixel width. Convert to metres with code such as result$sims.list$sigma * pixelWidth(mymask).

The AC  locations in the JAGS output will be in habMat units. Convert them back to the original coordinate system with the convertOutput function. For example, your code may look like this:

S <- result$sims.list$S
AC <- convertOutput(S, JAGSmask = mymask)

The plot on the right shows points from the posterior distribution of the ACs of two animals captured and several animals from the augmented portion of the data set, representing animals present but not caught.

Code to run an example

The following code will run a full example and generate the plots shown above. Here we have a small population and don't need to add too many zeros, so it runs quickly. For the example, we have only done a short run; for a real data analysis you would need a lot more iterations.


# Get the habitat polygon and traps data and plot them:
str(simSCR, 1)
MASS::eqscplot(patchDF, type='l')
points(traps, pch=3, col='red')

# Generate a 'secr' mask and plot:
secrmask <- secr::make.mask(traps, spacing=1000, type='polygon', poly=patchDF)
plot(traps, add=TRUE)

# Convert to a 'JAGSmask':
mymask <- convertMask(secrmask, traps)

# Generate a (sparse) population of Activity Centres
if(getRversion() >= '3.6.0')
  RNGkind(sample.kind="Round") # use the old (buggy) RNG
popn <- sim.popn (D=3e-05, core=traps, buffer = 1e5, poly = patchDF)
nrow(popn) # True number of animals in the polygon
MASS::eqscplot(patchDF, type='l')
plot(traps, add=TRUE)
points(popn, pch=17)

# Simulate capture histories
simCH <- sim.capthist(traps, popn, detectfn = 0,
  detectpar = list(g0 = 0.02, sigma = 2500), noccasions = 90)
plot(simCH, add=TRUE, tracks=TRUE, rad=500, icolours=palette())
plot(traps, add=TRUE)
points(popn, pch=17)

# JAGS code for the model vvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
JAGScode <- "
model {
  sigma ~ dunif(0, 10)     # set up the priors for the parameters
  alpha <- 1/(2*sigma^2)
  p0 ~ dbeta(1, 1)
  omega ~ dbeta(0.001, 1)

  for (i in 1:M){             # loop through the augmented population
    w[i] ~ dbern(omega)       # state of individual i (real or imaginary)
    S[i, 1] ~ dunif(1, upperLimit[1]) # priors for the activity centers for each individual
    S[i, 2] ~ dunif(1, upperLimit[2]) # lower x and y coordinates = 1
    pOK[i] <- habMat[trunc(S[i,1]), trunc(S[i,2])] # habitat check
    OK[i] ~ dbern(pOK[i])     # OK[i] = 1, the ones trick
    for(j in 1:nTraps) {      # loop through the camera trap locations
      Dsq[i,j] <- (S[i,1]-trapMat[j,1])^2 + (S[i,2]-trapMat[j,2])^2
      p[i,j] <- p0 * exp(-Dsq[i,j] * alpha)
      y[i,j] ~ dbin(p[i,j] * w[i], nOcc)
  N <- sum(w)   # derive number (check that N < M)
writeLines(JAGScode, "patch.jags")
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

# Organise capture histories:
# simCH is an array, 5 animals x 90 occasions x 89 traps
# Add up number of detections across occasions
yObs <- apply(simCH, c(1,3), sum)

# Data augmentation: add all-zero rows
M <- 25
nTraps <- nrow(traps)
y <- rbind(yObs, matrix(0, M - nrow(yObs), nTraps))

# organise data for JAGS
dataList <- c(mymask, list(y = y, M = M, nTraps = nTraps, OK = rep(1, M), nOcc = 90))

# Initial values - several options:
# (1) Quick and dirty, but works as {15, 15} is in good habitat:
inits <- function() list(w = rep(1, M), S = matrix(15, M, 2))

# (2) A sounder approach is to start all ACs at random locations in good habitat:
inits <- function() list(w = rep(1, M), S = randomPoints(M, mymask))
# ... but all-random AC locations will produce errors if the starting point for a captured
# animal is too far from the place it was captured.

# (3) Use starting values close to traps where the animal was caught.
Scapt <- matrix(NA, nrow(yObs), 2)
for(i in 1:nrow(yObs)) {
  captTraps <- which(yObs[i, ] > 0) # Which traps caught the animal
  captLocs <- mymask$trapMat[captTraps, , drop=FALSE] # Locations of the traps
  Scapt[i, ] <- colMeans(captLocs)
  # Check it's in good habitat (might not be if trap is on edge):
  stopifnot(mymask$habMat[Scapt[i, , drop=FALSE]] == 1)
head(Scapt, 7) # Check
inits <- function() list(w = rep(1, M), S = randomPoints(M, mymask, Scapt))
head(inits()$S, 7)
head(inits()$S, 7) # First 5 unchanged, last 2 different.

# Estimates to extract:
wanted <- c("N", "omega", "p0", "sigma", "S", "w")

# Run it, takes about 5 mins
result <- jags(dataList, inits, wanted, "patch.jags", DIC=FALSE,
  n.chains=3, n.adapt=1000, n.iter=5000, n.burnin=500, n.thin=10,
result # n.eff very small, but ok for this test.
# Check that data augmentation was adequate
max(result$sims.list$N) # Must be less than M

# Extract the chains for X and Y coordinates and w for 10 animals,
#   first 5 captured, next 5 not captured or phantoms.
S <- result$sims.list$S[, 1:10, ]
w <- result$sims.list$w[, 1:10]
# Remove phantom animals
S[w == 0] <- NA # Yes, this works fine, even though S has an extra dim:

# Convert to the original coordinate system
AC <- convertOutput(S, JAGSmask = mymask)

# Do plots with the original coordinates
MASS::eqscplot(patchDF, type='l')
for(i in 1:5)
  points(AC[, i, ], col=i)
points(AC[, 6:10, 1], AC[, 6:10, 2], col='grey')
points(traps, pch=3, col='red')
points(popn, pch=21, bg='white', cex=2) # true AC locations.

detach(simSCR) # clean up
Updated 2 July 2019 by Mike Meredith