Making a habitat mask for SECR in JAGS

library(githubinstall)
githubinstall("makeJAGSmask")

Note that you will need to have package abind installed to run makeJAGSmask .

Three years ago 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 .

Using secr::make.mask

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

library(makeJAGSmask)
# Get the habitat polygon and traps data and plot them:
data(patch)
data(traps)
MASS::eqscplot(patch, 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=2000, type='polygon', poly=patch)
plot(secrmask)
plot(traps, add=TRUE)
str(secrmask)

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[93:101, ]
secrmask2 <- secr::make.mask(t2, buffer=1e4, spacing=2000,
    type='trapbuffer', poly=patch)
MASS::eqscplot(patch, type='l')
plot(secrmask2, add=TRUE)
plot(t2, add=TRUE)

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.

Converting the mask to a habitat matrix

This is easily done with the convertMask function:

mymask <- convertMask(secrmask, traps)

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.
  
• 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 AC[i, 1] ~ dunif(1, upperLimit[1]) # priors for the activity centre for each individual
10 AC[i, 2] ~ dunif(1, upperLimit[2])
11 pOK[i] <- habMat[trunc(AC[i, 1]), trunc(AC[i, 2])] # habitat check
12 OK[i] ~ dbern(pOK[i]) # OK[i] = 1, the ones trick

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 AC 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:

ACpix <- result$sims.list$AC
ACm <- convertOutput(ACpix, 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 interations.

library(secr)
library(jagsUI)
library(makeJAGSmask)

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

# Generate a 'secr' mask and check that:
secrmask <- secr::make.mask(traps, spacing=2000, type='polygon', poly=patch)
str(secrmask)
plot(secrmask)
plot(traps, add=TRUE)

# Try this
t2 <- traps[93:101, ]
secrmask2 <- secr::make.mask(t2, buffer=1e4, spacing=2000,
type='trapbuffer', poly=patch)
MASS::eqscplot(patch, type='l')
plot(secrmask2, add=TRUE)
plot(t2, add=TRUE)

# Generate 'habMat':
mymask <- convertMask(secrmask, traps)
str(mymask)
image(mymask$habMat)

# Generate a (sparse) population of Activity Centres
set.seed(1)
popn <- sim.popn (D=3e-05, core=traps, buffer = 1e5, poly = patch)
nrow(popn) # Number of animals in the polygon
MASS::eqscplot(patch, 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(secrmask)
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 {
  sigma2 ~ dunif(0, 20) # set up the priors for the parameters
  sigma <- sqrt(sigma2 / 2)
  lam0 ~ dgamma(0.1,0.1)
  omega ~ dunif(0, 1)

  for (i in 1:M){ # loop through the augmented population
    w[i] ~ dbern(omega) # state of individual i (real or imaginary)
    # priors for the activity centers for each individual:
    AC[i, 1] ~ dunif(1, upperLimit[1])
    AC[i, 2] ~ dunif(1, upperLimit[2])
    pOK[i] <- habMat[trunc(AC[i, 1]), trunc(AC[i, 2])] # habitat check
    OK[i] ~ dbern(pOK[i]) # OK[i] = 1, the ones trick
    for(j in 1:J) { # loop through the J camera trap locations
      Dsq[i,j] <- pow(AC[i,1]-trapMat[j,1], 2) + pow(AC[i,2]-trapMat[j,2],2)
      g[i,j] <- exp(-Dsq[i,j]/sigma2)
      lambda[i,j] <- K * g[i,j] * lam0 * w[i]
      y[i,j] ~ dpois(lambda[i,j])
    }
  }
  N <- sum(w[1:M]) # derive number (check that N << M)
  D <- N / (area / 1e6) # derive density
}"
writeLines(JAGScode, "patch_model.txt")
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

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

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

# organise data for JAGS
dataList <- mymask
dataList$y <- y
dataList$M <- M
dataList$J <- J
dataList$OK <- rep(1, M)
dataList$K <- 90
str(dataList)

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

# A sounder approach is to start all ACs at random locations in good habitat:
inits <- function() list(w = rep(1, M), AC = randomPoints(M, mymask))

# All-random AC locations will produce errors if the starting point for a captured
# animal is too far from the place it was captured. In that case, use starting 
# values close to traps where the animal was caught.
ACcapt <- 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
  ACcapt[i, ] <- colMeans(captLocs)
  # Check it's in good habitat (might not be if trap is on edge):
  stopifnot(mymask$habMat[ACcapt[i, , drop=FALSE]] == 1)
}
head(ACcapt, 7) # Check
inits <- function() list(w = rep(1, M), AC = randomPoints(M, mymask, ACcapt))
head(inits()$AC, 7) 
head(inits()$AC, 7) # First 5 unchanged, last 2 different.

# Estimates to extract:
wanted <- c("N", "D", "omega", "lam0", "sigma", "AC", "w")

# Run it, takes about 2 mins; there's also precooked output.
result <- jags(dataList, inits, wanted, "patch_model.txt",
  n.chains=3, n.adapt=1000, n.iter=5000, n.burnin=500, n.thin=10,
  parallel=TRUE)
result # n.eff very small, but ok for this test.
# Check that data augmentation was adequate
max(result$sims.list$N) # Must be clearly less than M
M

# Extract sigma and convert to metres
sigma <- result$sims.list$sigma
plot(density(sigma)) # units are pixel widths
sigmaM <- sigma * pixelWidth(mymask)
plot(density(sigmaM)) # units are metres

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

# To use the precooked output, do this:
# data(JAGSoutput)
# AC <- JAGSoutput

# Do some plots on the habMat scale to check:
plot(AC[, 1, ]) # Animal #1
points(AC[, 1, ], pch=19, col=rep(2:4, each=450)) # 3 chains in colours
plot(AC[, 2, ])
plot(AC[, 6:10, 1], AC[, 6:10, 2]) # real uncaught animals
points(mymask$trapMat, pch=3, col='red') # None inside the trap clusters

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

# Do plots with the original coordinates
MASS::eqscplot(patch, type='l')
points(myresult[, 1, ]) # Animal #1
points(myresult[, 2, ], col=2)
points(myresult[, 6:10, 1], myresult[, 6:10, 2], col='grey')
points(traps, pch=3, col='red')