Basic spatial capture-recapture models

We will begin with the simplest SCR model, where the same parameter values apply to all animals and all the detectors, where the detectors do not actually capture the animal, and where the detector array is in the middle of a large area of suitable habitat.

The data

For our example we will use the detection histories, stoatCH, from the stoatsDNA data set in Murray Efford’s secr package. The data are for stoats from 94 hair snags deployed for 7 days – see the details in the help file. You can load and examine the data with the following R code:

library(secr)
?stoatDNA
data(stoatDNA)
plot(stoatCH, rad=50, tracks=TRUE)

The hair tubes (+ in the plot) were laid out in a rectangular grid with minimum spacing 250m. Twenty individual animals were detected, several of them being picked up at multiple locations, as shown by the coloured lines linking detection locations.

Detector types

The formulation of the likelihood depends on the type of detector used:

  • Single-catch traps such as Sherman traps are used for small mammals. Each trap can only capture one animal in one capture occasion, and one animal can only be caught in one trap.
  • Multi-catch traps can catch more than one animal on one occasion, but an animal can only be caught in one trap. Pitfall traps and mist nets are examples.
  • Proximity detectors do not physically catch the animal but detect it as it passes. Camera traps and hair snags are examples. A detector can detect more than one animal on one occasion, and an animal can visit more than one detector. Our example data set uses proximity detectors.

Detection probability

In SCR we make explicit use of the locations where an animal was detected. At least some of the animals must be detected in multiple locations, and we work with a model of animal movement.

Each animal has an activity centre (AC) and mostly moves around close to its AC. The probability of detection depends on the distance between the detector and the AC, being a maximum when the distance is zero and declining monotonically as the distance increases. A variety of detection functions could be used, but a half-normal (half-Gaussian) function is the best in many cases. It is shown below:

This detection function has two parameters:
$p_0$ is the probability of detecting the animal on one occasion if the detector is placed at the activity centre; it will be less than 1, often very small;
$\sigma$ is the scale parameter, determining how detection probability declines with distance; for large animals such as tigers, $\sigma$ will be a couple of kilometers, but for voles less than 10m.

The state space

We need to specify the geographical area within which we will model the ACs. This is known as the state space. The simplest state space is a rectangle around the detector array, which works if suitable habitat extends well beyond the array.

The state space must be large enough to include all the animals which might be detected. With a half-normal detection function, the rule of thumb is that the state space should include a buffer with width $d = 4\sigma$ around the detector array. The probability of detection when that’s the distance between the detector and the AC is $p_0 \times 0.0003$, negligibly small.

Looking at the plot for the stoats data above, we see that animals are captured at locations up to 800m apart. We’ll use a buffer of 1000m and check later to see how this compares with $\sigma$.

In the next block of R code we extract the coordinates of the detectors and convert from m to km and from a data frame to a matrix. Then we define the minimum and maximum coordinates for our state space:

# Extract the detector locations
traps <- traps(stoatCH)
# Convert to a matrix and meters to km
Detlocs <- as.matrix(traps) / 1000
( nDetlocs <- nrow(Detlocs) )
# 94

# Decide on state space with a buffer around the detector array
buffer <- 1  # kilometres
# outer edges of the state space
( xmin <- min(Detlocs[,1]) - buffer )
# -2.5
( xmax <- max(Detlocs[,1]) + buffer )
# 2.5
( ymin <- min(Detlocs[,2]) - buffer )
# -2.5
( ymax <- max(Detlocs[,2]) + buffer )
# 2.5

# Plot it:
plot(Detlocs, xlim=c(xmin, xmax), ylim=c(ymin, ymax),
  bty='n', pch=3, col='red')
rect(xmin, ymin, xmax, ymax, border='red', lwd=3)

The SCR model

As usual, we start reading this from the bottom of the hierarchy. Each data point, $y_{ij}$, is the number of detections of animal $i$ at detector $j$. $n_j$ is the number of occasions when detector $j$ was operating; SCR models can deal with surveys where detectors were not operating for the whole study period, as often happens with camera traps.$$y_{ij} \sim {\rm Binomial} (p_{ij} \cdot w_i, n_j)$$

The probability of detection, $p_{ij} \cdot w_i$, will be zero if the animal is a phantom and $w_i = 0$. Otherwise it is calculated from the base detection probability, $p_0$, the probability of detection when the distance between the detector and the AC is zero, and the Gaussian detection function, $\exp(-d_{ij}^2 / (2\sigma^2))$, where $d_{ij}$ is the distance between detector $j$ and the AC of animal $i$.$$p_{ij} = p_0 \cdot \exp(-d_{ij}^2 / (2\sigma^2))$$

For $p_0$ we use a uniform Beta distribution, though we could use an informative prior if we wished. $$p_0 \sim {\rm Beta}(1, 1)$$

The squared distance, $d_{ij}^2$, between the AC of animal $i$ and detector $j$ is calculated from the differences in the x and y coordinates, the first and second columns of the $AC$ matrix and the matrix of detector locations, $D$.$$d_{ij}^2 = (AC_{i1} – D_{j1})^2 + (AC_{i2} – D_{j2})^2$$

We use uniform priors for the x and the y coordinates of the ACs.$$AC_{i1} \sim {\rm Uniform}(x_{min}, x_{max})$$ $$AC_{i2} \sim {\rm Uniform}(y_{min}, y_{max})$$

On the right hand side of the model diagram, we draw $w_i$ from a Bernoulli distribution with probability $\Omega$ that the animal is real, and we will try a uniform Beta prior for $\Omega$:$$w_i \sim {\rm Bernoulli}(\Omega)$$ $$\Omega \sim {\rm Beta}(1,1)$$

Finally, we derive our estimate of the number of real animals in the state space by summing the $w_i$: $$N = \Sigma{w_i}$$

The effect of using the Bernoulli and Beta distributions in this way, is to induce a discrete uniform prior for $N$:$$N \sim {\rm DiscreteUniform}(0, M)$$where $M$ is the number of rows in the augmented data set. As we discussed on the Data Augmentation page, this might not be the most sensible prior and it might not even converge. As we say there, we can choose $M$ and the parameters for the Beta distribution to implement an informative prior for $N$.

The JAGS code

# Model for proximity detectors and rectangular state space
# File name "SCRrect.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(xmin, xmax) # x-coord of activity centre
    AC[i, 2] ~ dunif(ymin, ymax) # y coord of activity centre
    for(j in 1:nDetlocs){           # Loop over all detectors
      d2[i,j] <- (AC[i,1] - Detlocs[j,1])^2 +
          (AC[i,2] - Detlocs[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, 3)       # Half-normal scale
  omega ~ dbeta(1, 1)       # Data augmentation parameter

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

Preparing the data and running the model

The capture history format for stoatCH is a 3-D 0/1 array, animals x occasions x detectors. We collapse that to an animals x detectors matrix with the number of detections for each animal at each detector. We augment the data with 200 all-zero rows. We know that the first 20 animals are real, so include that in the data as a vector for w with 20 ones and 200 NAs.

# Aggregate to an animals x detectors matrix
#   with number of occasions detected
dim(stoatCH)
# 20  7 94    # animals x occasions x detectors
y <- apply(stoatCH, c(1,3), sum)
dim(y)
# 20 94
table(y)
#    0    1    2
# 1852   26    2 
sum(y) == sum(stoatCH)  # check
# TRUE
# All the detectors were active for all 7 days
nOcc <- 7

# Augment the data
( nCaps <- nrow(y) )  # number of animals captured
# 20
nAug <- 200  # number of rows to add
yAug <- rbind(y, matrix(0, nAug, nDetlocs))
( M <- nrow(yAug) )
# 220

# Bundle data for JAGS
jagsData <- list(y = yAug, M = M, nOcc = nOcc,
    w = c(rep(1, nCaps), rep(NA, M-nCaps)),
    Detlocs = Detlocs, nDetlocs = nDetlocs,
    xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax)
str(jagsData)
# List of 10
# $ y       : num [1:220, 1:94] 0 0 0 0 0 0 0 0 0 0 ...
# $ M       : int 220
# $ nOcc    : num 7
# $ w       : num [1:220] 1 1 1 1 1 1 1 1 1 1 ...
# $ Detlocs : num [1:94, 1:2] -1.5 -1.5 -1.5 ...
# $ nDetlocs: int 94
# $ xmin    : num -2.5
# $ xmax    : num 2.5
# $ ymin    : num -2.5
# $ ymax    : num 2.5

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

library(jagsUI)
# model run takes approx. 45 mins.
out <- jags(jagsData, NULL, wanted, "SCRrect.jags", DIC=FALSE,
    n.chains=3, n.adapt=1000, n.iter=10000, parallel=TRUE)

We’ll use the functions in the mcmcOutput package to check for convergence of the main parameters and to look at their distributions. The summary table will include convergence statistics for all.

# Convert to an mcmcOutput object
library(mcmcOutput)
( mco <- mcmcOutput(out) )
# Object of class 'mcmcOutput'; approx. size 159.41 MB
# MCMC values from jagsUI object ‘out’
# Model definition file: SCRrect.jags
# Model run date: 2020-05-23 11:02:11
# MCMC chain generation took 42.83 mins
# Adaptation was sufficient.
# The output has 3 chains each with 10000 draws.
# It has 6 parameters with 664 nodes monitored:
#       nodes
# p0        1
# sigma     1
# omega     1
# N         1
# AC      440
# w       220

The print-out of mco gives some basic information about the object, but it doesn’t give any summary or diagnostic statistics. For that we need summary. We know there are 664 nodes, so we’ll use View instead of displaying everything in the Console – and we’ll only show the first part of the output from View here.

View(summary(mco))
# MCMC values from jagsUI object ‘out’
# The object has 664 nodes with 10000 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 1.00; NONE are greater than 1.10; 20 (3%) are NA.
# MCEpc is the Monte Carlo standard error as a percentage of the posterior SD:
#         largest is 3.0%; NONE are greater than 5%; 20 (3%) are NA.

An overview of the diagnostic statistics appears in the Console, and none of them gives any cause for concern. The 20 NAs correspond to the w‘s for the animals detected, which we know are 1, so Rhat and MCEpc are meaningless.

Next we’ll do the diagnostic plots and posterior plots for the four main parameters, and check that N keeps below M.

diagPlot(mco, params=1:4)
plot(mco, 1:4)
# check that N << M
max(mco$N)
# 202
M
# 220

This all looks fine, N is clearly less than M (and omega less than 1) and mean sigma is 258m, so our buffer of 1000m is ok.

Estimating density

We have a marginal posterior distribution for N, the number of activity centres in the state space, but that depends on a fairly arbitrary choice of buffer. What we really want is the density, the number of stoat activity centres per unit area. We get that by dividing N by the area of the state space:

( area <- (xmax - xmin) * (ymax - ymin) ) # km2
# 25
D <- mco$N / area
postPlot(D, xlab="Density per sq km", showCurve=TRUE)
mean(D)
# 2.657516
HDInterval::hdi(D)
# lower upper
#  1.40  4.28
# attr(,"credMass")
# [1] 0.95

A big advantage of MCMC and the “mountain of numbers” way of describing a posterior distribution is that you can do some simple (vectorized) arithmetic and get a new “mountain of numbers”: N / area produces a new MCMC chain for density, D, and we can get a posterior mean and credible interval.

D is a discrete variable but not integers: postPlot handles integers nicely, but histograms of discrete variables can be a mess, so we use a density plot with showCurve=TRUE.

The density of stoats is the result we were looking for, but we can get additional information from the JAGS output.

Home range radius

The value of sigma indicates the range of movement of animals in the population we are dealing with, which intuitively must be related to home range size. Home range size implies that the home range has a ‘hard’ boundary, while the SCR concept of activity center implies a fuzzy boundary – the probability of finding the animal declines with distance with no ‘hard’ limit. In 2 dimensions this means a bivariate normal distribution

Researchers tracking animals with telemetry are aware of the hard boundary problem, and a standard way of defining the home range is the “95% convex polygon”, the smallest convex polygon containing 95% of the relocations. If we apply this to the activity-center model with a half-normal detection function, we find that 95% of the probability mass lies within $2.45\sigma$ of the activity center:

HRR <- mco$sigma * 2.45
postPlot(HRR, xlab="home range radius, km")

Plotting the activity centers

We monitored the activity centre locations, AC, and the indicator, w, of whether an animal is real or a phantom. We will remove the AC coordinates for the phantoms, replacing them with NA. But first we use the window function to reduce the number of draws from 30,000 to 3,000 for plotting.

# reduce the number of iterations (30k is too many to plot)
mcox <- window(mco, thin=10)
# get AC coordinates, remove phantoms
ACall <- mcox$AC  #includes phantoms
w <- mcox$w
w[w==0] <- NA
ACs <- sweep(ACall, 1:2, w, "*")
str(ACs)
# num [1:3000, 1:220, 1:2] -1.71 -1.15 -1.74 -0.999 -1.48 ...
mean(is.na(ACs))  # phantoms have NA
# 0.6978864
# get posterior means of coordinates
post.means <- apply(ACs, 2:3, mean, na.rm=TRUE)
str(post.means)
# num [1:220, 1:2] -1.574 -1.543 -0.76 1.262 0.833 ...

Now we can do the plot: we’ll plot the AC locations for the 20 animals detected plus 4 animals we didn’t detect.

par(mar=rep(1,4))
plot(Detlocs, xlim=c(xmin, xmax), ylim=c(ymin, ymax),
    pch=3, col='red', axes=FALSE)
for(i in 1:24)
  points(ACs[, i, ], pch='.', col=i)
plotrix::boxed.labels(x=post.means[1:20, 1], 
      y=post.means[1:20, 2], labels=1:20)

There are lots of animals on the eastern and central parts of the study area. The plot is a bit cramped and the clouds of points for different animals overlap. Let’s do 4 plots and allocate animals to plots to minimize overlap. Each panel below has 5 animals that were detected plus 1 that was not detected.

par(mfrow=c(2,2), mar=rep(1,4))
toplot <- c(2,3,5,7,19,21,4,6,8,9,17,22,1,11,12,14,
    15,23,10,13,16,18,20,24)
colours <- rep(c(2:6,1), 4)
for(i in 1:24) {
  if(i %% 6 == 1)
    plot(Detlocs, xlim=c(xmin, xmax), ylim=c(ymin, ymax),
         pch=3, col='red', axes=FALSE)
  points(ACs[, toplot[i], ], pch='.', col=colours[i])
  if(i %% 6 !=0)
    plotrix::boxed.labels(x=post.means[toplot[i], 1],
        y=post.means[toplot[i], 2], labels=toplot[i])
}

The cloud of points represents what we know about the AC location for each animal. A small cloud means that we are more certain about the location than a large cloud: we are more certain about the AC location of #13 than #20 in the last plot. Note that the size of the cloud does not represent the size of the animal’s home range.

The clouds of points for the animals we didn’t detect cover the whole state space, though they are more dense away from the trap array; we think their ACs are probably outside the array, where the probability of detection is lower.

Abundance in a core area

The size of the state space is debated. We argued above that the state space must to large enough to include any animal with a non-negligible probability of detection; but then we are extrapolating our density estimate to areas for which we have negligible data. With a smaller state space, maybe $2\sigma$, our extrapolation is more cautious, but if we detect animals with activity centers outside this smaller state space we will over-estimate the density.

A good compromise is to fit the model with the larger state space, then estimate density for a smaller “core” area around the detector array. For our example data set, we might limit the core to an area within 250m of a detector, as shown below.

To get a posterior distribution for the number of activity centers inside the core, we check which are inside for each step of the MCMC chain, ignoring the phantom animals.

# allow a buffer of only 250m around the detector array
corebuf <- 0.25  # kilometres
cxmin <- min(Detlocs[,1]) - corebuf
cxmax <- max(Detlocs[,1]) + corebuf
cymin <- min(Detlocs[,2]) - corebuf
cymax <- max(Detlocs[,2]) + corebuf

# plot it
par(mar=rep(1,4))
plot(Detlocs, xlim=c(xmin, xmax), ylim=c(ymin, ymax),
    pch=3, col='red', axes=FALSE)
for(i in 1:24)
  points(ACs[, i, ], pch='.', col=i)
rect(cxmin, cymin, cxmax, cymax, lwd=2, border='blue')

# see which ACs are in the core, count them up
NinCore <- numeric(dim(ACs)[1])
for(i in seq_along(NinCore)) {
  draw <- ACs[i,,]
  inCore <- draw[,1] > cxmin & draw[,1] < cxmax & 
      draw[,2] > cymin & draw[,2] < cymax 
  NinCore[i] <- sum(inCore, na.rm=TRUE)
}
postPlot(NinCore, xlab="Number in core")

# convert to density
( Acore <- (cxmax - cxmin) * (cymax - cymin) ) # km2
# 12.25
Dcore <- NinCore / Acore
postPlot(Dcore, xlab="Density per sq km", showCurve=TRUE)

In this case we find that the result for the core area (2.64 activity centers per km2, HDI 1.39-4.08) is similar to the earlier result for the whole state space (2.66, 1.40-4.28).

Download a ZIP file with the code here.