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 capturerecapture
(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 barebones
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
devtools::install_github("mikemeredith/makeJAGSmask")
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:
library(makeJAGSmask)
# Get the example data and plot them:
data(simSCR)
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(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[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)
str(mymask)
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
pixelwidths from a false origin chosen so that the
lowerleft (southwest) corner of the mask is at {1, 1}.
trapMat : a
2column 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.
library(secr)
library(jagsUI)
library(makeJAGSmask)
# Get the habitat polygon and traps data and plot them:
data(simSCR)
str(simSCR, 1)
attach(simSCR)
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)
str(secrmask)
plot(secrmask)
plot(traps, add=TRUE)
# Convert to a 'JAGSmask':
mymask < convertMask(secrmask, traps)
str(mymask)
# Generate a (sparse) population of Activity Centres
if(getRversion() >= '3.6.0')
RNGkind(sample.kind="Round") # use the old (buggy) RNG
set.seed(2)
popn < sim.popn (D=3e05, 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(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 {
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:
dim(simCH)
# 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)
dim(yObs)
# Data augmentation: add allzero 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))
str(dataList)
# 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 allrandom 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,
parallel=TRUE)
result # n.eff very small, but ok for this test.
wiqid::diagPlot(result)
# Check that data augmentation was adequate
max(result$sims.list$N) # Must be less than M
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:
dim(S)
dim(w)
# Convert to the original coordinate system
AC < convertOutput(S, JAGSmask = mymask)
str(AC)
# 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
