Enhancements to 'makeJAGSmask' package

Back to home page I've added extra functions to the 'makeJAGSmask' package since it first appeared in August 2016 - see the old post here for details of the original functions. It's now possible to create a mask from a raster object and to use the values of the raster as a spatial covariate. You can also define a 'core' area smaller than the habitat mask and get an estimate of the number of activity centres in the core.

You can download and install the latest version (currently 0.1.1) of the package with
Note that you will need to have packages abind, plotrix and raster installed to run makeJAGSmask.

The idea of a 'core' area

Often suitable habitat extends outside the study area, and we can't ignore the possibility that animals with activity centres (ACs) outside the study area will appear in our traps. We usually reckon that animals with ACs less than 4 \(\sigma\) from any trap have some probability of capture, and work with a "state space" which includes all suitable habitat within that distance. But the abundance or density of animals in this larger state space is not the figure we are looking for, and it involves extrapolating estimates a long way from the traps, into areas where we have little information other than habitat suitability.

The solution is to produce an estimate of the number or density of ACs in the smaller area we are interested in, the core area. The same methods can be used to estimate populations in regions within the main study area, as in the example below.


We'll demonstrate this using example data from the makeJAGSmask package.

attach(simSCR) # remember to detach later

Take a look at the posterior distributions of the AC locations:

# Convert AC locations to the original units, remove phantoms:
AC <- convertOutput(sims.list$S, JAGSmask)
AC[sims.list$w == 0] <- NA
plotACs(which=1:30, AC=AC, traps=traps, Y=Y, hab=patchDF, show.labels=FALSE)

There is a river flowing through the study area, and the portion north of the river lies in the municipality of Arkadia.

polygon(Arkadia, border='red')

Here Arkadia is a 2-column data frame with the vertices of the polygon forming the boundary, but a SpatialPolygons or SpatialPolygonsDataFrame object imported from a GIS file could also be used.

We add Arkadia to the mask as a core area (yellow in the plot below), then use getACinCore to obtain an MCMC chain for the number, N, of ACs in Arkadia.

coremask <- addCore(JAGSmask, type="poly", poly=Arkadia)
N <- getACinCore(sims.list$S, sims.list$w, coremask)
plot(table(N)) # posterior probability distribution
mean(N)        # posterior mean
quantile(N, probs=c(0.025, 0.975))  # 95% credible interval

Of course it's possible to make more than one JAGSmask object with different core areas and get population estimates for different subsets of the state space.


Using rasters

The new version of makeJAGSmask has a function to generate a 1/0 habitat mask from an object of class RasterLayer or RasterStack. This ensures that the pixels of the 1/0 mask align exactly with those of the raster. The raster itself is converted to a matrix which can be used as a covariate in the analysis in JAGS.

We'll use the example data set above, where simSCR$patchRS is a RasterStack:


The "rand" layer is a random field while "edge" is the distance in metres from the edge of the patch. Non-habitat pixels have NA values.

We will use the distance-from-edge values as a spatial covariate, but first we need to standardise to mean 0, SD 1. To keep the code readable, we'll pull out the "edge" layer:

edge <- patchRS$edge
patchRS$edgeS <- (edge - mean(values(edge), na.rm=TRUE)) / sd(values(edge), na.rm=TRUE)

Now we create our JAGSmask object which will include the raster information as matrices:

mymask <- convertRaster(patchRS, traps)

Note that mymask has three covariate matrices with the original names from the raster stack, but only the first of these appears in the plot. If we had used a single raster layer, the covariate would have been named covMat. With only 5 animals captured (out of 6 total) we can't expect to be able to estimate coefficients properly, but this example will run quickly. We'll only use the edgeS covariate for the analysis below.

In the JAGS model, we prepare a matrix, probs, with the probability that an AC will occur in the corresponding pixel, which is a function of the distance-from-edge covariate and the coefficient beta. Unfortunately, the link function exp is not vectorised, so we have to do that step pixel by pixel. We ensure that the probability of landing in non-habitat is zero, and that the probabilities sum to one.

Candidate AC locations are generated as usual from a uniform prior, we look up the probability of landing in that pixel in the probs matrix, then use the 'zeros trick' (Lunn et al, 2013, §9.5.2, Kéry & Royle, 2016, §5.8) to incorporate this into the likelihood.

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)
  beta ~ dunif(-5, 5)
  # spatial covariate
  for(i in 1:(upperLimit[1]-1)) {
    for(j in 1:(upperLimit[2]-1)) {
      lam[i, j] <- exp(beta * edgeS[i, j]) # 'exp' ensures 'lam' is non-negative
  lam0 <- lam * habMat      # convert 'lam' to 0 for non-habitat
  probs <- lam0 / sum(lam0) # 'probs' must sum to 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]) # uniform priors for the activity centres for each individual
    S[i, 2] ~ dunif(1, upperLimit[2])
    # the zeros trick:
    negLogDen[i] <- -log(probs[trunc(S[i,1]), trunc(S[i,2])])
    zeros[i] ~ dpois(negLogDen[i])  # 'zeros' is a vector of zeros in the data
    # loop through the camera trap locations
    for(j in 1:nTraps) {         
      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)

We save the model in the file "patch_covar.jags", then prepare the other inputs and do a short run:

M <- 25  # size of the augmented data set
nTraps <- nrow(traps)
yAug <- rbind(Y, matrix(0, M - nrow(Y), nTraps))

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

# Initial values
inits <- function() list(w = rep(1, M), S = matrix(15, M, 2))

wanted <- c("N", "omega", "p0", "beta", "sigma", "S", "w")
result <- jags(dataList, inits, wanted, "patch_covar.jags", DIC=FALSE,
  n.chains=3, n.adapt=1000, n.iter=5500, n.burnin=500, n.thin=10, parallel=TRUE)

This takes about 4 mins to run and convergence is good, though n.eff's are small. The randomly-placed ACs did tend to be near the edge as can be seen in the plot in the previous section, so an apparent negative effect of distance-from-edge is expected.

        mean    sd    2.5%    25%    50%    75%  97.5%  Rhat n.eff
N      7.049  2.079  5.000  6.000  7.000  8.000 12.000 1.004   670
omega  0.271  0.115  0.098  0.185  0.258  0.337  0.544 1.002   728
p0     0.023  0.019  0.006  0.012  0.018  0.027  0.080 1.050   220
beta  -0.742  1.032 -3.474 -1.214 -0.544 -0.036  0.712 1.010  1500
sigma  2.922  0.664  1.996  2.432  2.809  3.266  4.586 1.001  1500


Updated 2 July 2019 by Mike Meredith