Introduction

Imperfect detection is one of the primary challenges to the estimation of the size of wild populations. Regardless of the data collection method employed, rarely, if ever, are all individuals in a population detected. This challenge is amplified with the increasingly widespread application of non-invasive sampling methods for making landscape-level assessments across time and space, such as camera trapping and non-invasive DNA sampling (Burton et al. 2015; Beng and Corlett 2020), which often trade off local sampling intensity for extent of spatial coverage (grain vs. extent; Chandler and Hepinstall-Cymerman 2016; Steenweg et al. 2018). Capture-recapture and, more recently, spatial capture-recapture (SCR) models estimate and account for imperfect detection, thereby producing robust estimates of the focal ecological parameters (Chao 2001; Efford 2004; Lukacs and Burnham 2005; Borchers and Efford 2008; Royle et al. 2014). SCR has become particularly popular as it exploits the information contained in the spatial configuration of detections and non-detections across the study area to produce spatially explicit estimates of abundance (i.e., density; Borchers and Efford 2008; Royle and Young 2008; Royle et al. 2018).

In most field studies, detection probability—the probability of detecting a species or individual when it is present—is not only imperfect but also variable. Detection probability can vary across individuals, time, and space (Gimenez et al. 2008; Kellner and Swihart 2014; Conn et al. 2017; Guélat and Kéry 2018). The implications and treatment of individual and temporal variability in detection probability have been extensively documented in the non-spatial capture-recapture literature (Chao 2001; Link 2003; Borchers et al. 2006; Gimenez et al. 2018a). SCR models are particularly well-suited to account for spatially variable detectability, as studies are usually configured into discrete detection locations referred to as detectors (or traps) that are distributed across the study area (Efford et al. 2013; Royle et al. 2014). Spatial variation in detection probability resulting from individual space use relative to detector locations, i.e. the declining probability of detection with increasing distance from an individual’s activity center (AC), is in fact exploited by SCR models to estimate the distribution of individual ACs (Borchers and Efford 2008; Royle and Young 2008). Other potential sources of spatial heterogeneity in detectability include those caused by the study itself, such as variable search effort, and local intrinsic and extrinsic factors that influence the detectability of animals. Known sources of variation in detection probability can be modeled and accounted for in SCR (Efford et al. 2013; Efford and Mowat 2014; Royle et al. 2014). This is the case when variable sampling effort is recorded, for example, during camera trapping or non-invasive DNA sampling (Royle et al. 2009; Efford et al. 2013), or when spatial covariates, such as habitat proxies for vulnerability to detection, are used (Bischof et al. 2017; Kendall et al. 2019).

In many wildlife monitoring studies, spatial heterogeneity in detection probability remains partially unknown. Unaccounted environmental factors may impact exposure to detectors; for example, site-specific characteristics may affect visibility, or local climate can influence genotyping success rate of non-invasively collected DNA samples (Efford et al. 2013; Kendall et al. 2019). Survey effort may also vary across the study area unbeknownst to the investigator. For example, many large-scale monitoring programs combine structured sampling with unstructured data collection methods to increase the spatial extent and sampling intensity, and involve members of the public in the process (Thompson et al. 2012; Conn et al. 2017; Altwegg and Nichols 2019; Bischof et al. 2020a; Sicacha-Parada et al. 2021). Data from unstructured sources introduce unknown spatial heterogeneity in detection probability in SCR studies. In the analysis of monitoring data, ignoring the variability in detection can seriously degrade population inferences (Nichols and Williams 2006; Gimenez et al. 2008; Gerber and Parmenter 2015). However, it has also been previously shown that spatially random variation in detection probability does not seem to be a major source of bias in parameter estimates from SCR analysis (Bischof et al. 2017). In such situations, most individuals in the sampling grid will be detected by at least one detector, if detectors are spaced closely enough (high grain) and cover a large enough area (large extent) relative to home range sizes in the study population (Sollman et al. 2012; Efford and Fewster 2013).

Spatial autocorrelation in detection probability—when detectability is more similar among neighboring than distant detectors—is common in ecological studies (Guélat and Kéry 2018). Observed and unobserved spatially autocorrelated variation in detection probability could have many causes (Gaspard et al. 2019), divided into two broad categories. On the one hand is the nature of the data collection; for example, regional differences in the mobilization of volunteers for non-invasive DNA collections, variation in camera trap efficiency due to inadvertent scent contamination at a cluster of sites, or reduced physical capture success in traps installed by a less-experienced operator in their designated area (Kristensen and Kovach 2018; Bischof et al. 2020a; Tourani et al. 2020a). On the other hand are characteristics of the study species and its environment, such as spatial variation in site attractiveness or individual behaviors (e.g., shyness) that lead to variable detectability (Efford et al. 2013; Howe et al. 2013; Stevenson et al. 2021). Variable detection probability that is highly autocorrelated may disproportionately affect the overall detectability of certain individuals in the population based on the location of their ACs. When areas with virtually zero probability of detection are known, for example in clustered sampling designs, they can be specified as detector-free regions (or holes) in the sampling grid and be treated like the unsampled habitat buffer in SCR (Efford and Fewster 2013; Royle et al. 2014). In extreme cases, however, large swaths of the study area, and the individuals inhabiting them, may be left unsampled, unbeknownst to the investigator. This can occur when, for example, SCR data are obtained opportunistically (e.g., contributed by citizen scientists), where the incomplete information about the detection process makes it unclear whether areas without detections were sampled or not (i.e., true vs. false negatives; Thompson et al. 2012; Bird et al. 2014; Bischof et al. 2020a). We expect flawed inferences from SCR studies if autocorrelated heterogeneity in detection probability remains cryptic, and thus unaccounted for. Specifically, when detection probability is not correlated with animal density (Clark 2019; Paterson et al. 2019), we predict that an SCR model ignoring contiguous areas with low or zero detection probability would produce positively biased estimates of average detection probability as it is inferred from areas with detections, and therefore negatively biased estimates of density (Royle et al. 2013; Guélat and Kéry 2018).

Most, if not all, SCR analyses of empirical data that are collected across large spatial extents, ignore some sources of spatial variation in detection probability. The extent to which the misspecification of the detection process in the presence of variable and spatially autocorrelated detection probability may affect the estimates of focal parameters in SCR studies has not yet been systematically investigated. Using simulations with an envelope that includes scenarios encountered during wildlife monitoring, we quantify how unmodeled, spatially variable detection probability influences parameter estimates obtained via SCR. Specifically, our objective in this study is to identify scenarios where unmodeled heterogeneity in detection probability can lead to problematic inferences in SCR analysis. We do so with a focus on both the magnitude of variability in detection probability and the autocorrelation therein.

Methods

Spatial capture-recapture model

For the purpose of this study, we used a standard, single-session SCR model in a Bayesian framework (Royle et al. 2014; see the model code in Online Appendix 1). Our model is composed of two hierarchical levels:

  1. 1.

    The ecological sub-model reflects the underlying ecological process of interest and describes the distribution of individuals in space, i.e. density. Following a homogeneous point process (Royle et al. 2014), we assume every individual i in the population has a fixed AC (or home range center; \(s_{i}\)) and that these individual ACs are randomly distributed across the habitat S according to a uniform distribution:

    $$\begin{aligned} s_{i} \sim Uniform(S) \end{aligned}$$
    (1)

    We used a data augmentation approach to account for those individuals in the population that are not detected (Royle et al. 2007). Detected and augmented individuals make up the super-population of size M (see below). A latent state variable \(z_{i}\) describes inclusion of individual i in the population, governed by the inclusion probability \(\psi\):

    $$\begin{aligned} z_{i} \sim Bernoulli(\psi ) \end{aligned}$$
    (2)

    where \(z_{i}\) takes value 1 if individual i is a member of the true population and 0 otherwise. Population size N is therefore:

    $$\begin{aligned} N = \sum ^{M}_{i=1}{z_{i}} \end{aligned}$$
    (3)
  2. 2.

    The observation sub-model describes the individual and detector-specific probability of detection \(p_{ij}\) as a function of Euclidean distance \(d_{ij}\) between the location of an individual AC and a given detector j. We used the half-normal detection function to model \(p_{ij}\) (Borchers and Efford 2008; Royle et al. 2014):

    $$p_{{ij}} = p_{0} {\mkern 1mu} \exp \left( {{{ - d_{{ij}}^{2} } \mathord{\left/ {\vphantom {{ - d_{{ij}}^{2} } {2\sigma ^{2} }}} \right. \kern-\nulldelimiterspace} {2\sigma ^{2} }}} \right)$$
    (4)

    where \(p_{0}\) is the baseline detection probability (or magnitude of the detection function). In the half-normal model, the detection probability \(p_{ij}\) decreases monotonically with distance \(d_{ij}\) from \(s_{i}\) (Borchers and Efford 2008; Royle et al. 2014). The spatial scale parameter \(\sigma\) defines the rate of decline in \(p_{ij}\) with distance \(d_{ij}\) from detector j to the \(s_{i}\). By including the spatial information through an explicit model for detection, SCR accounts for one important source of spatial variation in detection probability: the location of an individual’s AC relative to the detectors (Borchers and Efford 2008; Royle et al. 2014).

    We considered the observations of individuals at detectors as the outcome of a Bernoulli process (detections [\(y_{ij} = 1\)] and non-detections [\(y_{ij} = 0\)]) with probability \(p_{ij}\) and conditional on the state \(z_i\) of individual i:

    $$y_{{ij}} \sim Bernoulli\left( {p_{{ij}} {\mkern 1mu} z_{i} } \right)$$
    (5)

    Importantly, this observation model assumes constant baseline detection probability among detectors, and thus does not account for additional detector-specific variation in detectability.

Simulation

General approach

To evaluate the consequences of unmodeled spatial heterogeneity in detectability, we generated SCR data sets with varying patterns of spatial heterogeneity in the baseline detection probability (i.e., detector-specific \(p_{0} = p_{0_j}\)), before fitting SCR models assuming constant baseline detection probability across detectors. We considered two types of scenarios: continuous and categorical spatial variation in detectability. For each scenario, we varied both the level of spatial autocorrelation and the magnitude of the spatial variation to resemble sampling configurations and intensities that may occur in real-life studies (Fig. 1). We also included a reference scenario without spatial heterogeneity in detection probability; thus, the observation sub-model was not misspecified.

Fig. 1
figure 1

Examples of spatially variable and autocorrelated baseline detection probability (higher = darker blue shading) in grid of detectors (gray dots) centered in a habitat (entire area surrounded by the blue line with rounded corners). Shown in rows, spatial variation may be continuous or categorical (with different proportion of area in the lower detectability category). Shown in columns, spatial autocorration may vary from high (Moran’s \(I \approx\) 1) to low (Moran’s \(I \approx\) 0). For a detailed description of each scenario, see the main text.

Full size image

Set up

For all simulations, we used a \(20\times 20\)-distance unit (du) square grid of 400 detectors with 1 du inter-detector spacing. The habitat S included the region covered by detectors and a 4.5-du wide buffer around it (i.e., three times the simulated value for \(\sigma\); Efford 2011) for a total area of \(29\times 29\) du2 (Fig. 1). The buffer allows individuals with ACs located outside but near the detector area to be detected within. We fixed the values for the true population size to \(N = 250\) and the spatial scale parameter of the half-normal detection function to \(\sigma = 1.5\) du across all simulation scenarios for simplicity. We set the size of the augmented population size to be 2.5 times the simulated number of ACs (\(M = 625\)). With this set up, we detected, on average, about \(40\%\) (range = 15–60%) of the true population size N under any given scenario (Table S2.1, Online Appendix 2).

Detector-level covariates

We generated spatially autocorrelated covariates of the baseline detection probability encompassing the extent of the \(20\times 20\) du2 detector grid with the same spatial resolution (1 du) using a function developed by Guélat (2013), with minor modifications (Online Appendix 1). The detector-specific, spatial covariate \({\mathbf {X}}\) was generated using a multivariate normal distribution:

$${\mathbf{X}}\sim MVN\left( {0,\,\Sigma } \right)$$
(6)

where the covariance matrix \(\Sigma\) determines the spatial association between detectors. \(\Sigma\) is calculated using a function D representing the decay in correlation between pairs of detectors with distance. We followed Guélat (2013) and modeled the covariance of \({\mathbf {X}}\) at two detectors j and \(j^{\prime }\) as an exponential decay with distance:

$$D\left( {\delta _{{jj^{\prime}}} } \right) = \exp \left( { - \phi {\mkern 1mu} \delta _{{jj^{\prime}}} } \right)$$
(7)

where \(\delta _{jj^{\prime }}\) is the distance between detectors j and \(j^{\prime }\), and \(\phi\) is the rate determining how rapidly correlation declines with distance. We varied \(\phi\) to simulate covariates with low (\(\phi = 1000\)), intermediate (\(\phi = 1\)), or high (\(\phi = 0.001\)) spatial autocorrelation (Fig. 1). We randomly generated 100 covariate surfaces for each simulation scenario and scaled the resulting values. We then extracted the spatial covariate values for each detector grid cell \({\mathbf {X}}_{j}\) (but see below the extra step for simulating continuous spatial variation in detectability). We quantified the realized spatial autocorrelation using the Moran’s index of global spatial autocorrelation, Moran’s I (Moran 1950; Sokal and Oden 1978; Lichstein et al. 2002) for each simulated spatial covariate with the R package ‘raster’ (Hijmans 2019). I ranges from \(-1\) (perfectly negatively correlated) to 0 (no correlation) and 1 (perfectly positively correlated).

Simulation scenarios

  1. 1.

    Continuous, detector-level variation in detectability—Continuous spatial variation in detection probability may arise in situations where an underlying habitat covariate, such as elevation, forest cover, or distance from a feature (e.g., roads or human settlements) linearly affects baseline detection probability at detectors (Fig. 1) but remains unmodeled. We modeled the detector-specific baseline detection probability \(p_{0_j}\) as a linear function of the simulated spatial covariate \({\mathbf {X}}\):

    $$logit\left( {p_{{0_{j} }} } \right) = \beta _{0} + \beta _{{\mathbf{X}}} {\mkern 1mu} {\mathbf{X}}_{j}$$
    (8)

    where \(\beta _{0}\) is the intercept value of the simulated baseline detection probability \(p_0\) on the logit scale and \(\beta _{{\mathbf {X}}}\) is the regression coefficient of the covariate effect. We kept \(\beta _{0}\) constant across simulations, which corresponds to an intercept value of 0.15 for \(p_0\). We used two values of \(\beta _{{\mathbf {X}}}\) to generate low (\(\beta _{{\mathbf {X}}} = -0.5\)) or high (\(\beta _{{\mathbf {X}}} = -2.0\)) amounts of spatial variation in detection probability. Generating spatial covariates directly as spatially autocorrelated rasters as described above, leads to less tractable outcomes as the level of spatial autocorrelation \(\phi\) influences not only the spatial distribution, but also the density distribution of the covariate. To ensure that the density distribution of spatial covariates on detection probability remained comparable across simulations, regardless of the level of autocorrelation, we mapped a uniformly-distributed spatial covariate with values in the range between \(- 1.96\) and 1.96 onto the spatially autocorrelated similarity raster created as described above (see the code in Online Appendix 1). As a result, the mean covariate value and variance were constant across simulations, regardless of the level of spatial autocorrelation. In other words, we sought stationarity when all simulations of a given scenario are considered jointly, but not within a given simulated detector grid. Following Eq. 8, the detector-specific baseline detection probability \(p_{0_j}\) was between 0.06 and 0.32 (median = 0.15) when the variation in detectability was low (\(\beta _{{\mathbf {X}}} = - 0.5\)). By increasing \(\beta _{{\mathbf {X}}}\) to \(- 2.0\), \(p_{0_j}\) was between 0.003 and 0.9 (median = 0.15; Table S2.1, Online Appendix 2). \(p_{ij}\) was then calculated by reformulating Eq. 4:

    $$p_{{ij}} = p_{{0_{j} }} {\mkern 1mu} \exp \left( {{{ - d_{{ij}}^{2} } \mathord{\left/ {\vphantom {{ - d_{{ij}}^{2} } {2\sigma ^{2} }}} \right. \kern-\nulldelimiterspace} {2\sigma ^{2} }}} \right)$$
    (9)

    Finally, the SCR data \(y_{ij}\) was generated by realizing the detection process following Eq. 5.

    Note that such large differences in \(p_{0_j}\) do not translate into a correspondingly large differences in overall detectability. The probability of detecting any given individual at least once, depends not only on \(p_{0_j}\), but also on the scale parameter of the detection function \(\sigma\), and the position of its AC relative to the entire detector grid. Comparatively extensive spatial differences in baseline detection probabilities have been reported for three large carnivore species in Scandinavia (Bischof et al. 2020a). For example, \(p_{0}\) estimates (after conversion from binomial to Bernoulli) for wolves Canis lupus ranged from near 0 to 0.8, depending on administrative region.

  1. 2.

    Categorical, detector-level variation in detectability—In real-world monitoring studies, this situation arises when there are at least two regions with different sampling intensity across the study area; for example, two regions with varying sampling effort (or different sampling protocols that lead to varying sampling intensity), or two contrasting landscapes (e.g., forest vs. grassland) or detector types (e.g., camera traps on and off trails) across the study area that influence the detection probability (Fig. 1). To represent this situation, we transformed the underlying continuous spatial covariate into a binary one (\({\mathbf {X}}_{j} = 0\) or 1) before simulating scenarios with two classes of detector-specific baseline detection probability \(p_{0_j}\) using Eq. 8. The R code in Online Appendix 1 shows the method used to define the cut-off value to discretize the spatial covariates. We used the same values of \(\beta _{{\mathbf {X}}}\) to generate low or high amount of spatial variation in detectability between the discrete classes of detectors. With this set-up, the baseline detection probability at a group of detectors (\({\mathbf {X}}_{j} = 0\)) was equal to the intercept value (\(p_{0_j} = \beta _{0} = 0.15\)), and \(p_{0_j}\) for the remaining detectors with lower detectability (\({\mathbf {X}}_{j} = 1\)) was either 0.1 (when \(\beta _{{\mathbf {X}}} = -0.5\)) or 0.02 (\(\beta _{{\mathbf {X}}} = -2.0\)) depending on the simulated amount of spatial variation in detectability (Table S2.1, Online Appendix 2). In addition, we considered an extreme case, where a portion of the study area remained entirely unsampled unbeknownst to the investigator, and therefore unaccounted for in the model. This situation could be encountered in, for example, volunteer-based monitoring programs, where SCR data are collected opportunistically but no or only limited spatial information exists about the spatial configuration of sampling. Logistic issues, such as systematic equipment failure (e.g., trap malfunction or damage) or human error, if unreported, may also result in clusters of detectors that remain inactive without the investigator’s knowledge. We set \(\beta _{{\mathbf {X}}}\) to \(-10000\), so that \(p_{0_j} \approx 0\) for detectors with a covariate value of \({\mathbf {X}}_{j} = 1\), and therefore no detection could occur at these inactive detectors.

    To evaluate the potential effect of the proportion of detectors with a different baseline detection probability, we varied the proportion of detectors assigned to each class of the discrete covariate simulated (Fig. 1). We simulated SCR data sets with 25% (n = 100 detectors), 50% (n = 200), or 75% (n = 300) of the detectors assigned to the lower detectability class (\({\mathbf {X}}_{j} = 1\)).

Simulation realization—In total, we generated six simulation scenarios of continuous spatial variation in baseline detection probability with the combinations of two values of \(\beta _{{\mathbf {X}}}\) and three levels of spatial autocorrelation \(\phi\). For the discrete categorical variation in baseline detection probability, we simulated 27 scenarios from all possible combinations of the three values of \(\beta _{{\mathbf {X}}}\), three different proportions of detectors with lower detectability, and three levels of spatial autocorrelation. Finally, we included a scenario of constant detection probability across detectors for reference, i.e. the observation sub-model was correctly specified. We repeated the data simulation procedure 100 times for each combination of parameters, resulting in a total of 3400 simulated SCR data sets (Table S2.1, Online Appendix 2).

Additional simulations—Another common formulation of the detection process in SCR is the Poisson distribution. This models count-based data arising from sampling methods in which individuals can be detected multiple times at a detector during a single occasion, such as in camera trapping studies (Royle et al. 2014). To evaluate the consequences of similar model misspecification in SCR models with a Poisson detection process, we simulated additional count-based SCR data sets and repeated the analysis. We simulated SCR data with counts as observations and fitted a standard, single-session, Poisson-distributed SCR model assuming homogeneous baseline detection rate. A description of the data simulation procedure and fitting of the misspecified Poisson SCR model is provided in Online Appendix 3.

SCR model fitting and evaluation of model performance

We fitted the SCR model described earlier, which does not account for spatial variability in baseline detection probability, to the simulated data sets using NIMBLE (version 0.8.0; de Valpine et al. 2017), nimbleSCR (Bischof et al. 2020b), and R 3.6.1 (R Core Team 2019). We chose vague priors for all primary parameters (\(\psi\), \(\sigma\), \(p_0\), and \(\lambda _0\) as the baseline detection rate in the Poisson SCR model; Table S1.1, Online Appendix 1). We used a local evaluation approach to reduce computation time (Turek et al. 2021). We drew from 3 chains, 15000 Markov chain Monte Carlo (MCMC) samples each, and discarded the initial 5000 samples as burn-in. We visually inspected the mixing of the chains using trace-plots and considered models as converged when all parameters had a potential scale reduction value (Rhat) \(<1.10\) (Brooks and Gelman 1998). We removed from further analysis simulation runs that had failed to reach convergence. R code exemplifying the SCR data simulation for each scenario and model fitting are provided in Online Appendix 1.

To quantify the consequences of unaccounted spatial heterogeneity in detection probability for SCR models, we calculated the relative bias (RB), coefficient of variation (CV), and coverage probability of the 95% credible intervals (hereafter, coverage; Walther and Moore 2005) of the estimates of population size (\({\hat{N}}\)) and spatial scale parameter (\({\hat{\sigma }}\)). RB was calculated as:

$$\widehat{{RB}}\left( \theta \right) = {{\left( {\hat{\theta } - \theta _{0} } \right)} \mathord{\left/ {\vphantom {{\left( {\hat{\theta } - \theta _{0} } \right)} {\theta _{0} }}} \right. \kern-\nulldelimiterspace} {\theta _{0} }}$$
(10)

where \(\hat{\theta }\) is the posterior mean estimate for the parameter of interest and \(\theta _0\) is the true (simulated) value of that parameter. CV was calculated to assess the precision of each parameter estimate as:

$$\widehat{{CV}}\left( \theta \right) = {{\widehat{{SD}}\left( \theta \right)} \mathord{\left/ {\vphantom {{\widehat{{SD}}\left( \theta \right)} {\hat{\theta }}}} \right. \kern-\nulldelimiterspace} {\hat{\theta }}}$$
(11)

where \({\widehat{SD}}(\theta )\) is the standard deviation of the MCMC posterior samples of that parameter. Further, we calculated the coverage as a metric of model fit, which was computed as the proportion of simulation runs in which the 95% credible interval of the parameter estimate included the simulated value of that parameter (Walther and Moore 2005).

Results

Overview

Of the 3400 simulation runs, 3384 (99.5%) reached convergence and were retained: 599 (99.8%) and 2685 (99.4%) for the continuous and categorical scenarios of spatial variation in detectability, respectively (Table S2.2, Online Appendix 2). Our results indicate that the misspecified SCR model was robust to unmodelled spatial heterogeneity in detection probability for most of the scenarios explored. Population size (\({\hat{N}}\)) and spatial scale parameter (\({\hat{\sigma }}\)) estimators exhibited pronounced bias only in the presence of high spatial autocorrelation and high amount of variation in detectability among detectors (Fig. 2). Precision (Fig. 3) and coverage (Fig. 4) were also affected in our extreme scenarios when the detector-specific variation in detection probability remained unmodeled. The consequences of spatial autocorrelation were amplified with increasing proportions of detectors with lower detectability in the categorical scenarios and, in certain situations, when the amount of variation in detection probability was high (Table S2.2, Online Appendix 2). Qualitatively, misspecification of the observation sub-model for Bernoulli- and Poisson-distributed data had comparable consequences. Results for the simulations with the Poisson SCR model are provided in Online Appendix 3.

Fig. 2
figure 2

Relative bias (RB) for population size (\({\hat{N}}\)) and the spatial scale parameter of the half-normal detection function (\({\hat{\sigma }}\)) estimated by a spatial capture-recapture model fitted to simulated data sets. Spatial variation in detection was simulated but not accounted for in the estimation model. Numbers provided on the x-axes refer to the unique ID of the simulation scenario and link with the more detailed information provided for each simulation in Online Appendix 2. Violins show the biasing effects of spatial autocorrelation at the detector-level covariates simulated (decreasing from dark [High] to light colors [Low]) in different scenarios of continuous (top row) and discrete categorical (bottom three rows) variation in baseline detection probability. For the categorical scenarios, the proportions of detectors that belong to a group of detectors with lower detectability increases from top (25%) to bottom (75%). Background colors correspond to three different sets of simulation scenarios with similar values of the magnitude of the covariate effect: Low (\(\beta _{{\mathbf {X}}} = -0.5\); light gray), High (\(\beta _{{\mathbf {X}}} = -2.0\); white), and Extreme (\(\beta _{{\mathbf {X}}} = -10000\); dark gray). Yellow violins labeled as Reference in the top row show the results for a baseline scenario without heterogeneity in detection probability

Full size image

Continuous, detector-level variation in detectability

We observed increasingly negative bias in \({\hat{N}}\) with increasing spatial autocorrelation (Fig. 2; Table S2.2, Online Appendix 2). The magnitude of RB in \({\hat{N}}\) obtained with high spatial autocorrelation (simulation scenarios 1 and 4: median Moran’s I = 0.96; median RB(\({\hat{N}}\)) = \(16\%\)) was substantially greater than that for scenarios with low spatial autocorrelation (scenarios 3 and 6: median Moran’s I = 0; median RB(\({\hat{N}}\)) = \(0.2\%\)). This pattern was amplified when the amount of variation in detectability was high (scenario 4: median RB(\({\hat{N}}\)) = \(30\%\)). In contrast, we detected no noticeable bias in \({\hat{\sigma }}\) (Fig. 2), even when the amount of variation in detectability was high and spatial autocorrelation was high or intermediate (scenarios 4 and 5: median RB(\({\hat{\sigma }}\)) = \(4\%\) and \(7\%\), respectively).

The pattern in precision of \({\hat{N}}\) and \({\hat{\sigma }}\) were almost identical across the scenarios considered (Fig. 3). When the amount of variation in detectability was low (scenarios 1–3), precision of the parameter estimates were comparable to those of the reference scenario without spatial heterogeneity in detection probability (scenario 7: median CV(\({\hat{N}}\)) = \(7\%\) and median CV(\({\hat{\sigma }}\)) = \(5\%\)). Precision of \({\hat{N}}\) and \({\hat{\sigma }}\) were inflated when the magnitude of variation in detectability was high (scenarios 4–6: median CV(\({\hat{N}}\)) = \(5\%\) and median CV(\({\hat{\sigma }}\)) = \(3\%\)), where increasing the spatial autocorrelation slightly decreased the precision of \({\hat{N}}\) (scenario 4; Table S2.2, Online Appendix 2).

Coverage of both \({\hat{N}}\) and, to a lesser extent, \({\hat{\sigma }}\) were drastically impacted by spatial autocorrelation (Fig. 4; Table S2.2, Online Appendix 2). In situations of high spatial autocorrelation and low variability in detection probability, we observed a \(14\%\) reduction in coverage of \({\hat{N}}\), compared to low spatial autocorrelation (from Coverage(\({\hat{N}}\)) = \(93\%\) in scenario 3 to Coverage(\({\hat{N}}\)) = \(80\%\) in scenario 1). The combination of high spatial autocorrelation and high variation in detectability led to coverage of \({\hat{N}}\) as low as \(3\%\) (scenario 4; Fig. 4). The pattern was less pronounced for \({\hat{\sigma }}\) (scenario 4: Coverage(\({\hat{\sigma }}\)) = \(70\%\)). However, we detected a drastic decrease in the coverage of \({\hat{\sigma }}\) for the scenario of intermediate spatial autocorrelation and high amount of variation in detectability (scenario 5: median Moran’s \(I = 0.63\), Coverage(\({\hat{\sigma }}\)) = \(26\%\); Fig. 4).

Fig. 3
figure 3

Coefficient of variation (CV) for population size (\({\hat{N}}\)) and the spatial scale parameter of the half-normal detection function (\({\hat{\sigma }}\)) estimated by a spatial capture-recapture model fitted to simulated data sets. Spatial variation in detection was simulated but not accounted for in the estimation model. Numbers provided on the x-axes refer to the unique ID of the simulation scenario and link with the more detailed information provided for each simulation in Online Appendix 2. Violins show the effects of spatial autocorrelation in the detector-level covariates simulated (decreasing from dark [High] to light colors [Low]) in different scenarios of continuous (top row) and discrete categorical (bottom three rows) variation in baseline detection probability. For the categorical scenarios, the proportions of detectors that belong to a group of detectors with lower detectability increases from top (25%) to bottom (75%). Background colors correspond to three different sets of simulation scenarios with similar values of the magnitude of the covariate effect: Low (\(\beta _{{\mathbf {X}}} = -0.5\); light gray), High (\(\beta _{{\mathbf {X}}} = -2.0\); white), and Extreme (\(\beta _{{\mathbf {X}}} = -10000\); dark gray). Yellow violins labaled as Reference in the top row show the results for a baseline scenario without heterogeneity in detection probability, which was included for reference. The dashed lines show median CV for the parameter estimates achieved in the reference scenario

Full size image

Categorical, detector-level variation in detectability

We observed increasing negative bias in \({\hat{N}}\) with increasing spatial autocorrelation and increasing proportion of detectors with lower detectability (Fig. 2; Table S2.2, Online Appendix 2). The bias was particularly pronounced in scenarios of extreme, spatially autocorrelated variation in detectability (median Moran’s I = 0.86), where 50% or more of detectors were inactive, ultimately leading to an entire region with little or no sampling (scenarios 23 and 32: median RB(\({\hat{N}}\)) = \(49\%\)). However, when the variation in detectability was low (scenarios 8–10, 17–19, and 26–28), \({\hat{N}}\) was minimally affected in terms of bias regardless of the spatial autocorrelation and proportions of detectors with lower detectability (Table S2.2, Online Appendix 2). Similarly, when spatial autocorrelation in detectability was low, \({\hat{N}}\) remained relatively unbiased even in the extreme scenarios of a portion of detectors being inactive. We observed no systematic bias in \({\hat{\sigma }}\) under any of the scenarios tested (Fig. 2; Table S2.2, Online Appendix 2).

Precision of \({\hat{N}}\) and \({\hat{\sigma }}\) decreased with increasing proportions of detectors with lower detectability and increasing variation in detectability (Fig. 3; Table S2.2, Online Appendix 2). However, precision increased with spatial autocorrelation. For the extreme scenario of a partially-sampled area, when 75% of the detectors were inactive and spatial autocorrelation was low (scenario 34; Fig. 3), median CV(\({\hat{N}}\)) and median CV(\({\hat{\sigma }}\)) increased by \(249\%\) and \(242\%\), respectively, compared to the reference scenario. In contrast, the increase in median CV(\({\hat{N}}\)) and median CV(\({\hat{\sigma }}\)) was \(129\%\) and \(145\%\), respectively, for the scenario of high spatial autocorrelation (scenario 32; Table S2.2, Online Appendix 2).

The coverage of \({\hat{N}}\) drastically decreased with increasing spatial autocorrelation (Fig. 4). With high spatial autocorrelation, the coverage of \({\hat{N}}\) dropped to between \(39\%\) and \(44\%\) in the presence of high variation in detectability and higher proportions of detectors with lower detectability (scenarios 29 and 20). When more than 50% of the detectors were inactive, coverage of \({\hat{N}}\) was near 0 (scenarios 23 and 32; Fig. 4). Coverage was nominal for \({\hat{\sigma }}\) under low and high levels of spatial autocorrelation, and only decreased to between \(79\%\) and \(87\%\) when spatial autocorrelation was intermediate (scenarios 15, 24, and 33: median Moran’s I = 0.42; Table S2.2, Online Appendix 2).

Fig. 4
figure 4

Coverage probability of the 95% credible intervals of population size (\({\hat{N}}\)) and the spatial scale parameter of the half-normal detection function (\({\hat{\sigma }}\)) estimated by a spatial capture-recapture model fitted to simulated data sets. Spatial variation in detection was simulated but not accounted for in the estimation model. Numbers provided on the x-axes refer to the unique ID of the simulation scenario and link with the more detailed information provided for each simulation in Online Appendix 2. Points show the effects of spatial autocorrelation of continuous (top row) or discrete categorical (bottom three rows) baseline detection probabilities (decreasing from dark [High] to light colors [Low]). For the categorical scenarios, the proportions of detectors that belong to a group of detectors with lower detectability increases from top (25%) to bottom (75%). Background colors correspond to three different levels of variation in the baseline detection probability: Low (\(\beta _{{\mathbf {X}}} = -0.5\); light gray), High (\(\beta _{{\mathbf {X}}} = -2.0\); white), and Extreme (\(\beta _{{\mathbf {X}}} = -10000\); dark gray). Yellow points labeled as Reference in the top row show results for a baseline scenario without heterogeneity in detection probability

Full size image

Discussion

Our study revealed that unmodeled spatial variation in detection probability, which is ubiquitous in wildlife monitoring studies, can have pronounced consequences for inferences from SCR analyses. The critical factor is spatial autocorrelation in detection probability: highly autocorrelated and highly variable detectability leads to pronounced bias and reduction in precision and coverage of SCR estimates of population size, the main parameter of interest in such studies (Royle et al. 2014, 2018). Conversely, at low levels of spatial autocorrelation, SCR model estimates remained robust even to extremely high levels of unmodeled spatial heterogeneity in detection probability. Estimates of the spatial scale parameter of the half-normal detection function were, however, unbiased and coverage remained comparatively high for most scenarios of spatially autocorrelated detection probability. Nonetheless, the pattern of reduction in precision of the estimates of the scale parameter was similar to those in the estimates of population size.

SCR is a powerful analytical tool for spatially explicit inference on the ecology of wild populations while accounting for imperfect detection (Royle et al. 2018). One of the advantages of SCR models over non-spatial capture-recapture models is the ability to address the effects of known spatial variation in detectability as a result of the spatial configuration of ACs relative to detectors and detector-specific characteristics; thus, providing a more realistic model of the detection process (Efford et al. 2013; Royle et al. 2014). Yet, despite the spatially explicit nature of SCR data collection and analysis, it is liable to be afflicted by unknown or undocumented variation in detection probability; for example, due to inadvertent (e.g., equipment failure, within-individual variations) or design-induced factors (e.g., non-random or preferential sampling), or as the result of pooling or integrating data across multiple survey methods or detectors (Efford et al. 2013; Howe et al. 2013; Efford and Mowat 2014; Royle et al. 2013, 2014; Gerber and Parmenter 2015). Consistent with our predictions, failure to adequately account for variable and spatially autocorrelated detection probability may bias SCR parameter estimates, impact estimates of precision, and generally lead to erroneous inferences. In extreme cases among the ones we considered, this led to up to \(65\%\) bias and to virtually zero probability that the credible interval contains the true value of population size (Table S2.2, Online Appendix 2). Comparable results are reported when individual heterogeneity remains unaccounted for in both spatial and non-spatial capture-recapture (Royle et al. 2013; Gimenez et al. 2018a, Stevenson et al. 2021). Thus, a misspecified detection model may not only impact inferences, but also reduce the chance of meeting the management and conservation objectives that often motivate these studies.

We observed the strongest systematic bias in cases where a substantial and distinct region of the detector space was unsampled, essentially leaving a hole in the detector grid that remained unaccounted for in the model. The pronounced negative bias in \({\hat{N}}\) that arises in these situations can be explained if we consider that individuals with ACs within such a hole have a high chance of being missed entirely during sampling, essentially creating a latent class of individuals with very low or zero detectability. This is the case when a portion of detectors are assumed to be active or searched when they are not. As a consequence, detection probability estimates are biased high as they are based primarily on recaptures of individuals detected in active regions of the detector grid, and estimates of abundance and density are consequently biased low. This pattern is amplified by the magnitude of variation in detection probability between regions of lower and higher detectability, and by the size of the region with lower detection probability (Table S2.2, Online Appendix 2).

Situations where entire regions of the study area remain unsampled are fairly common in monitoring studies, where the data are gathered opportunistically, either in part or as a whole (Conn et al. 2017; Altwegg and Nichols 2019; Sicacha-Parada et al. 2021). Unstructured sampling can augment information, and thus improve population inferences about rare or elusive species (Thompson et al. 2012; Tenan et al. 2017; Sun et al. 2019; Bischof et al. 2020a). Opportunistically-collected data obtained as part of public surveys (e.g., citizen-science data) are sometimes integrated into monitoring programs as this allows investigators to sample areas at unprecedented scales with lower costs and the added benefit of public involvement in management and conservation practices (Altwegg and Nichols 2019; Bischof et al. 2020a; Sicacha-Parada et al. 2021). To minimize and account for variation in detection probability in such sampling schemes, volunteers should be encouraged to visit all habitats within the study area, reduce variability in observer proficiency by providing standardized training, and collect data on potentially relevant covariates (Altwegg and Nichols 2019).

We considered a wide range of scenarios in terms of the magnitude and spatial configuration of variability in detection probability. Our extreme scenarios of unknown and unmodeled, high spatial variation and high autocorrelation in \(p_{0}\) are probably less common, but plausible, in real-life SCR monitoring studies. These scenarios were specifically motivated by our previous work on landscape-scale monitoring of large carnivores, where spatial variation in effort can be challenging to quantify because the data are collected in a structured fashion by field personnel associated with multiple jurisdictions and opportunistically by volunteer members of the public (Bischof et al. 2016, 2017, 2020a; Milleret et al. 2020). These SCR studies revealed significant differences in baseline detection probability across administrative units (i.e., counties in Sweden and large carnivore management regions in Norway). Even after controlling for spatio-temporal variations in sampling effort by including detector-level and individual covariates on the baseline detection probability, spatial variation in detectability remained substantial, presumably linked with unmeasured regional differences in search effort by volunteers and sampling configuration implemented by the authorities (Bischof et al. 2020a). Nevertheless, we expect the consequences of such an unmodeled spatial heterogeneity to be relaxed when detection probability is relatively high. In such situations, unknown and unmodeled variation in detectability is likely to be less problematic as individuals from the population are more likely to be detected at multiple detectors (i.e., increase in both the proportion of detected individuals from the population and spatial recaptures). The effects of spatial heterogeneity in detectability may also be mediated by species space-use characteristics. Here, we simulated a target population with intermediate home-range overlap (\(k = \sigma \sqrt{Density} = 0.83\); Efford et al. 2016) with a density of 0.3 AC/du2, which was constant across simulation scenarios. For species with small home-range sizes, where the spatial scale parameter is small relative to the distance between detectors, unmodeled spatially autocorrelated detection probability may lead to more problematic inferences due to the increased risk of individuals going completely undetected in the resulting “holes” in the detector grid. In addition, we assumed homogeneous population density across the study area, which is rarely, if ever, the case in wild populations. Previous studies have suggested that the consequences of unmodeled heterogeneity in detection probability could be amplified when detectability is correlated with density (e.g., Clark 2019; Paterson et al. 2019). As a positive correlation between animal density and detection probability results in a positive bias in estimates of population size, we expect that the biasing effects we detected would be comparatively less pronounced when the unmodeled detection probability is spatially autocorrelated but the effects of unmodeled spatial heterogeneity would not be mitigated. The misspecification of the observation sub-model may also lead to confounding effects between variation in detection probability and animal density that must be accounted for in order to avoid flawed inferences (Efford et al. 2013).

When the drivers of spatial variation in detection probability are known or suspected, the variation can be accounted for by, for example, including relevant covariates (e.g., transect length, camera trap nights, habitat or line-of-sight visibility) in the observation sub-model of SCR models (Royle et al. 2009, 2014; Efford et al. 2013; Efford and Mowat 2014; Bischof et al. 2017, 2020a). In the absence of proxies that can serve as fixed or random effects, how should investigators deal with unknown spatial heterogeneity in detectability in SCR analyses? One solution would be to discard affected data, such as observations contributed by members of the public without corresponding measures of search effort. This, however, would lead to a loss of potentially valuable information, both in terms of number of individuals detected and spatial recaptures (Marques et al. 2011; Sollmann et al. 2012; Tourani et al. 2020b). Explicitly modeling spatially autocorrelated detectability, for example using mixture methods (e.g., by drawing baseline detection probability at each detector from a finite mixture of distributions; Royle 2006) or conditional autoregressive (CAR) models, may offer a solution. CAR models have been implemented in other hierarchical modeling frameworks, such as non-spatial capture-recapture and spatial distribution modeling (e.g., Johnson et al. 2013; Chen and Ficetola 2019; Nicolau et al. 2020), but we are not aware of similar extensions in SCR studies. The added complexity resulting from the inclusion of an autoregressive component on a latent variable like detection probability could pose a significant computation barrier to implementation in Bayesian SCR at large spatial scales. However, recent advances in both software (de Valpine et al. 2017; Bischof et al. 2020b) and Bayesian SCR model formulations (Milleret et al. 2018; Turek et al. 2021) have improved model fitting efficiency, allowing fitting SCR models of unprecedented complexity and spatial scales (Bischof et al. 2020a); thus, opening new possibilities.

Goodness-of-fit tests could help practitioners diagnose potential violations of model assumptions, including unmodeled spatial variation in detection probability, and determine whether there is a need to account for it in the model in the first place. Such diagnostics should be an integral part of any ecological modeling exercise (Conn et al. 2018), and Bayesian p-values (Gelman et al. 1996) have been proposed as a general framework for goodness-of-fit testing of Bayesian SCR models (Royle et al. 2014). However, goodness-of-fit diagnostics for SCR is a developing field of research and specific tests and recommendations, such as those that were developed for non-spatial capture-recapture models (e.g., Gimenez et al. 2018b) are still lacking. Thus, testing and accounting for possible violations, as well as correcting model estimates, is still challenging in SCR.

Conclusions

Unmodeled spatially autocorrelated variation in detection probability can noticeably impact the reliability of inferences derived using SCR, when variability in detection and spatial autocorrelation is high. This specifically affects estimates of abundance and density, primary parameters in wildlife monitoring studies, and can therefore have severe consequences for wildlife management and conservation. Unobserved spatial variability in detection probability is likely ubiquitous in real-life SCR studies, and we encourage research to develop approaches that help practitioners diagnose and account for it.