SCR with an irregular habitat patch

With the stoats example, the detectors were all far from the edges of the suitable stoat habitat, so we used a rectangular state space extending 4 km from the nearest trap. The activity centres (ACs) for the animals we detected could be anywhere within this rectangle. Often this would be unrealistic, a simple rectangle around the detector array would include unsuitable terrain: our study site might be an island surrounded by inhospitable land, or it may just be close to the edge of good habitat.

We need a way to keep ACs within the habitat patch or, more technically, to ensure that the probability of an activity centre being outside the habitat patch is zero. We show how to do that in this unit.

The data

To explore how this works, we’ll use some simulated data from the makeJAGSmask package. These may be data for a species where individuals can be identified from pelage patterns – let’s call them “tigs” – and we used camera traps to collect the data.

The makeJAGSmask package is not on CRAN, but you can install it with the following code: devtools::install_github("mikemeredith/makeJAGSmask")

The habitat patch map

We’ll begin with the habitat defined by a raster, though it could also be defined as a polygon. Let’s access the data set and plot the habitat and the trap locations.


library(raster)
library(makeJAGSmask)
data(simSCR)
str(simSCR, 1)
?simSCR
# Check and plot the habitat raster
simSCR$patchR
# class      : RasterLayer
# dimensions : 59, 65, 3835  (nrow, ncol, ncell)
# resolution : 1000, 1000  (x, y)
# extent     : 124358, 189358, 762937, 821937  (xmin, xmax, ymin, ymax)
# ...
plot(simSCR$patchR)
# Add traps to the plot
points(simSCR$traps, pch=3, col='red')

The raster has pixels 1 x 1 km, and it is 59 pixels wide by 65 pixels tall. The variable shown is the distance from the edge of the patch. The traps are arranged in clusters of up to 8 traps, with clusters systematically distributed over the whole habitat area. Traps were deployed at 89 locations and all worked flawlessly for 90 days.

summary(simSCR$patchR) # The NAs are non-habitat
#             layer
# Min.        0.000
# 1st Qu.  2000.000
# Median   4123.106
# 3rd Qu.  7071.068
# Max.    13601.471
# NA's     1839.000

sum(!is.na(values(simSCR$patchR)))  # each pixel is 1 x 1 km
# [1] 1996

To begin with, we’ll ignore the values of the covariate and just use the fact that pixels outside the habitat area have NA values. There are 1996 pixels which are not NA, so the patch has an area of nearly 2000 sq km.

The detection data

The simulated detection data are in the object simSCR$Y.


dim(simSCR$Y)
# [1]  5 89  # 5 animals captured x 89 detectors
table(simSCR$Y)
#   0   1   3
# 431  13   1 
# How many captures for each animal?
rowSums(simSCR$Y)
# 1 2 3 4 5
# 3 7 1 2 3 
# How many capture locations for each animal?
rowSums(simSCR$Y > 0)
# 1 2 3 4 5
# 3 5 1 2 3

Five animals were detected. One animal (#2) was detected 3 times in the same location, but all except #3 were detected at multiple locations.

Limiting locations to good habitat

Recall how you select random locations in an irregular study area. You draw values for x and y coordinates from a uniform distribution, then check whether the point lies inside the study area. If it’s outside, you reject the candidate point and draw new x and y coordinates. We use the same idea to ensure that ACs do not stray outside the habitat area.

In our JAGS model, we keep the uniform priors for x and y, but we add a step which ensures that points in non-habitat are rejected. We look up the candidate AC location on a “map” to see if the habitat is ok. If it’s not ok, the probability of the AC being there is zero.

The habitat look-up

The map we use is simply the habitat raster converted to a matrix of 1s and 0s, with 1 indicating suitable habitat (called habMat in the code).

The look-up process is performed for each animal in the augmented population for each iteration of the MCMC code: if you have just 50 animals and run 20,000 iterations, that’s 1 million look-ups! So we need to make the look-up process in JAGS really simple, even if it means a bit of work for us before and after the JAGS run.

Instead of working with the “real” x and y coordinates for the ACs, we draw values on the scale of the raster. The raster is 59 pixels wide; we draw values for x between 1 and 60; we then simply truncate the value to find which row of the matrix to check. Similarly, we draw y values between 1 and 66. So if we draw x = 13.45 and y = 28.32, we can check cell [13, 28] of the matrix to see what kind of habitat the point is in. Heres the new JAGS code to do this:

ac[i, 1] ~ dunif(1, upperLimit[1])    # x-coord of activity centers
ac[i, 2] ~ dunif(1, upperLimit[2])    # y coord of activity centers
hab[i] <- habMat[trunc(ac[i, 1]), trunc(ac[i, 2])] # habitat look-up

Of course we need to convert the detector locations to the new raster-based scale before the JAGS run, and convert the AC locations back to the real coordinates afterwards, but functions in the makeJAGSmask package will do that for us.

ACs in non-habitat have zero probability

Code to ensure that the result of the look-up is 1 is a bit tricky: it’s called “the ones trick” in the literature. If our code says A ~ dbern(B) and A = 1, B cannot be 0.

We have seen this before, in the unit about willow tit occupancy with detection covariates. There we had Y[i, j] ~ dbern(p[i, j] * z[i]). Y[i, j] is 1 if we detected the species at site i on occasion j, p[i, j] is the probability of detecting the species at site i on occasion j if it’s present, and z[i] is 1 if the species is present. If Y[i, j] is 1, z[i] cannot be 0; z[i] has to be 1 for all the sites where the species was detected.

We plug hab[i] obtained from the look-up step in as B and plug in a vector of ones (which we supply as data) as A. Here’s the line of code:

ones[i] ~ dbern(hab[i])                # ones[i] = 1, the ones trick

The overall model

The model is very similar to the model for stoats, with the addition of the code to ensure the AC is in good habitat. In the JAGS code below, the new code is highlighted purple. We’re also using ac (lower case) for the coordinates of the activity centre in pixel units. The prior for sigma must be in pixel units; in this case the pixel width is 1 km, so that’s easy.

# Model for proximity detectors and irregular state space
# File name "SCRpatch.jags"

model{
  # Likelihood
  for(i in 1:M){             # Loop over all M individuals
    w[i] ~ dbern(omega)      # w = 1 if animal is real/present
    ac[i, 1] ~ dunif(1, upperLimit[1])    # x-coord of activity centers
    ac[i, 2] ~ dunif(1, upperLimit[2])    # y coord of activity centers
    hab[i] <- habMat[trunc(ac[i, 1]), trunc(ac[i, 2])] # habitat look-up
    ones[i] ~ dbern(hab[i])                # ones[i] = 1, the ones trick

    for(j in 1:nTraps){           # Loop over all detectors
      d2[i,j] <- (ac[i,1] - trapMat[j,1])^2 +
        (ac[i,2] - trapMat[j,2])^2           # distance^2
      p[i,j] <- p0 * exp(- d2[i,j]/(2*sigma^2)) # Detection prob
      Y[i,j] ~ dbin(p[i,j] * w[i], nOcc)     # The observed data
    }
  }

  # Priors
  p0 ~ dbeta(1, 1)      # Baseline detection probability
  sigma ~ dunif(0, 10)  # Half-normal scale
  omega ~ dbeta(1, 1)   # Data augmentation parameter

  # Derived quantities
  N <- sum(w)          # Population size in state-space (=area)
}

Preparing the data

Convert the raster and trap coordinates

We need to convert the raster to the proper matrix for the look-up step in JAGS, and we have to convert the trap coordinates from metres to pixel units. The function convertRaster in the makeJAGSmask package will do that for us:

# Convert the raster to JAGS format
JAGSmask <- convertRaster(simSCR$patchR, simSCR$traps)
str(JAGSmask)
# List of 6
#  $ covMat    : num [1:65, 1:59] 0 0 0 0 0 0 0 0 0 0 ...
#  $ habMat    : num [1:65, 1:59] 0 0 0 0 0 0 0 0 0 0 ...
#  $ trapMat   : num [1:89, 1:2] 6.75 6.75 8.25 9.75 9.75 ...
#   ..- attr(*, "dimnames")=List of 2
#   .. ..$ : chr [1:89] "1.2" "1.3" "1.4" "1.5" ...
#   .. ..$ : chr [1:2] "x" "y"
#  $ upperLimit: Named num [1:2] 66 60
#   ..- attr(*, "names")= chr [1:2] "x" "y"
#  $ pixelWidth: num 1000
#  $ area      : num 2e+09
#  - attr(*, "boundingbox")= num [1:4, 1:2] 124358 189358 189358 ...
#  - attr(*, "origin")= num [1:2] 123358 761937
#  - attr(*, "class")= chr "JAGSmask"

The output from convertRaster contains the habMat matrix we need, and the trap locations in trapMat, plus the upperLimits for x and y, the pixelWidth and the area in square metres, which we can use to calculate the density within JAGS if we wish. The covMat matrix has the values in the original raster for use as a covariate, but we will not use that in this module.

Data augmentation

We augment the data set to 30 individuals by adding 25 extra rows of zeros to the capture history and bundle up the data for JAGS. Note that this includes the elements we need from the JAGSmask list as well as a vector of ones.

# Augment the data
M <- 30                 # size of augmented data set
Nx <- M - nCaps         # number of all-zero capture histories
Yaug <- rbind(simSCR$Y, matrix(0, Nx, nTraps))
w <- c(rep(1, nCaps), rep(NA, Nx)) # 1 for captured animals, NA for others

# Bundle data for JAGS
JAGSdata <- with(JAGSmask, list(habMat=habMat, trapMat=trapMat,
    upperLimit=upperLimit, Y=Yaug, M=M, nOcc=90, nTraps=nTraps,
    w=w, ones=rep(1, M)))
str(JAGSdata)
# List of 9
#  $ habMat    : num [1:65, 1:59] 0 0 0 0 0 0 0 0 0 0 ...
#  $ trapMat   : num [1:89, 1:2] 6.75 6.75 8.25 9.75 9.75 ...
#   ..- attr(*, "dimnames")=List of 2
#   .. ..$ : chr [1:89] "1.2" "1.3" "1.4" "1.5" ...
#   .. ..$ : chr [1:2] "x" "y"
#  $ upperLimit: Named num [1:2] 66 60
#   ..- attr(*, "names")= chr [1:2] "x" "y"
#  $ Y         : num [1:30, 1:89] 0 0 0 0 0 0 0 0 0 0 ...
#   ..- attr(*, "dimnames")=List of 2
#   .. ..$ : chr [1:30] "1" "2" "3" "4" ...
#   .. ..$ : chr [1:89] "1" "2" "3" "4" ...
#  $ M         : num 30
#  $ nOcc      : num 90
#  $ nTraps    : int 89
#  $ w         : num [1:30] 1 1 1 1 1 NA NA NA NA NA ...
#  $ ones      : num [1:30] 1 1 1 1 1 1 1 1 1 1 ...

Run the model fit

Now we specify the parameters we want to monitor and run the model, then check the output.

wanted <- c("sigma", "p0", "omega", "N", "ac", "w")

out <- jags(JAGSdata, NULL, wanted, "SCRpatch.jags", DIC=FALSE,
  n.chains=3, n.adapt=1000, n.burn=0, n.iter=10000, n.thin=10,
  parallel=TRUE) # 7 mins

( mco <- mcmcOutput(out) )
# Object of class 'mcmcOutput'; approx. size 2.27 MB
# ...
# The output has 3 chains each with 1000 draws.
# It has 6 parameters with 94 nodes monitored:
      # nodes
# sigma     1
# p0        1
# omega     1
# N         1
# ac       60
# w        30

# Check augmentation
max(mco$N)
# 29
M
# 30  # just ok

View(summary(mco))
# MCMC values from jagsUI object ‘out’
# The object has 94 nodes with 1000 draws for each of 3 chains.
# l95 and u95 are the limits of a 95% Highest Density Credible Interval.
# Rhat is the estimated potential scale reduction factor:
#         largest is 3.00; 7 (7%) are greater than 1.10; 6 (6%) are NA.
# MCEpc is the Monte Carlo error as a percentage of the posterior SD:
#         largest is 3.1%; NONE are greater than 5%; 5 (5%) are NA.

diagPlot(mco, 1:4)
plot(mco, 1:4)

The plots above were made from the output of a longer run of the model with 10 times the number of iterations and taking 1 hour. These all appear fine. It turns out that the large Rhat values are coming from instances of w with almost exactly 90% zeros, which causes a problem with the non-parametric estimation of Rhat.

Since this is a simulation, we know the true values. The actual population in the habitat patch is 6, and the detection data were simulated with p0 = 0.02 and sigma = 2.5 km. Our estimates are rather high for all those parameters, but they lie well within our credible intervals.

Plot the activity centres

The coordinates in the JAGS output object are on the pixel scale, and we should convert them to metres for plotting. We do that with the convertOutput function in the makeJAGSmask package:

AC <- convertOutput(mco$ac, JAGSmask)
dim(AC)
# [1] 3000   30    2

We have 3000 draws for 30 animals (including many phantoms) with two coordinates, x and y.

For a nice plot, we’ll prepare a raster with a single colour to use as the background, and prepare the data needed to plot the capture locations for each of the animals captured.

# Generate a 1-colour raster for habitat
patch1 <- simSCR$patchR / simSCR$patchR
plot(patch1) # not shown

# Get capture locations for each tig (in metres)
cloc <- vector("list", nCaps)
for(i in 1:nCaps)
  cloc[[i]] <- simSCR$traps[which(simSCR$Y[i, ] > 0), , drop=FALSE]

We’ll do separate plots for each of the 5 animals captured, plus 1 plot for an animal that was not captured. For that one, we only include draws when the animal was a real, uncaught tig, not when it was a phantom.

op <- par(mfrow=c(2,3), mar=c(1,1,1,1))
for(i in 1:5) {
  plot(patch1, legend=FALSE, axes=FALSE, box=FALSE)
  points(AC[, i, ], cex=0.01, col="skyblue")
  points(cloc[[i]], pch=16, cex=2)
  points(simSCR$traps, pch=3, col='red')
  legend('topleft', legend=rownames(simSCR$Y)[i], bty='n', cex=2)
}
real <- mco$w[, 6] == 1
plot(patch1, legend=FALSE, axes=FALSE, box=FALSE)
points(AC[real, 6, ], pch=16, cex=0.5, col="skyblue")
points(simSCR$traps, pch=3, col='red')
legend('topleft', legend="not captured", bty='n', cex=1.5)
par(op)

We have a good idea about the locations of the ACs for the captured animals: remember that the size of the blue blob shows our uncertainty about the AC location, not the size of the animal’s home range. The plots for all uncaptured animals look much the same; they could be anywhere in the habitat patch, but probability is lower in and around each cluster of detectors.

Download a ZIP file with the code here.