Back to home page 
I have just 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.
You can install it by opening R and running
install.packages("makeJAGSmask", repos="http://mikemeredith.net/R")
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 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 .
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
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.
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=3e05, 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 allzero 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))
# Allrandom 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')
