A local evaluation of the individual state‐space to scale up Bayesian spatial capture–recapture

Abstract Spatial capture–recapture models (SCR) are used to estimate animal density and to investigate a range of problems in spatial ecology that cannot be addressed with traditional nonspatial methods. Bayesian approaches in particular offer tremendous flexibility for SCR modeling. Increasingly, SCR data are being collected over very large spatial extents making analysis computational intensive, sometimes prohibitively so. To mitigate the computational burden of large‐scale SCR models, we developed an improved formulation of the Bayesian SCR model that uses local evaluation of the individual state‐space (LESS). Based on prior knowledge about a species’ home range size, we created square evaluation windows that restrict the spatial domain in which an individual's detection probability (detector window) and activity center location (AC window) are estimated. We used simulations and empirical data analyses to assess the performance and bias of SCR with LESS. LESS produced unbiased estimates of SCR parameters when the AC window width was ≥5σ (σ: the scale parameter of the half‐normal detection function), and when the detector window extended beyond the edge of the AC window by 2σ. Importantly, LESS considerably decreased the computation time needed for fitting SCR models. In our simulations, LESS increased the computation speed of SCR models up to 57‐fold. We demonstrate the power of this new approach by mapping the density of an elusive large carnivore—the wolverine (Gulo gulo)—with an unprecedented resolution and across the species’ entire range in Norway (> 200,000 km2). Our approach helps overcome a major computational obstacle to population and landscape‐level SCR analyses. The LESS implementation in a Bayesian framework makes the customization and fitting of SCR accessible for practitioners working at scales that are relevant for conservation and management.


| INTRODUC TI ON
Spatial capture-recapture (SCR) models are now routinely used in ecological studies to estimate density of animal populations using encounter history data from recognizable individuals (Borchers & Efford, 2008;Efford, 2004;Royle & Young, 2008). More recently, extensions of SCR have been used to investigate many aspects of spatial ecology (Royle, Fuller, & Sutherland, 2018), including dispersal and survival (Ergon & Gardner, 2014;Schaub & Royle, 2014), landscape connectivity (Sutherland, Fuller, & Royle, 2015), habitat fragmentation (Bischof, Steyaert, & Kindberg, 2017), landscape conservation management (Morin, Fuller, Royle, & Sutherland, 2017), and epidemiology (Muneza et al., 2017). The basic principle of SCR models is the use of spatial patterns of individual detection/ non-detections to estimate detection probability as a function of the distance from individual activity centers. The activity center positions are themselves latent variables which are estimated using a spatial point process model.
A growing number of studies employ noninvasive methods, such as camera trapping and noninvasive genetic sampling, for monitoring wildlife. These approaches allow the collection of data suitable for SCR analysis at the level of populations and landscapes (Bischof, Brøseth, & Gimenez, 2016). The greatest hurdle to scaling up inferences to increasingly large spatial domains is the computational burden associated with fitting spatially explicit model and individual-based model to data collected over a large spatial extent. Point process models offer flexibility and are an important strength of the state-space formulation of SCR; however, they can quickly become computationally intensive or even prohibitive with increasing size of the spatial domain .
Here, we motivate the need for improved computational efficiency of SCR models using our own challenge attempting to estimate population densities of three large carnivores species (wolves Canis lupus, Brown bears Ursus arctos, and wolverine Gulo gulo) using noninvasive genetic sampling (NGS) data collected across their entire range in two countries (Norway and Sweden). For example, fitting a classical formulation of the SCR model to wolverines NGS data collected over its entire range in Norway using a Bayesian analysis (Royle & Young, 2008) was not possible on a standard desktop computer due to memory constraints. A typical SCR model assumes that the study area is large enough to contain several individuals representing the population of interest.
This implies that, during the study period, individuals occur within smaller subsets (i.e., windows) of the study area defined by the home range size. Therefore, the spatial pattern of observations of known individuals will be restricted to areas defined by their space use and ranging behavior over the study duration. This is explicitly accounted for in SCR by modeling detection probability as a decreasing function of the distance between an individual's activity center (AC) and a detector (e.g., device: camera and hair snare; observer: NGS transects). As the size of the study area (hereafter named spatial domain) increases relative to the species home range size, the fraction of the spatial domain not used by a given individual (i.e., where its detection probability is effectively 0) increases. In the current SCR formulation, all detectors are included in the likelihood optimization when fitting models, which is both inefficient and useless because null detection probability is non-informative. Therefore, removing likelihood calculations for an individual in parts of the spatial domain where its detection probability is close to 0, which we refer to as "local evaluation," should result in improved computational efficiency without affecting SCR parameter estimates.
Prior information on the study species' ranging behavior can be used to inform the spatial scale of the localized likelihood evaluation, which essentially produces individual level state-spaces that are nested within the entire state-space of SCR models. In other words, local evaluation restricts the spatial domain available to each individual. While the local evaluation approach is implemented in oSCR (Sutherland, Royle, & Linden, 2017), a package for fitting SCR models using integrated likelihood, the approach has not been implemented for more time-intensive and widely applied Bayesian analyses of SCR models. In addition, frequentist implementation of SCR does not offer the same flexibility and accessibility to build complex models that have made Bayesian SCR models popular with many ecologists.
Here, we describe an approach for local evaluation of the individual state-space (LESS) for computationally efficient fitting of large-scale SCR models in a Bayesian framework. We use simulations to assess the performance of the LESS approach, both in terms of computation speed and the performance of the estimators. We then use the LESS approach to estimate wolverine density throughout Norway using noninvasive genetic sampling data from the national monitoring program ( Figure 1) on a standard desktop computer, a task that we were unable to execute using SCR without LESS. We provide example R codes and functions for implementing our approach in JAGS (Plummer, 2003).

| Bayesian formulation of SCR models
In SCR, the spatial locations of detections and non-detections at a set of detectors are used to estimate the latent locations of individual activity centers (ACs). SCR models are hierarchical state-space F I G U R E 1 Photography of a wolverine (Gulo gulo) © Roy Andersen, Rovdata models that combine a spatial point process model describing the spatial distribution of individual ACs, and an observation model, describing the relationship between individual detection probability and distance to its AC. The classical half-normal detection model assumes that the probability p of detecting individual i at detector j decreases with distance between the detector and the AC (D ij ): where p 0 and σ are the magnitude and scale parameters, respectively.
While p 0 represents the detection probability at the location of the AC, σ is directly related to the width of the utilization distribution, that is, a function of home range size (Royle, Chandler, Sollmann, & Gardner, 2013). More generally, σ is related to the extent of space used over the period of study (see Royle et al. (2013) to relate σ to area used).
Standard SCR models assume that location of ACs (s i ) is uniformly distributed across the spatial domain (S) which is modeled as a homogeneous point process: In order to constrain the location of ACs within suitable habitat (e.g., irregularly shaped study area), we used the "ones tricks" (Chandler, 2015;Meredith, 2016;Spiegelhalter, Thomas, Best, & Lunn, 2003). This trick consists of determining whether an individual AC was located within suitable habitat or not pOK i (1/0) and rejecting the proposed AC location if pOK i = 0 using: where OK is a set of ones.
Bayesian analysis by Markov Chain Monte Carlo (MCMC) with data augmentation can be used to analyze SCR models (Royle, Dorazio, & Link, 2007;Royle, Karanth, Gopalaswamy, & Kumar, 2009). Data augmentation involves (a) augmentation of the dataset with an arbitrary (but sufficiently large) number of individuals that were never detected, and (b) use of a latent binary variable z i reflecting the realization of the inclusion probability ψ, which is an estimate of the proportion of individuals truly present in the population.
Estimated abundance (N) is then derived by summing over the inclusion vector z, and density d is derived by dividing N by the area of S. The local evaluation of the state-space was performed using an inner square evaluation window assigned to each individual to restrict the placement of their AC (hereafter the AC evaluation window) and a larger one to restrict the window within which their associated detection can occur (hereafter the detector evaluation window; Figure 3a).

| Local evaluation of the state-space
For each detected individual, we centered the AC and detector evaluation windows on the centroid of all its detections (Figure 3a (1) F I G U R E 2 (a) Illustration of the half-normal detection function (Eq. (1) representing a declining detection probability (p) of an individual with increasing distance from its activity center (AC) in spatial capture-recapture models. Values for σ and p 0 were set to 2 and 0.07, respectively. A restriction window with a width of 9σ centered on the AC covers 99.9993% of the total range of p. (b) Spatial illustration of the detection probability distribution of an individual using the half-normal detection function. Black shading illustrates the increasing detection probability of an individual around its AC (white dot). The white circle corresponds to a circular area with a radius equal to σ. The white square represents a restriction window with a width of 9σ centered on the AC Distance from activity center In contrast to the frequentist local evaluation approach implemented by Sutherland et al. (2017), Bayesian formulation of SCR involves data augmentation (Royle et al., 2007(Royle et al., , 2009. Therefore, it is necessary to define AC and detector evaluation windows for the unobserved (augmented) individuals. To do so, we juxtaposed AC evaluation windows of augmented individuals over the full extent of the study area (Figure 3c), thus generating a constant density of augmented ACs within the spatial domain. As a result, the width of the AC window relative to the size of the study area dictates the number of augmented individuals. If a larger number of augmented individuals is needed, it is possible to superimpose layers of juxtaposed squared windows ( Figure 3d). Depending on their location, AC windows of both detected and augmented individuals may include non-suitable habitat (i.e., rejection of AC location if in non-suitable habitat; Figure 3b-c). In such cases, the size of the area available to locate ACs differs among individuals. We therefore defined the inclusion parameter ψ for each individual i as a function of the proportion of suitable habitat (prop.habitat) within their AC window to obtain unbiased density estimates:

Where
The formulation of a Bayesian SCR model with LESS relies on proper indexing to specify the extent of AC and detector windows ( Figure 4). R code for generating the indexes necessary to implement a SCR model with a LESS is provided in the Supporting Information Appendix S1 and S2.

| Population and survey characteristics
We used a 50 × 50 detector array (2,500 detectors located one distance unit (du) apart from each other). The spatial domain was defined as the area including the detectors and a buffer of 4 du (2σ, see below). We simulated two different population sizes (50 and 100 individuals) and assumed demographic closure (i.e., no recruitment emigration/immigration or death). Individual AC locations were drawn as a uniform random sample from the spatial domain.
Following the recommendation by (Sun, Fuller, & Royle, 2014) for the choice of detector spacing relative to σ, we simulated binary detection data y using Eq. (6) with σ = 2, and p 0 = 0.07 leading to an average detection rate of 66% (range: 57%-78%) of alive individuals (N): Illustration of the local evaluation of the individual state-space (LESS) procedure in spatial capture-recapture (SCR). SCR with LESS involves using evaluation windows to restrict the area of the spatial domain used by each individual. The black area and small white dots represent the spatial domain and detectors, respectively. (a) Example of the placement of the activity center (AC) evaluation window (inner dark red square) and detector (outer paler red square) evaluation window centered on the centroid (triangle) of all detections of one individual. "Extension" represents the distance by which the detector window extends beyond the edge of the AC window. (b) Darker inner square evaluation windows delineate the areas within which ACs are located and lighter outer evaluation windows delineate the additional areas within which detections are considered. (c) A grid of AC and associated detector evaluation windows is established to accommodate a layer of augmented individuals. AC windows and associated detector window are shown in lighter shading for three example individuals. (d) Additional layers of AC and detector windows grids are superimposed to reach a desired level of augmentation

| Local evaluation of the state-space
We investigated the performance (see below) of SCR with LESS using different AC and detector evaluation windows widths. We used 3σ (6 du) and 5σ (10 du) as the width of the AC window ( Figure 5). We then sequentially increased the width of the detector evaluation window, so that the detector window extended beyond the edge of the AC window by 1, 2, 3, and 4 times σ. This resulted in detector windows of different width according to the AC window chosen ( Figure 5). We adjusted the number of layers of AC windows (corresponding to augmented individuals) to obtain >4× number of simulated individuals.

| SCR key parameters
For reference, we fitted an SCR model without local evaluation.
We then ran 100 simulations of each scenario and calculated the relative bias � and the precision of N , ̂ and p 0 for each simulation scenario using the coefficient of variation CV = SD(̂) × 100 (Walther & Moore, 2005), where n is the number of iterations, SD is the standard deviation, θ is the true parameter value, and ̂ the estimate of the parameter obtained from MCMC.
In addition, we calculated the 95% credible interval coverage as the percentage of simulations where the 95% credible interval contained the true parameter value. We also recorded the computing time for each model fit (i.e., time necessary to run three MCMC chains in series on a single computer core). Because we used cores with different characteristics to run the different SCR models, computing times reported only serve as a crude estimate of the time used for model fitting.

| Density predictions
One of the main advantage of SCR methods is to derive spatially explicit estimates of density. We therefore quantified the deviation between the simulated and realized density maps produced using SCR at the habitat grid cell level (1 × 1 du). To derive the predicted average space used, we used the uncertainty around the estimates in the posterior samples (location of ACs and σ) obtained from the Bayesian Markov chain Monte Carlo (MCMC) and summed predicted average space used of all individuals that were considered as part the population (z = 1). To compare the deviation between the shape of true and predicted density maps, we scaled both maps to sum to 1 and calculated mean relative error (MRE) in density as: where D j is the predicted and D j the simulated density at cell j in a habitat raster consisting of h cells. We excluded the buffer area when calculating MRE. We then summarized the average MRE across all simulations for any given scenario.

| Application to the wolverine data
We fitted an SCR model with LESS to NGS data from the national monitoring program of wolverines in Norway (see description in We used AC evaluation windows of 5σ (the largest LESS restriction used in our simulation, σ≈ 6 km for females wolverines ) and detector evaluation windows that extended by 2σ beyond the edges of the AC windows. The goal of the empirical analysis was to put the LESS approach to the test on a large-scale SCR problem, and we used a simple SCR model for this task (as described in section 2.1). Comprehensive and reliable density estimation of wolverine density would require a model that, among other things, accounts for heterogeneity in density and detectability, which is beyond the scope of the analysis presented here.
To avoid misinterpretation of density and abundance estimates from our simplified analysis as actionable results with respect to population management, we do not provide absolute estimates, only report computation time and a map of relative density, to be interpreted with caution.

| Model fitting
We fitted Bayesian SCR models using Markov chain Monte Carlo (MCMC) with JAGS (Plummer, 2003) and rjags (Plummer, 2016) in R version 3.3.3 (R Core Team, 2017). After an adaptive phase of 1,000 iterations, we ran 3,000 iterations of three chains thinned by three. We considered models as converged when the Gelman-Rubin Illustration of the different AC (in rows) and detector (in columns) evaluation windows widths used in the simulations to test the local evaluation of the individual state-space (LESS) in SCR model. The darker and lighter black areas represent the detector and buffer areas of the spatial domain. The darker and lighter red square represents the AC and detector evaluation windows, respectively. All values are expressed relative to σ. See Figure 3 for a description of the windows AC Window width (σ) statistics (gelman.diag function in coda package (Plummer, Best, Cowles, & Vines, 2006)) (Brooks & Gelman, 1998;Gelman & Rubin, 1992) was ≤1.1 and after visually inspecting trace plots for all monitored parameters. R and JAGS code for the different SCR models and simulations used are provided in Supporting Information Appendix S1 and S2, and list of priors used in Supporting Information Appendix S3.

| Simulations
All 1800 models (with and without LESS) reached convergence after 3,000 iterations. Relative bias of SCR models with LESS generally decreased with increasing width of the AC evaluation window ( Figure 6, Table 1). When using a detector evaluation window with the smallest extension beyond the AC evaluation window (1σ), N and ̂ tended to be slightly overestimated and p 0 underestimated, regardless of the width of the AC window ( Figure 6, Table 1). However, as the extension increased (>1σ), ̂ tended to be underestimated and p 0 overestimated. Regardless of the scenario (N = 50 or 100), relative bias in N (<0.005), ̂ (≈−0.01), and p 0 (<0.1) were lowest for an AC window that was equal to 5σ (10 du) and when detector evaluation window extended 2σ beyond the edges of AC window (width of the detector window = 18 du). MRE between true and predicted density surfaces was similar for SCR models with and without LESS (Table 1, Figure 7). With our survey and population characteristics, fitting SCR model with LESS (AC region = 5σ, extension between the detector and AC regions = 2σ) was between 37-57 times faster than fitting SCR models without LESS (Table 1).

| Wolverines
Using the LESS approach, we were able to fit a SCR model to wolverine NGS data over the entire range of the species in Norway in less than 22 hr. The resulting map of relative density is shown in Figure 8 (Supporting Information Appendix S4). We reemphasize that this analysis is a proof of concept serving the sole purpose of demonstrating the feasibility of using LESS to enable analysis of very large state-space problems using Bayesian methods.

| D ISCUSS I ON
Spatial capture-recapture methods are now widely used to investigate a range of important spatial ecological processes, and to inform conservation and management actions .
The emergence of efficient inventory techniques such as noninvasive genetic sampling and remote camera trapping has resulted in a substantial increase in the spatial scales of some monitoring projects. The challenge that has emerged is computational in nature: finding tractable modeling procedures that allow SCR methods to be applied over large spatial domains and make inferences about population-level processes. Our local evaluation method provides a solution that produces unbiased inferences, while vastly reducing the computational burden of enormous spatial models, and in doing so allows large-scale Bayesian SCR models to be applied using a standard desktop in a fraction of the time required in the absence of LESS.
When using LESS in a Bayesian framework with data augmentation, our results suggest that there was virtually no bias in focal parameters when using an AC evaluation window width of 5σ and a detector evaluation window that extended by 2σ beyond the edge of the AC window (detector window width = 9σ). A detector window of 9σ covers most of the half-normal distribution defining the detection function ( Figure 2). A 2σ extension of the detector window beyond the edges of F I G U R E 6 Relative bias of spatial capture-recapture (SCR) model parameters with and without local evaluation of the individual state-space (LESS) for a simulated population size of 100 individuals. Relative bias plotted against the different widths of the AC evaluation windows and different detector extensions (i.e. the distance by which the detector window extends beyond the edges of the AC window, see Figure 5). AC window width and detector window extension are expressed relative to σ. Violins show kernel density (dots: mean) of the relative bias computed from the posterior distribution of three key SCR parameters (N, ̂, and p 0 ) based on 100 simulations the AC window corresponds to the general recommendations of habitat buffer size in SCR models (Royle et al., 2013). Using those settings in LESS, we measured up to 57 times gains in computing speed.
In SCR models, the computational burden increases rapidly with the number of detectors used . The LESS approach, by using detector evaluation windows, reduces the number of detectors over which detection probability has to be estimated for each individual. In our simulated example, individual detection probability had to be estimated for each of the 2,500 detectors when a SCR model was fitted without LESS. However, using the LESS approach and our recommended AC and detector evaluation windows widths (AC: 5σ, detector: 9σ), the maximum number of detectors for which individual detection probability had to be estimated was 324. Although the location of individual activity centers Table 1 Average relative bias, coefficient of variation, coverage of key parameters (N, ̂, p 0 ), and average running time of SCR models with and without local evaluation of the state-space (LESS) from the simulation study Notes. Results are presented for different AC evaluation window width and different detector evaluation window extensions (see Figure 5), expressed relative to σ. *MRE are expressed ×10 −4 .
is also constrained to a smaller area using the AC evaluation window, such a reduction in the number of detectors used for each individual likely explains the major speed gains with LESS. This means that increasing the size of the problem (e.g., higher number of detector, larger spatial extent, or greater population size) will amplify the gains even further as demonstrated in the wolverine example. While we were unable to fit the wolverine SCR model for the entire Norwegian range of the species (>200,000 km 2 ), the LESS approach allowed us Implementing the LESS in a Bayesian SCR framework is straightforward, with proper indexing to define the windows requiring the most care (Figure 4). Although the half-normal detection function represents a circular utilization of the space by individuals (Figure 2), we chose a square evaluation window for ease of indexing (Figure 3). This simplification did not seem to introduce additional bias.
Local evaluation is already possible in a maximum likelihood framework (oSCR, ) but remains rarely used or described (but see Sun, Fuller, Hare, & Hurst, 2017) and has yet to gain the attention it deserves. Conceptually, the approach used in oSCR is identical to our "2-window" approach, but instead of using square evaluation windows, oSCR uses circular evaluation windows.
First, a "trim" value defines the width of a buffer around individual detections that identifies where ACs can be located. Then, an additional buffer (double the width of the first buffer) is used to identify the detection ranges of all possible ACs for a given individual.
The SCR model with LESS that we used to demonstrate local evaluation did not include any covariates to account for individual/ spatial heterogeneity in detection and space use. If detector-specific covariates were to be included, similar indexing as on the detection function should be employed. Ultimately, the width of the windows should be specific to the study system; an optimal trade-off between computation speed and reliability of prediction can be determined using sensitivity analysis. Customizable R code for simulation and model fitting is provided in the Supplementary information Appendix S1.
The use of data augmentation to estimate abundance in SCR is specific to Bayesian inferences (Royle et al., 2007(Royle et al., , 2009. Although proven to be a robust tool to estimate abundance, dealing with data augmentation can be a challenge with more complex model formulation. For example, the correction that we applied to the inclusion parameter (ψ i ) as a function of the size of available habitat for each individual is specific to augmented individuals. In cases where individuals show a high degree of spatial aggregation within the study area, density of augmented individuals should be larger than the maximum observed density, to ensure that enough AC windows are available for augmented individuals to estimate density in areas with high individual density. Additionally, in the case of open population SCR models , steps would have to be taken to deal with individual movement. Under an open population SCR model where a dispersal kernel is explicitly defined, or for nonstationary activity centers (Royle, Fuller, & Sutherland, 2016), there exists a spatial parameter similar to σ (dispersal σ) to redistribute potential ACs. The general rule of the LESS approach that restricts calculation where detection probability is >0 would also apply to where the relocation probability of an individual's AC is >0 (i.e., based on dispersal kernel function). As it stands, our approach does not accommodate for moving individuals because we defined static individual evaluation windows.
Ecological studies have an increasing need for methodological developments to analyze the ever-increasing size of spatially explicit datasets and produce landscape, regional, and even global inferences about population state variables (Sutherland, Brambilla, Pedrini, & Tenan, 2016;Tenan, Brambilla, Pedrini, & Sutherland, 2017). The LESS approach presented here is a pragmatic redundancy reduction approach, as it avoids calculation that are unnecessary (i.e., where detection probability is close to 0). Although some studies have used such reduction to increase computational efficiency of spatial models (e.g., Latimer, Banerjee, Sang, Mosher, andSilander (2009), Gramacy (2016)), our approach could stimulate the development of more efficient model formulation that have a built-in correlation structure.

| CON CLUS ION
SCR is a powerful ), yet computationally demanding tool . Incorporating LESS in SCR model opens the door to SCR application at large spatial scales, a prerequisite for comprehensive population-level conservation and management plans .
Because the computational burden associated with the spatial component of SCR models can be considerably reduced using LESS, it will also be possible to develop computationally efficient models that include relevant ecological processes (e.g., survival, recruitment, movement; Chandler, Hepinstall-Cymerman, Merker, Abernathy-Conners, & Cooper, 2018) at large spatial scales. Generally, we encourage practitioners to use local evaluation when performing SCR at large spatial scales, whether using Bayesian or maximum likelihood inference .

ACK N OWLED G M ENT
This work was funded by the Norwegian Environment Agency (Miljødirektoratet) and the Swedish Environmental Protection Agency (Naturvårdsverket). We thank the field staff from the State Nature Inspectorate and members of the public that collected the agency's large carnivore database Rovbase3.0 (rovbase.no).

CO N FLI C T O F I NTE R E S T
None Declared.