A spatially explicit capture–recapture estimator for single‐catch traps

Abstract Single‐catch traps are frequently used in live‐trapping studies of small mammals. Thus far, a likelihood for single‐catch traps has proven elusive and usually the likelihood for multicatch traps is used for spatially explicit capture–recapture (SECR) analyses of such data. Previous work found the multicatch likelihood to provide a robust estimator of average density. We build on a recently developed continuous‐time model for SECR to derive a likelihood for single‐catch traps. We use this to develop an estimator based on observed capture times and compare its performance by simulation to that of the multicatch estimator for various scenarios with nonconstant density surfaces. While the multicatch estimator is found to be a surprisingly robust estimator of average density, its performance deteriorates with high trap saturation and increasing density gradients. Moreover, it is found to be a poor estimator of the height of the detection function. By contrast, the single‐catch estimators of density, distribution, and detection function parameters are found to be unbiased or nearly unbiased in all scenarios considered. This gain comes at the cost of higher variance. If there is no interest in interpreting the detection function parameters themselves, and if density is expected to be fairly constant over the survey region, then the multicatch estimator performs well with single‐catch traps. However if accurate estimation of the detection function is of interest, or if density is expected to vary substantially in space, then there is merit in using the single‐catch estimator when trap saturation is above about 60%. The estimator's performance is improved if care is taken to place traps so as to span the range of variables that affect animal distribution. As a single‐catch likelihood with unknown capture times remains intractable for now, researchers using single‐catch traps should aim to incorporate timing devices with their traps.


Introduction
Animal density is a crucial parameter in wildlife management and conservation (Buckland et al. 1993;Marques et al. 2013) and there is often interest in understanding how and why density varies in space (Gaston 2003). Spatially explicit capture-recapture (SECR) models provide a tool for investigating this as they incorporate spatial information on where captures are made (Efford et al. 2009;Royle et al. 2009;Gerber et al. 2012;Noss et al. 2012).
A variety of different detectors or traps are used in capture-recapture or SECR studies. The majority of studies of small mammals use single-catch traps that catch and hold a single animal at a time (Efford et al. 2009;Krebs et al. 2011;Gerber and Parmenter 2015). Multicatch traps also hold an individual animal until it is released but are able to simultaneously hold multiple individuals. Examples include mist nets for birds and pitfall traps for lizards. Proximity detectors are devices that record the presence of an individual without actually holding it, and unlike the previous two detector types allow an individual to be detected at more than one detector during an occasion. Camera traps, acoustic devices, and hair snares are all examples of proximity detectors.
The characteristics of the type of trap determine the specification of the detection process component of the SECR model (Efford et al. 2009). Capture in either a multicatch or single-catch trap precludes capture in any other trap during that occasion. The competition between traps for individuals leads to a competing risks formulation for multicatch traps, but single-catch traps have the additional complexity that once they are full, they are effectively unable to catch any other individuals. A suitable capture model for single-catch traps therefore needs to account for a second kind of competing risk, that of competition among individuals for traps (Efford et al. 2009). The construction of a suitable likelihood for single-catch traps is considerably more complicated than for multicatch traps, and to date no likelihood function for single-catch traps currently exists (Efford et al. 2009;Royle et al. 2013). Consequently, the multicatch trap estimator is typically used for the analysis of single-catch trap data.
Trap saturation can be calculated as the average proportion of traps that are occupied at the end of an occasion. As explained above, the extent to which the multicatch estimator assumption that traps do not fill up after catching an individual is violated depends on the degree of trap saturation. The multicatch estimator is therefore expected to perform well for low levels of trap saturation. Efford et al. (2009) conducted simulations that explored the performance of the multicatch trap estimator when applied to single-catch trap data. Three distributions for the activity centers were considered: a homogenous Poisson distribution, a Neyman-Scott distribution (with clustered centers), and an inhomogeneous Poisson distribution with an east-west linear gradient in density. The fitted model assumed constant density and a halfnormal detection function that uses two parameters (g 0 , which determines detection function height, and r, which determines its range).
They reported that in all cases, even at high levels of trap saturation (of around 86%), the multicatch estimator of both the density and r parameters performed well. There was negative bias in the g 0 parameter that increased with increasing trap saturation. The only scenario that exhibited slight bias in density (of around À5%) was that with a gradient in the density of activity centers and a high degree of trap saturation. The tentative conclusion was that the multicatch estimator may be sufficiently robust to use with single-catch traps as long as extreme trap saturation is avoided.
Traditionally, data from live-trapping studies do not contain actual capture times. However, devices that record times when a trap is triggered are available and have been used by Cowan and Forrester (2012) to study the activity patterns of possums. A continuous-time SECR model for proximity detectors that record exact capture times was developed by Borchers et al. (2014). With slight modifications, it can be used to obtain a single-catch trap likelihood. This study presents a single-catch trap likelihood for situations in which capture times are recorded and uses simulation to compare the performance of the associated likelihood-based estimator with that of the multicatch estimator under various scenarios.

Materials and Methods
We assume that the actual times of capture in singlecatch traps are available, and model the process generating detections as a competing hazards survival process (Borchers and Efford 2008) in which "death" corresponds to being caught and all individuals become "alive" again after release. Each individual is exposed to trap-specific hazards that we assume are at any time independent of the individual's capture history up to that time (although the model is easily extended to estimate different hazard levels before and after first capture as per model M b ).
The likelihood for single-catch traps needs to account for the consequences of a trap catching and holding an individual. The first consequence is that the trapped individual cannot be caught at any other trap until it is released, that is, the individual's exposure to detection by all other traps is zero for the remaining period of capture (the competing hazards formulation takes care of this). The second consequence is that the trap in which the individual is held cannot catch any other individuals until the time of release, that is, exposure to that trap for all other individuals is zero.
If we were dealing with proximity detectors, it would be straightforward to handle latent times of capture by integrating times out of the likelihood as done by Barker et al. (2014) (although their model is for abundance rather than density and does not include both time and space). However, the fact that single-catch traps induce a dependence between individuals complicates matters and means that a high-dimensional integral would need to be solved.

Notation
There are n unique individuals caught over a survey of duration T with an array of K traps. If release times are the same for all traps, then this leads to a natural definition of occasion (for discrete SECR models), and the survey duration T can be divided into L occasions.
As is typical for SECR models, it is assumed that the individuals have fixed activity center locations for the duration of the survey period: x i for the ith individual, which is a distance d k ðx i Þ from trap k. Detection probability is a decreasing function of d k ðx i Þ. The number of times the ith individual is caught at the kth detector is denoted by x ik , and instead of a capture history of length L, we have the capture times of the x ik L captures t ik ¼ ðt ik1 ; . . .; t ikx ik Þ at trap k, and t ¼ ft ik g (i = 1,. . ., n; k = 1,. . .,K) denotes the set of all detection times.
The hazard function (representing the mean capture rate per unit time) for the ith individual and the kth trap at time t is denoted as h k ðt; x i ; hÞ and can depend on both space (in terms of the distance from the trap to the activity center x i ) and time. h is an unknown vector of hazard function parameters. In the absence of other traps and other individuals, the "survivor function" for individual i at trap k over the whole survey (the probability of individual i not being caught in the trap by time T) is detection hazard over all traps at time t is h Á ðt; x i ; hÞ ¼ P K k¼1 h k ðt; x i ; hÞ, and the overall probability of detection in (0, T) over all detectors is p Á ðx In addition to h, / is the vector of parameters of the Nonhomogeneous Poisson Process (NHPP) governing animal density and D(x; /) indicates that the density at a point in space depends on both the / parameters and the spatial coordinate x. For example, if density follows an exponential east-west gradient then Dðx; /Þ ¼ exp b 0 þ b 1 Â x coordinate and / ¼ ðb 0 ; b 1 Þ.

A continuous-time likelihood for singlecatch traps
The likelihood for / and h is the joint distribution of the number of individuals captured n, and the density of the outcomes "x ik events, at times t ik1 \ . . . \ t ikx ikr ", for all i and k. With single-catch traps, the survival function term needs to take account of traps having been taken out of action by catching other individuals. Exposure to any particular trap falls to zero as soon as that trap catches any individual, and once an individual is caught in a particular trap, it cannot be caught in any other trap until it is released.
To construct a likelihood with these features, we define an indicator variable a k ðtÞ that is 1 if trap k is unoccupied at time t and zero otherwise (k = 1,. . ., K), and we define another indicator variable v i ðtÞ to be 1 if individual i is not in a trap at time t, and zero otherwise (i = 1,. . ., n). (These variables are readily calculated from the observed capture and release times of individuals at each trap.) The hazard function for individual i for trap k at time t is then conveniently written as v i ðtÞa k ðtÞ h k ðt; x i ; hÞ. The survivor function for individual i to time t is defined to be S Á ðt; The likelihood for / and h for single-catch SECR surveys then becomes: hÞ dx, and the integral is over all possible activity center locations that could have led to a detection on the survey. The term p Á ðx; hÞ is the overall probability of being caught during the survey, which depends on the combined detection hazard h Á ðt; x; hÞ over the duration of the survey. This in turn depends on a k ðtÞ (k = 1, . . ., K), which depend on random variables (the times of capture in each trap). Calculating p Á ðx; hÞ requires taking expectation over these K random variablessomething that is prohibitively computationally expensive.
Our estimator therefore involves maximizing the above likelihood equation with k(/, h) replaced bŷ which depends on the observed a k ðtÞ (k = 1, ..., K).
Consequently, the proposed estimator may not be an MLE and may not enjoy the asymptotic properties of MLEs. We evaluate the bias of the estimator and the coverage of a confidence interval estimator based on the observed information, by simulation.

Simulations
As stated in the introduction, Efford et al. (2009) found that the multicatch model estimator exhibited slight bias when there was a gradient in density, although an estimator with a constant density model was used in those cases and hence both the detection and density components of the model were misspecified. The simulations conducted here use the same form of density model in simulating and estimating and contrast the performance of the multicatch estimator with that of the single-catch estimator for other kinds of nonconstant density surfaces. We consider a range of NHPPs, with either exponential or quadratic rate parameters as a function of distance east. Table 1 and Figure 1 provide details of the scenarios used in the simulations. Note that scenario 3 in the quadratic simulations is similar to scenario 2 but with the maximum density being shifted from the center of the trap array to the right-hand side of the array area.
Except where stated otherwise, all simulations are over 5 9 24-h occasions (i.e., all trapped individuals are released simultaneously after each 24-h period) with a 5 9 4 array of traps and use a r of 100 m, trap spacings of 100 m, and a g 0 of 0.2. For all scenarios, single-catch trap data with observed capture times are simulated and two estimators (namely the discrete time SECR multicatch trap estimator and the single-catch trap estimator proposed in this study) used to estimate the parameters of interest. In both cases, the estimators use the correct form of NHPP rate parameter (exponential or quadratic). As explained in Borchers et al. (2014), the hazard function can be specified in a way that links it with the discrete time model to allow the same detection function to be fitted when the performance of the two models are compared. The model parameters are estimated using an integration area constructed with a buffer equal to 4 9 r, but the estimated density surfaces are evaluated within the convex hull of the trap array with a buffer of width 2 9 r added.
The approach used to simulate single-catch trap detection times is adapted from a method for simulating competing risks data ( (Beyersmann et al. 2009). Individuals compete for traps, and hence, the capture of one individual changes the relative hazards of capture elsewhere for all other individuals. For this reason, the simulation cannot generate capture times for each individual in isolation and needs to move forward with time rather than loop over individuals.
The steps of the simulation are summarized below: 1. A population of individuals from the given NHPP is simulated. Function sim.popn from the R package secr (Efford 2013) was used for this step. 2. The total hazard across all traps for each individual is calculated and used to generate a vector of capture times (one for each individual). We assume a constant hazard through time leading to the density of capture times following an exponential distribution. 3. The minimum capture time from this vector is taken and the rest discarded. If this time is greater than the end of the study the simulation ends, if not the time is taken to be the capture time. The time of release is also calculated and is based on the assumption that all traps are checked and reset on a set time each day (08:00 used in these simulations). 4. The particular trap where the capture event took place is then drawn from a multinomial distribution using the relative hazard at each as yet unfilled trap as the appropriate vector of probabilities, where the relative hazard for the kth trap is ðh k = P K h k Þ , and the sum is over all unfilled traps at the given capture time. 5. The total hazard from the remaining traps and the revised trap-specific relative hazards are recalculated.
A new vector of capture times is simulated and the minimum of these times added to the last capture time. If this new capture time exceeds the release time from step 3, the time is discarded and step 2 restarted from the release time, if not it becomes the next capture time and this step is repeated.
The statistical computing language R (R Core Team 2013) is used for the analysis and the R package secr (Efford 2013) used to fit the multicatch models. Computations are performed using facilities provided by the University of Cape Town's ICTS High Performance Computing team (http://hpc.uct.ac.za).

Model evaluation
The performance of the estimators is evaluated in a variety of ways. Firstly, the relative biases of the predicted mean density over the area (D) and of the detection function parameters (ĝ 0 andr) for both exponential and quadratic simulations, and of the density slope parameter (D slope ) for the exponential simulations are calculated. The estimated parameters of the quadratic coefficients are not reported as they are correlated and are more difficult to interpret than the slope parameter of the exponential rate parameter. Secondly, two measures of overall model performance that are based on predicted density at each point in space are calculated and reported, namely the root-mean-squared prediction error (RMSPE), and the Table 1. Details of four different exponential and quadratic scenarios used in the simulations. D Max is the maximum density (at 4r from the trap array for the exponential simulations), D S refers to the density at the start of the trap array (where the x coordinate is equal to zero), D is mean density, "Unique" is the mean number of unique individuals captured, and "Trap %" refers to trap saturation and is the proportion of traps occupied at the end of each occasion. Means have been rounded off. root-mean-squared bias (RMSB). These measures of model performance are calculated for two different areas: the "full" area which extends 2r beyond the trap array, and the "reduced" area which is defined as the area encompassed by the convex hull of the trap array. The RMSPE and RMSB are calculated as follows: where D m is the mean estimated density at the mth point in space (m = 1,. . .,M) averaged over the R simulations.

Results
With single-catch traps, there is an upper limit on the total number of captures over the survey, which is equal to the number of traps multiplied by the number of occa-sions. With 20 traps and 5 occasions, there are a maximum of 100 captures, and consequently, the mean number of captures for these simulations is equal to the mean percentage trap saturation and only the latter is reported. Figure 2 presents a set of plots from the exponential simulations that show the estimated density surface from each simulation overlaid on the true density surface. It is apparent that at high levels of trap saturation the multicatch estimator has a tendency to flatten out the estimated density surface. Table 2 and Figure 3 show that the multicatch model underestimates the slope parameter and that the extent of underestimation varies with trap saturation. Table 2 also shows that the relative bias in mean density is similar for the two models. Consistent with the results from the simulations performed by Efford et al. (2009), the g 0 parameter is negatively biased with the multicatch estimator and again depends on trap saturation while estimates of r are unbiased. Figure 4 shows that the single-catch model has lower bias in all cases except scenario 4, and there is not much difference between the two models in terms of RMSPE particularly when not extrapolating beyond the trap array.
For the quadratic simulations, Figure 5 presents the estimated density surface plots, Figure 6 the RMSPEs and RMSBs, and Table 3  function parameters and the mean density estimates. The results are similar to those from the exponential simulations although, in addition to g 0 being negatively biased, estimates of r from the multicatch estimator tend to be slightly positively biased in scenarios 2 and 3 where the gradient in density is steeper than the other two scenar- ios. The RMSPE for the single-catch estimator is noticeably higher than that for the multicatch estimator for scenario 1 which has the flattest quadratic bump, and for scenario 3 when the peak in density is shifted to the side of the trap array. Note that while in general the multicatch estimator again performs worse for higher levels of trap saturation, scenario 2 (and scenario 3 over the reduced area) has more bias than scenario 1 despite having lower trap saturation (81% vs. 96%).
On average, the single-catch estimator accurately estimates the true density surface although the occasional replicate overestimates the gradient in density. In some cases, a slight discrepancy between the average estimated density from the single-catch estimator and the true density can be seen at the edges of the density surface where no sampling occurs. It should be noted that the sample sizes produced by these simulations are not large. The simulations are rerun with 10 occasions, and the results in Figures 4 and 6 confirm that both the RMSPE and the slight bias reduces with larger sample sizes.

Comparing the estimators
When density is constant, the multicatch estimator is unbiased or nearly so even when the g 0 parameter is badly negatively biased (Efford et al. 2009). The multicatch estimator ignores the fact that occupied traps are out of action until they are reset. The estimator appears to compensate for this by underestimating the g 0 parameter. As stated by Efford et al. (2009), this compensator mechanism results in a surprisingly robust estimator of density although the incorrect estimation of g 0 would still have implications if used in movement or space-use models.
A nonconstant density surface can lead to high trap saturation in areas of high density but not in low-density areas. The assumption implicit in the multicatch trap model that traps continue to operate after catching an individual can therefore be badly violated in the highdensity areas leading to density being underestimated in those areas. The consequent underestimation of g 0 also gets applied to the traps in low-density areas where trap saturation may not be high, resulting in density being overestimated in low-density areas.
When density follows an east-west exponential gradient, the multicatch estimator therefore overestimates density where density is low in the west and underestimates density where it is high in the east. These two errors tend to cancel each other out and the estimator of mean density is nearly unbiased although it underestimates density in high-density areas and overestimates it in low-density areas. A slight negative bias is evident when evaluating density over the full area due to the steep exponential increase in the eastern part of the true density surface.
A similar thing happens with a quadratic bump in density whereby the multicatch estimator underestimates density around its peak and overestimates density at the edges resulting in a reasonably unbiased estimate of mean density over the full area. However, the estimator is Table 2. Simulation of bias in density and detection parameters estimated by the SECR multicatch estimator and the proposed single-catch estimator when data are from single-catch traps with 5 and 10 occasions and density follows an exponential gradient. Relative % bias is shown for each parameter followed by the standard error in parentheses. RB(D) is the relative bias in mean density over the area, F refers to the full area (with 2 9 r) and R to the area spanned by the convex hull of the trap array. In all cases, 500 replications were run and converged.

Scenario
Model RB(D slope ) RB(ĝ 0 ) RB(r) RB(D F ) RB(D R ) negatively biased for mean density if evaluated over the reduced area since then the error from underestimating density dominates the corresponding error from overestimating density at the edges. This bias is worst when the peak in density is centered on the trap array as in scenarios 1 and 2.
As expected, trap saturation affects the extent that the multicatch estimator underestimates g 0 . However, the steepness of the change in density also plays an important part and can lead to overestimation in r and to the deterioration of the robustness of the multicatch estimator. When the gradient is slight (as with scenario 4 in both types of simulations), the multicatch estimator performs well.
The single-catch estimator is approximately unbiased for the parameters of interest, and the confidence interval estimator has reasonably good coverage for the parameters of interest (see Table 4). It is clear that the singlecatch trap estimator has lower bias than the multicatch estimator for trap saturations above about 60%, and the estimators have similar RMSPEs.

Implications for trap design
If the multicatch estimator is used in a single-catch study, a trap design that lays traps out with trap density roughly proportional to expected animal density in space may avoid higher trap saturation in areas of high density.
Because the single-catch estimator sometimes estimates density with substantial positive bias when extrapolating beyond the range of explanatory variables spanned by the traps (the variable x in our simulations), it is important that traps adequately span the range of any covariate that is included in the density model. Furthermore, the variance in the single-catch estimator for density seems to increase when one extrapolates in this way. For example, the RMSPE from the single-catch estimator is worse compared to the multicatch estimator in the 1st and 3rd quadratic scenarios with 5 occasions. Both these scenarios Standard errors are calculated using the Delta method for the root-mean-square prediction error (RMSPE) and bootstrapping for the root-meansquare bias (RMSB), and error bars are plotted using two standard errors. The x-axis is ordered by trap saturation (94, 80, 64, and 60%).

Conclusion
If the focus of interest is only overall density or if density is reasonably constant, then the multicatch estimator should perform well. However, this performance deteriorates with high trap saturation and increasing density gradients. Furthermore, the multicatch estimator is poor at estimating the height (but not range) of the detection function and the detection function parameters may be of interest in their own right (for example to inform models of animal movement).
By contrast, the single-catch estimators of density, distribution, and detection function parameters are found to be unbiased or nearly unbiased in all scenarios considered. If accurate estimation of the detection function is of interest, or if density is expected to vary substantially in space, then there is merit in using the single-catch estimator.
In the absence of a single-catch trap likelihood that does not require observed capture times, we recommend that where possible researchers who are using singlecatch traps and are interested in modeling variation in density in space incorporate timing devices and use a single-catch trap estimator when trap saturation is expected to be above about 60%. Table 3. Simulation of bias in density and detection parameters estimated by the SECR multicatch estimator and the proposed single-catch estimator when data are from single-catch traps with 5 and 10 occasions and density follows a quadratic shape. Relative % bias is shown for each parameter followed by the standard error in parentheses. RB(D) is the relative bias in mean density over the area, F refers to the full area (with 29r) and R to the area spanned by the trap array. Any replications that did not converge or estimated negative or very large variances for any parameter were excluded, "Reps" reports the number of included replications.