SCR with multi-catch and single-catch traps

So far we have only dealt with detectors which do not physically catch and detain the animal. Animals can be detected at several locations on the same occasion and detectors can record multiple animals on the same occasion. On this page we deal mainly with multi-catch traps, which do detain the animal – it can only be recorded at one detector on one occasion – but can catch several animals on the same occasion. A section at the end will talk about single-catch traps.

Examples of multi-catch traps are mist nets for birds, harp traps for bats, or pitfall traps for insects and herps. The animal is not just detected, it is captured and detained until it is released at the end of the sampling occasion. Each animal can only be caught in one trap. So animals have to “chose” their trap or, looked at from the other side, traps must “compete” for animals. This has a big effect on the model we use.

The example data

Our example data come from a simulation we do in our Boot Camps. Each participant has a counter representing a gecko. They start from a randomly-selected square on the board, their activity centre, and move N, S, E or W according to the result of a spinner. If they pass or land in a square with a trap (the yellow squares in the photo), they are captured; they record their ID, the trap ID and the occasion number. The game continues for 25 occasions. We know the true number of geckos, ie, the number of participants, even though some escape capture.

The data are in two files, one with the capture details and one with the trap locations. You can download and check the data with the following code:

capthist <- read.csv("http://bcss.org.my/data/geckoCapt.csv",
    comment="#")
str(capthist)
# 'data.frame':   44 obs. of  3 variables:
#  $ animalid: int  1 1 2 3 3 3 3 3 3 3 ...
#  $ occasion: int  5 11 17 1 3 4 9 13 15 24 ...
#  $ trapid  : int  18 11 10 4 4 4 4 5 4 4 ...

traps <- read.csv("http://bcss.org.my/data/geckoTraps.csv",
    comment="#")
str(traps)
# 'data.frame':   36 obs. of  2 variables:
#  $ x: int  45 95 125 175 205 255 55 85 135 165 ...
#  $ y: int  45 55 45 55 45 55 95 85 95 85 ...
traps <- as.matrix(traps)

The multinomial distribution

We’ve used the binomial distribution a lot to handle binary events with 2 possible outcomes; now we need to extend the idea to events with more than 2 possibilities and discuss the multinomial distribution.

For the “sock in the box” data we used for the logistic regression page, we had 105 people trying to throw a sock into a box from different distances. Throwing from 4m with their good hand, 37 succeeded in getting the sock in the box. We use the binomial distribution to get the likelihood, ie. the probability of 37 successes in 105 trials with a given probability of success, p: dbinom(37, 105, p), and the maximum likelihood corresponds to p = 0.352.

Now suppose we added 2 new boxes, placing them left and right of the target box. Many of the 70 throws which missed the target would land in one of those, but some would still miss. The data might look like this:

leftcentrerightmissed
18372228

We use the multinomial distribution to get the likelihood. Instead of using successes and trials as the input, we use a vector with all 4 values. And instead of a single value for p, we use a vector of 4 values, which must add up to 1, as we have included all possible outcomes. The probability of getting the results in the table, y <- c(18, 37, 22, 28) with probabilities p <- c(0.2, 0.3, 0.2, 0.3) is found with dmultinom(y, prob=p).

For each gecko, we need a vector with the number of times the animal was caught in each trap, one value for each of the 36 traps, plus the number of times it was not caught. And for probability of detection we also need a vector with 37 values adding up to 1.

The multinomial logit link

We will need to link the probability of capture in each trap to the distance between the trap and the AC. For the sock-in-the-box analysis we used a logit link with a linear predictor which included an intercept and terms for distance, hand, gender, etc.

For the geckos, our linear predictor will include an intercept and a term for the squared distance between the AC and the trap, so for animal i and trap j we have:$$\mu_{i,j} = \alpha_0 – \alpha_1 d_{i,j}^2$$ where $\alpha_1 = 1/(2 \sigma^2)$.

The multinomial logit, or “mlogit”, function looks like this: $$p_{i,j} = \frac {e^{\mu_{i,j}}} {1 + \sum_j{e^{\mu_{i,j}}}}$$

Because the $p_{i,j}$ must add to 1, we can’t have a linear predictor for all of them. That works well for our model: we have a linear predictor for each of the traps, then the last value in the probability vector, $p_{i,37}$, the probability of not being captured, is calculated by subtracting the sum of the capture probabilities from 1: $$p_{i,37} = 1 – \sum_{j=1}^{36}{p_{i,j}}$$

The model

The JAGS model code is shown below. It is very similar to the code used for the stoats example, with differences highlighted in purple.

# Model for multi-catch traps and rectangular state space
# File name "SCR_multinom.jags"

model{
  for(i in 1:M){
    w[i] ~ dbern(omega)
    AC[i,1] ~ dunif(xlim[1], xlim[2])
    AC[i,2] ~ dunif(ylim[1], ylim[2])
    for(j in 1:ntraps){
      d2[i,j] <- (AC[i,1] - traps[j,1])^2 + 
          (AC[i,2] - traps[j,2])^2
      expmu[i,j] <- exp(alpha0 - alpha1*d2[i,j]) * w[i]
      p[i,j] <- expmu[i,j] / (1+sum(expmu[i, ]))
    }
    p[i, ntraps+1] <- 1 - sum(p[i,1:ntraps])
      # last p = not captured
    Y[i, ] ~ dmulti(p[i, ], nocc)
  }
  
  # Priors
  omega ~ dbeta(1, 1)
  alpha0 ~ dunif(-5, 5)
  sigma ~ dunif(0, 100)
  alpha1 <- 1/(2*sigma^2)

  # Derived quantities
  N <- sum(w)
}

The code is the same up to the calculation of expmu[i,j], the exponentiated linear predictor for capture of animal i in trap j. We include w[i] here to ensure that all the capture probabilities are zero for phantoms. Probability of capture for each of the 36 traps, p[i,j], is calculated using the mlogit formula, then we calculate the last value, p[i,37], the probability of not being captured, by subtraction.

Finally we deal with the data. On the left hand side of the tilde we have a vector with the number of captures in each of the 36 traps and the number of occasions the animal was not captured. The dmulti distribution needs the vector of 37 probabilities, and the number of survey occasions.

Running the model

The state space will be a rectangle, and we’ll allow a 100m buffer around the traps.

buffer <- 100
xlim <- c(min(traps[,1]) - buffer, max(traps[,1]) + buffer)
ylim <- c(min(traps[,2]) - buffer, max(traps[,2]) + buffer)

We build the Y matrix. This is an animals x traps matrix with the number of times animal i was caught in trap j, plus a last column with the number of occasions animal i was not captured.

( nind <- max(raw$animalid) )
( ntraps <- nrow(traps) )
nocc <- 25

# create and fill the matrix
Y <- matrix(0, nrow=nind, ncol=ntraps)
for(i in 1:nrow(caphist))
  Y[caphist[i,1], caphist[i,3]] <- 
      Y[caphist[i,1], caphist[i,3]] + 1
sum(Y) == nrow(caphist)
# TRUE

# add last column
noncap <- nocc - rowSums(Y)
Y <- cbind(Y, noncap)
sum(Y) == nind * nocc
# TRUE

Now we can do the data augmentation, and lots of extra undetected animals seem to be necessary. The uncaptured animals will have rows with 25 in the last column and everything else zero.

# data augmentation
nx <- 100
Y0 <- matrix(0, nx, ntraps+1)
Y0[, ntraps+1] <- nocc
Yaug <- rbind(Y, Y0)

Now we bundle the data for JAGS, run the analysis (takes about 7 mins), and check the results:

# data augmentation
nx <- 100
Y0 <- matrix(0, nx, ntraps+1)
Y0[, ntraps+1] <- nocc
Yaug <- rbind(Y, Y0)

bdata <- list(Y = Yaug, M = nind+nx, traps = traps, ntraps=ntraps,
    nocc=nocc, xlim=xlim, ylim=ylim, 
    w = c(rep(1, nind), rep(NA, nx)))
str(bdata)
# List of 8
#  $ Y     : num [1:117, 1:37] 0 0 0 0 0 0 0 0 0 0 ...
#  $ M     : num 117
#  $ traps : int [1:36, 1:2] 45 95 125 175 205 255 55 85 135 165 ...
#  $ ntraps: int 36
#  $ nocc  : num 25
#  $ xlim  : num [1:2] -55 355
#  $ ylim  : num [1:2] -55 355
#  $ w     : num [1:117] 1 1 1 1 1 1 1 1 1 1 ...

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

library(jagsUI)
out <- jags(bdata, NULL, wanted, "SCR_multinom.jags", DIC=FALSE,
    n.chains=3, n.iter=1e4, parallel=TRUE)

library(mcmcOutput)
( mco <- mcmcOutput(out) )
# Object of class 'mcmcOutput'; approx. size 85.47 MB
# MCMC values from jagsUI object ‘out’
# Model definition file: SCR_multinom.jags
# Model run date: 2020-07-01 13:39:00
# MCMC chain generation took 7.08 mins
# Adaptation was sufficient.
# The output has 3 chains each with 10000 draws.
# It has 7 parameters with 355 nodes monitored:
#          nodes
# omega        1
# sigma        1
# p0           1
# N            1
# AC         234
# w          117

View(summary(mco))
# MCMC values from jagsUI object ‘out’ 
# The object has 355 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.01; NONE are greater than 1.10; 17 (5%) are NA.
# MCEpc is the Monte Carlo SE as a percentage of the posterior SD:
#    largest is 3.8%; NONE are greater than 5%; 17 (5%) are NA.

# check that aumentation is adequate
max(mco$N)
# 102
nind+nx
# 117

# do the plots
plot(mco)
diagPlot(mco)

We only show results for the main parameters, but convergence and precision are ok, as shown by the Rhat and MCEpc values. Data augmentation with 100 extra all-zero capture histories is adequate. The estimate of $\sigma$ is 19m, so our 100m buffer is more than $4\sigma$ and is adequate. The diagnostic plots show that the posterior distributions for $\sigma$ and $\alpha_0$ are not constrained by our arbitrary uniform priors.

Estimating density

We get the density by dividing the number of activity centres in the state space by its area. Working with the MCMC chain for $N$, we get a new MCMC chain for density, and can plot the posterior distribution and get posterior mean and credible interval.


( area <- diff(xlim) * diff(ylim)/10000 ) # ha
# 16.81
D <- mco$N / area
postPlot(D, xlab="Density per ha", showCurve=TRUE)
mean(D)
# 2.798437
HDInterval::hdi(D)
#    lower    upper 
# 1.725164 3.866746 
# attr(,"credMass")
# [1] 0.95

For these simulated data we know the true density of “geckos”: There were 24 participants with activity centres randomly placed within an square 300m x 300m = 9ha, giving a true density of 24 / 9 = 2.67 ACs/ha. The posterior mean is thus a very good estimate, especially considering that we only caught 17 of the 24 “geckos” present.

Estimating home range radius

If our simple bivariate normal model of gecko movement is plausible, we can calculate the home range radius as $2.45\sigma$. Again we can get an MCMC chain and a posterior distribution for this:

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

Plotting ACs

With Bayesian SCR analysis we can get posterior distributions for the AC locations for the animals captured. And for those not captured, though that’s not so useful. Here’s the code to prepare and plot the results:

# 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:117, 1:2] 261 251 238 218 248 ...
mean(is.na(ACs))  # phantoms have NA
# 0.5975385

# We use posterior means to place labels on the plot
post.means <- apply(ACs, 2:3, mean, na.rm=TRUE)
str(post.means)
# num [1:117, 1:2] 240 166.1 178.8 63.2 83.3 ...

# One big plot
op <- par(mar=rep(1,4))
plot(traps, xlim=bdata$xlim, ylim=bdata$ylim, pch=3, col='red',
        axes=FALSE)
for(i in 1:17)
  points(ACs[, i, ], pch='.', col=i)
plotrix::boxed.labels(x=post.means[1:17, 1], 
    y=post.means[1:17, 2], labels=1:17)
par(op)

That’s a bit messy, with overlapping clouds of points. We’ll divide them into 4 plots:

# 4 smaller plots
op <- par(mfrow=c(2,2), mar=rep(1,4))
toplot <- c(2,9, 14, 1, 12, 7, 11, 16, 3, NA, 15, 10, 5, 17,
    NA, 6, 13, 8, 4, 18, 19)
colours <- rep(c(2:5,1), 4)
for(i in 1:20) {
  if(i %% 5 == 1)
    plot(traps, xlim=bdata$xlim, ylim=bdata$ylim, pch=3,
        col='red', axes=FALSE)
  points(ACs[, toplot[i], ], pch='.', col=colours[i])
  if(i !=20)
    plotrix::boxed.labels(x=post.means[toplot[i], 1],
        y=post.means[toplot[i], 2], labels=toplot[i])
}
par(op)

Remember that the size of the cloud of points represents the uncertainty about the AC location, not the size of the home range. The last plot includes one animal that was not caught; we don’t have much idea of where its AC could be, but probably not in the middle of the trap array.

Single-catch traps

Most small-mammal studies use single catch traps which close when one animal enters and thus can only capture one animal per trapping occasion. Multi-catch traps are limited to one trap per animal, but now we have an additional constraint, one animal per trap.

No one has yet devised an algorithm to cope with both limitations at the same time1but see Distiller & Borchers (2015, Ecol.Evol. 5, 5075-5087) for a method that works when time of capture is recorded.

Fortunately, the multi-catch method works pretty well for single-catch traps, provided trap saturation (the proportion of traps occupied) is not too great. A simple work-around is to place multiple single-catch traps at each point on the trap array, and then to treat the cluster of traps as a multi-catch trap.

Efford et al (2009)2Efford, M.G., Borchers, D.L., & Byrom, A.E. (2009). Density estimation by spatially explicit capture-recapture: likelihood-based methods. In Modeling demographic processes in marked populations (eds D.L. Thomson, E.G. Cooch & M.J. Conroy), pp. 255-269. Springer investigated the use of the multi-catch algorithm for simulated single-catch data, and found no bias in density estimates up to trap saturations as high as 85%. Estimates of the base-line capture probability, $p_0$, were biased low, and they concluded that moderate levels of trap saturation merely led to lower probability of capture.

Download a ZIP file with the code here.

Leave a Reply

Your email address will not be published. Required fields are marked *