# 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: 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

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
# 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.

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: