A hierarchical model for estimating the spatial distribution and abundance of animals detected by continuous-time recorders

Abstract

Motivation

Several spatial capture-recapture (SCR) models have been developed to estimate animal abundance by analyzing the detections of individuals in a spatial array of traps. Most of these models do not use the actual dates and times of detection, even though this information is readily available when using continuous-time recorders, such as microphones or motion-activated cameras. Instead most SCR models either partition the period of trap operation into a set of subjectively chosen discrete intervals and ignore multiple detections of the same individual within each interval, or they simply use the frequency of detections during the period of trap operation and ignore the observed times of detection. Both practices make inefficient use of potentially important information in the data.

Model and data analysis

We developed a hierarchical SCR model to estimate the spatial distribution and abundance of animals detected with continuous-time recorders. Our model includes two kinds of point processes: a spatial process to specify the distribution of latent activity centers of individuals within the region of sampling and a temporal process to specify temporal patterns in the detections of individuals. We illustrated this SCR model by analyzing spatial and temporal patterns evident in the camera-trap detections of tigers living in and around the Nagarahole Tiger Reserve in India. We also conducted a simulation study to examine the performance of our model when analyzing data sets of greater complexity than the tiger data.

Benefits

Our approach provides three important benefits: First, it exploits all of the information in SCR data obtained using continuous-time recorders. Second, it is sufficiently versatile to allow the effects of both space use and behavior of animals to be specified as functions of covariates that vary over space and time. Third, it allows both the spatial distribution and abundance of individuals to be estimated, effectively providing a species distribution model, even in cases where spatial covariates of abundance are unknown or unavailable. We illustrated these benefits in the analysis of our data, which allowed us to quantify differences between nocturnal and diurnal activities of tigers and to estimate their spatial distribution and abundance across the study area. Our continuous-time SCR model allows an analyst to specify many of the ecological processes thought to be involved in the distribution, movement, and behavior of animals detected in a spatial trapping array of continuous-time recorders. We plan to extend this model to estimate the population dynamics of animals detected during multiple years of SCR surveys.

About the Authors:

Robert M. Dorazio

* E-mail: [email protected]

Affiliation: Wetland and Aquatic Research Center, U.S. Geological Survey, Gainesville, Florida, United States of America

ORCID http://orcid.org/0000-0003-2663-0468

K. Ullas Karanth

Affiliations Wildlife Conservation Society, Bronx, New York, United States of America, National Centre for Biological Sciences, Bangalore, Karnataka, India, Centre for Wildlife Studies, Bangalore, Karnataka, India

Introduction

Spatial capture-recapture (SCR) models provide a statistical framework for estimating the density of individuals in a population, that is, the number of individuals per unit area in a finite region of two-dimensional, Euclidean space [1, 2]. In contrast to non-spatial models, SCR models use the detection histories and the locations of marked individuals. An important component of SCR models is the spatial point process that specifies the location of each individual’s center of activity (or home range in some cases). The activity center of an individual is not directly observable, but it is assumed to be fixed during the period of sampling. This latter assumption can be regarded as a form of demographic closure for SCR models.

Some of the earliest SCR models were developed for the analysis of detections of naturally marked animals at a set of fixed locations [3, 4]. Natural markings allow individuals to be identified uniquely and come in a variety of forms owing to technological advances in non-invasive sampling methods. For example, body markings (stripe or spot patterns) can be used to distinguish individual tigers or leopards, and acoustic profiles can be used to distinguish the vocalizations of individual songbirds or whales (at least for some species). Alternatively, DNA obtained from scats or hair snares can be used to identify individuals.

Different types of natural markings obviously require different types of detectors. Here we focus on SCR surveys that satisfy two design requirements: (1) individuals can potentially be observed at more than one of several fixed locations during the period of sampling and (2) the date and time of each observation is recorded, as is common when using microphones or motion-activated cameras as detectors [5, 6]. To satisfy the first requirement, detectors should be placed at locations based on the likely range of an individual’s movements. These placements induce heterogeneity in the rates of detection among individuals because animals whose activity centers lie within the interior of a spatial array of detectors are more likely to be observed than animals whose activity centers lie near the edge of the array.

Although several SCR models have been developed with the first design requirement in mind, most these models have ignored the continuous period of operation of the detectors. In some models the total period of sampling is partitioned into discrete intervals and a binary observation of 0 (no detections of an individual) or 1 (one or more detections of an individual) is used for each interval ([4] and [2], chapter 5). This discretizaton of time is subjective and potentially uses data inefficiently because the same individual can be detected multiple times during a single interval. Discretization into 24-hour intervals does not solve the problem because two detections of the same individual may be counted either once or twice depending on the choice of cut-off time (the so-called “midnight problem” [7]). Other SCR models use the total number of detections of an individual by treating the total period of sampling as a single interval ([3] and [2], chapter 9); however, these models do not use information associated with the exact dates and times of detection.

To our knowledge, only one SCR model has been formulated for times of detection [6], and this model was developed primarily for the estimation of population density. Here, we propose a hierarchical Bayesian model which includes many of the assumptions made by [6] but also allows inferences to be computed for the spatial distribution of individual activity centers. Our approach may be used to construct species distribution models. It is useful even when the sources of spatial variation in population density are either unobserved or unknown.

We begin by describing the hierarchical model, which includes a spatial point process for the distribution of individual activity centers and a temporal point process for the distribution of detection times of each individual. We show that this model is a logical extension of count-based SCR models that do not account for temporal heterogenity in detection rates. We illustrate the benefits of our approach by analyzing the detection times of individual tigers observed in a spatial array of camera traps.

Materials and methods

Model description

Assume a population of N individually identifiable animals living within a finite region B of Euclidean space (specifically, ) is surveyed using K camera traps whose locations are selected to be representative of B. The trap locations need not be random but their spacing should permit any animal to be detected by more than one camera as a consequence of its movements. We assume that each animal living within B moves around a fixed, center of activity, which we denote using s. Our goal is to infer the spatial distribution and abundance of the activity centers of all animals living within B using only the detection locations and times of animals observed during the camera-trap survey.

The data.

To state the inference problem more concretely, we first describe the kinds of data observed in the camera-trap surveys. Suppose each of K distinct cameras is used to detect animals in region B during a period of time when the activity center of each animal can be assumed to be fixed. This assumption implies a spatial form of demographic closure for the population of animals living within region B. Ideally, all cameras operate during the same period; however, the time required to deploy all K cameras usually induces some differences among their actual periods of operation. Therefore, let the time interval (0, T_k] denote the continuous period of operation of the kth camera, and let x_k ∈ B denote the fixed location of this camera (k = 1, …, K).

Each camera is equipped with a motion sensor so that a photograph (detection) is only obtained during an encounter with one or more animals. The time and date of the encounter are normally recorded with each photograph; therefore, it is possible to observe a sequence of detection times for each marked animal and trap. Let denote the observed sequence of detection times of the ith individual at trap k, where t_ikj ∈ (0, T_k] if y_ik > 0; otherwise, t_ik = ∅. The subscript i ∈ {1, …, n} indexes the ith of n distinct animals detected at least once during the camera-trap survey. (The prime symbol indicates the transpose of a matrix or vector.)

The observed data can be summarized using two arrays: one is a n × K matrix that contains the observed number of detections y_ik of the ith individual at trap k; the other is a n × K matrix that contains the observed sequence of detection times t_ik for the ith individual at trap k. All that is known of the N − n unobserved animals is that each individual’s detection frequency is zero (Fig 1).

[Figure omitted. See PDF.]

Fig 1. Notation used for the data observed during a camera-trap survey.

https://doi.org/10.1371/journal.pone.0176966.g001

If spatially varying covariates are thought to influence the density of individual activity centers over B, we assume that a (possibly vector-valued) set of covariate measurements v(s) is available for every location s in region B. In the absence of such covariates, a model of constant density can be specified using v(s) = 1 as a regressor (see below). Similarly, if temporally varying covariates are thought to influence the temporal pattern of detection times recorded at the kth camera, we assume that a (possibly vector-valued) set of covariate measurements z(t_k) is available for every time t_k ∈ (0, T_k] that the kth camera is operational.

In some surveys one or more covariates thought to be informative of an individual’s baseline rate of detection is available for each of the K cameras. Let w_k denote the (possibly vector-valued) set of covariate measurements associated with the camera located at x_k. These covariates are location dependent and can be interpreted as predictors of the baseline rate at which individuals encounter the kth camera during the time interval (0, T_k]. In the absence of such covariates, the baseline rate of encounter can be specified using w_k = 1 (∀k) as a regressor (see below).

Spatial point-process model of individual activity centers.

Following previous spatial capture-recapture models [8–10], we assume that the spatial distribution of activity centers is a realization of a Poisson process parameterized by first-order intensity function λ(s), which specifies the limiting, expected density of activity centers (number per unit area) at location s ∈ B [11]. Thus, the total number N(B) of individuals in region B is a Poisson random variable that depends on the intensity of the process integrated over B: Λ(B) = ∫_B λ(s) ds. We model the expected density λ(s) as a log-linear function of unknown parameters β and location-specific regressors v(s) as follows: log[λ(s)] = β′v(s).

Let N denote the unknown total number of individuals living in region B. We have established that(1)Furthermore, conditional on N, the joint probability density of the N individual activity centers is(2)[12, 13]. (Note that we use bracket notation [14] to specify probability density functions; thus, [x, y] denotes the joint density of random variables X and Y, [x|y] denotes the conditional density of X given Y = y, and [x] denotes the unconditional (marginal) density of X.) The goal of our analysis is to estimate N and s_i (i = 1, …, N) using only the observed data—that is, using the detection times of the n ≤ N individuals photographed by the cameras and their corresponding capture locations.

Temporal point-process model of detection times.

In this section we describe a model of the observed data (i.e., the trap-specific detection times of each individual) that depends on the latent activity center of each individual. As with most spatially explicit models of camera-trap surveys, we assume that pairwise distances between each individual’s activity center and the set of camera locations induces heterogeneity in detection rates among the animals living in B. For example, animals whose activity centers are located in the interior of the trapping array are more likely to be detected than animals whose activity centers are located outside the trapping array. However, our model of detection times includes additional components that can influence the detection rates of individual animals. To be specific, we assume that the sequence of detection times observed for the ith individual at camera trap k is a realization of a Poisson process parameterized by first-order, intensity function(3)for t ∈ (0, T_k], where ∥s_i − x_k∥ specifies the Euclidean distance between locations s_i and x_k and σ (> 0) is an unknown parameter whose value increases with the spatial extent an animal’s movements about its activity center. To specify the baseline rate ψ_k (> 0) at which individuals encounter the kth camera trap during time interval (0, T_k], we use a log-linear model log(ψ_k) = α′w_k, which includes a vector of unknown parameters α and an observed set of covariate measurements w_k associated with the kth trap. In addition, the function γ(t) may be used to specify sources of time-dependence in the sequence of each animal’s detection times. For example, if animals are primarily active at night (as with tigers in India), one would expect encounters of animals with traps to be higher during the night than during periods of daylight. In this case a time-dependent indicator of daylight, say z(t), might be used to model γ(t) in terms of an unknown parameter ξ as follows: log(γ(t)) = ξ z(t). Of course, this approach may be generalized to specify the effects ξ of several temporal covariates z(t). An important caveat is that the value of each temporal covariate must be available for the entire period of each trap’s operation.

We have not yet described the relationship between the observed detection times of an individual and the Poisson intensity given in Eq (3). Conditional on activity center s_i, the number of detections Y_ik(T_k) of the ith individual at trap k during time interval (0, T_k] has a Poisson distribution as follows:where is the expected number of detections of the ith individual at trap k during its period of operation. Furthermore, conditional on the observed number of detections y_ik and s_i, the joint probability density of the observed detection times isNote that unlike the spatial Poisson process, there is only one possible ordering of the observed detection times in t_ik (i.e., sequential order).

Likelihood functions for fitting models to data.

We have described all of the components needed to derive the likelihood function of the model. To clarify the derivation of this function, for the moment suppose (counterfactually) that the activity centers of all N individuals living in region B were observable and that we wanted only to estimate the unknown parameters θ = (β′, α′, σ, ξ′)′. In this case the likelihood function for θ is(4)where y_(1:N) and t_(1:N) are N × K matrices with elements y_ik and t_ik, respectively. The joint density of y_ik and t_ik (conditional on s_i) is [y_ik, t_ik|s_i] = [y_ik|s_i][t_ik|y_ik, s_i].

Of course, in reality only n of the N individuals living in region B are detected in the array of camera traps, and none of the individual activity centers is observable. All that is known about the n₀ = N − n unobserved individuals is that each animal was not detected at any of the K traps, that is y_ik = 0 (∀k, i ∈ {n + 1, …, N}). The likelihood function of the observed data requires the following modification of Eq (4) to account for the ways of observing n individuals in the camera-trap surveys:(5)where y_(1:n) and t_(1:n) are the n × K matrices of observed data and whereThe probability that an individual living in region B was not detected in any of the K camera traps during the survey period isNote that the likelihood function in Eq (5) includes some quantities that are not observed (n₀ and s_i). If desired, a marginal likelihood function for θ can be computed by integrating n₀ and s_i from Eq (5) (see Dorazio 2013 [10] for example). However, for our purposes this calculation is counterproductive because we want to estimate the activity centers of both observed and unobserved individuals living in B. We therefore adopt the Bayesian approach to inference and fit the model using the following unnormalized posterior density function:(6)where [θ] denotes the prior density of θ.

Restricted models for a homogenous temporal point process.

In the absence of time dependence in the detection process (i.e., γ(t) = γ, ∀t), the spatial distribution of individual activity centers can be estimated without using the detection times t_(1:n) of the n observed individuals. This is possible because the temporal point process is homogeneous and the intensity function in Eq (3) simplifies towhich is not a function of time t.

In this case the expected number of detections of the ith individual at trap k is proportional to the duration T_k of that trap’s operation, that is, Φ(T_k, s_i, x_k) = T_k ϕ(s_i, x_k). Furthermore, it is easily shown that the conditional distribution of observed detection times is Uniform(0, T_k). This result implies a considerable simplification of the likelihood function because the joint density of y_ik and t_ik (conditional on s_i) does not depend on the observed sequence t_ik of detection times:Note that the parameters used to specify ψ_k and γ (i.e., α and ξ, respectively) are no longer identifiable; thus, we can ignore γ—recognizing that it cannot be identifed from the parameter α₀—and formulate a likelihood function [y_(1:n), n, n₀, s₁, …, s_n|θ] of the observed matrix of counts y_(1:n) to estimate θ = (β′, α′, σ)′. As before, we adopt the Bayesian approach to inference and fit this restricted class of models using the following unnormalized posterior density function:(7)

Posterior inference and predicting activity centers of unobserved individuals.

We used Markov chain Monte Carlo (MCMC) methods [15] to estimate the parameters of models with posterior density functions given in Eqs (6) and (7) and to estimate maps of the spatial distribution of activity centers. Once a model was fitted, we estimated the activity centers of the N − n unobserved individuals by simulating draws from a common distribution whose conditional density function isBy simulating draws of s for each element of the Markov chain, we averaged over the posterior uncertainty of the model’s parameters to compute summaries (including maps) of the posterior predictions of these activity centers. For example, we estimated the spatial distribution of abundance of individuals by partitioning the spatial domain B into a set of disjoint (nonoverlapping) grid cells (B_m for the mth cell) and by computing the number of individual activity centers in each grid cell as follows: , where I(ε) denotes an indicator function that equals one if expression ε is true and equals zero otherwise. Abundance N(B_m) was computed for each element of the Markov chain so that posterior uncertainty of all model parameters was accounted for while estimating each cell’s abundance. Details about our MCMC algorithm and prior distributions are described in S1 Appendix.

Camera-trap survey of tigers

We surveyed the population of tigers (Panthera tigris) living within and adjacent to the Nagarahole Tiger Reserve of Karnataka, India using a spatial array of 162 motion-activated camera traps (Fig 2). Ecological descriptions of the study area and species were provided by [16, 17] and [18]. Each camera trap included two inward facing cameras so that both flanks of an individual tiger could be photographed [19]. Camera trap locations were selected to increase detections using prior knowledge of the behavior and range of movements of individual tigers [4, 17]. Each camera trap operated continuously during the 45-day period of sampling (26 Nov 2014 to 13 Jan 2015), recording the date and time when a tiger was photographed. Individual tigers were identified from photographs using their unique stripe patterns [16], as shown in Fig 3. All field work was carried out under research permits issued by the Karnataka Forest Department and the Kerala Forest Department, in forest areas administered by these agencies. No other approvals or permits were required because our non-invasive, camera-trap surveys did not involve any handling, collection, sampling or sacrifice of animals, or experimental manipulation of any type.

[Figure omitted. See PDF.]

Fig 2. Map of Nagarahole Tiger Reserve and adjacent areas in the Western Ghat region of India (shown in green).

Locations of camera traps are indicated by plus signs. The buffer region outside the study area (shown in blue) also contains suitable tiger habitat.

https://doi.org/10.1371/journal.pone.0176966.g002

[Figure omitted. See PDF.]

Fig 3. Photograph of tiger detected at a camera trap within the Nagarahole Tiger Reserve.

https://doi.org/10.1371/journal.pone.0176966.g003

Camera traps were placed 2–3 km apart throughout the 860 km² area to provide adequate exposure to all tigers. The spatial domain used to fit our SCR models also included a 10-km wide buffer where suitable habitat existed, thus yielding a total area of 1130 km² (Fig 2). The width of this buffer was based on previous radio-telemetry and camera-trap studies of tiger movements and home ranges [17, 18]. Areas outside of this domain contained land-use formations unlikely to be occupied by resident tigers.

We fitted two SCR models to the detection times of individual tigers observed in the spatial array of cameras. One model used a time-dependent binary indicator z(t) to specify whether each detection occurred during the daytime (i.e., between sunrise and sunset (z(t) = 1)) or during the nighttime (z(t) = 0). We obtained sunrise and sunset times for each day of sampling from an online application http://aa.usno.navy.mil/data/docs/RS_OneYear.php developed by the Astronomical Applications Department of the United States Naval Observatory. We also fitted a restricted model wherein detection rate was assumed to be constant at each camera location, as described in Model description. This model did not require the detection times and was fitted using only the matrix of observed detection frequencies. Neither of these SCR models included spatial covariates of tiger density or trap-specific covariates of tiger detection rate; however, we computed estimates of the spatial distribution of individual activity centers and tiger abundances using the full (unrestricted) model and a 500 m × 500 m grid of the spatial domain. The R source code [20] and data needed to fit these models are available in S2 Appendix.

Simulation study

We conducted a simulation study to examine the performance of our model when analyzing data sets of greater complexity than the tiger data. In this study we simulated values of one spatially varying covariate v(s) of individual density and one spatially varying covariate w(s) of baseline detection rate over a square region B (Fig 4). S3 Appendix contains a description of how values of each covariate were computed as a function of location s. The expected density of individuals at location s was computed using λ(s) = exp(1.4 + 0.8v(s)). Given these parameter values, the expected number of individuals in region B was 150.2 (= Λ(B)).

[Figure omitted. See PDF.]

Fig 4. Maps of covariates used in simulation study.

Upper panel = covariate of individual density. Lower panel = covariate of baseline detection rate. Trap locations are indicated by circles.

https://doi.org/10.1371/journal.pone.0176966.g004

We assumed that individuals living in B may have been detected within a square array of 100 evenly-spaced traps located at the center of the region (Fig 4). The baseline detection rate of individuals at trap location x_k was computed using ψ_k = exp(−0.7 + 1.0w(x_k)). Given these parameter values, baseline detection rates ranged from 0.41 to 0.61 detections per unit time among the K = 100 trap locations. We assumed each trap operated continuously for 30 units of time. We also assumed equal amounts of daytime and nighttime during this period and that each trap’s baseline detection rate during the day was about 37% of the rate at night. In other words, we assumed ξ = −1. We assumed the parameter σ was 0.4 to ensure that individuals in proximity to the trapping array were likely to be detected at multiple traps.

We fitted the full model and the restricted model to each of 300 data sets that were simulated independently. Each model was fitted to a data set by running the MCMC algorithm (S1 Appendix) for 1200 iterations and by estimating posterior summaries from the final 1000 iterations of the Markov chain. We examined the performance of these models by comparing the average of the estimated posterior mean of each parameter to its true (data-generating) value. We also compared the estimated coverage of a 95% credible interval for each parameter to the nominal 95% level.

Results

Camera-trap survey of tigers

A total of 86 individual tigers were detected during the period of sampling, which included a per-camera average of 21.2 days of daytime and 23.4 days of nighttime. Tigers were detected in a total of 355 capture events. Most of these detections occurred during the nighttime (Fig 5), providing evidence of nocturnal activity patterns of tigers.

[Figure omitted. See PDF.]

Fig 5. Temporal distribution of observed detection times for photo-captured tigers in a 24-hour period.

https://doi.org/10.1371/journal.pone.0176966.g005

Table 1 contains estimates of the parameters of both full and restricted SCR models. Substantial differences in parameter estimates were obtained only among the parameters that specify encounter rates (α₀ and ξ). In particular the estimate of ξ was significantly negative, indicating that baseline detection rates of tigers were about five times lower during the daytime (0.009 detections per day) than during the night (0.045 detections per day).

[Figure omitted. See PDF.]

Table 1. Estimates of posterior summary statistics for the parameters of the SCR models fitted to the tiger data.

https://doi.org/10.1371/journal.pone.0176966.t001

The expected density of tigers throughout our spatial domain was 11.3 individuals per 100 km² (95% credible interval = 9.1–13.9). Although spatially varying covariates of tiger density were not used to fit our SCR model, the estimated spatial distribution of tiger abundance (Fig 6) revealed some locations where tiger densities were much higher than average.

[Figure omitted. See PDF.]

Fig 6. Map of estimated distribution of tiger abundance in Nagarahole Tiger Reserve and adjacent areas in the Western Ghat region of India.

https://doi.org/10.1371/journal.pone.0176966.g006

The estimate of the scale parameter σ was 1.70 km (95% credible interval = 1.60–1.82). Using a Gausian kernel, we estimate that each tiger occupied an area of about 54 km² (= ), which slightly exceeds the home range expected for resident male tigers (50 km²) [18].

Simulation study

In our simulated data sets the average number of individuals living in region B (150.3) was very close to the expected population size (150.2). The average number of individuals detected at least once in the trapping array was 124.0. Table 2 contains averages of each parameter’s posterior mean and coverage estimated by fitting the full and the restricted SCR models to each simulated data set. The average of the estimated posterior mean was close to the true (data-generating) value of every parameter with one exception. The average of the logit-scale, mean baseline encounter rate α_o estimated under the restricted model was substantially lower than the true value of this parameter. Similarly, the coverage of credible intervals of α₀ obtained by fitting the restricted model was substantially lower than the nominal 95% level. The reason for these discrepancies are obvious: under the assumptions of the data-generating model, mean baseline encounter rates differed during daytime and nighttime, whereas under the assumptions of the restricted model these rates were assumed to be identical. Therefore, under the restricted model estimates of α₀ correspond to the logit-scale, mean baseline encounter rate averaged over daytime and nighttime periods. Apart from this difference, the operating characteristics of full and restricted models were nearly identical.

[Figure omitted. See PDF.]

Table 2. Results of simulation study.

https://doi.org/10.1371/journal.pone.0176966.t002

Discussion

In this article we described a class of hierarchical models for estimating the spatial distribution and abundance of animals detected with continuous-time recorders. In these models the spatial distribution of individual activity centers is assumed to be a realization of a spatial point process parameterized to include the effects of spatially varying covariates (e.g., habitat measurements). Similarly, the sequence of detection times observed for each unique individual is assumed to follow a temporal point process during the continuous period of each camera’s operation. This process is parameterized to include the effects of trap-specific covariates on each individual’s baseline rate of detection so that factors influencing animal movements or behaviors can potentially be inferred. The effects of behavior also can be specified using time-specific covariates, provided measurements of these covariates are available during the entire period of each trap’s operation. These novel aspects of our continuous-time SCR model provide a more realistic description of the nocturnal activities and movements of tigers and other species with similar habits.

Several components of our SCR model are also present in other SCR models. For example, the effects of spatially varying covariates on density of activity centers has been modeled in some SCR models using a multinomial distribution to specify a discrete representation of the locations expected to be occupied by each individual ([2, 21], chapter 11). A similar approach was used to specify each individual’s utilization of space with a multinomial resource selection function ([2, 22], chapter 13). This allowed [22] to establish a technical equivalence between an individual’s resource selection and encounter probability. Specifically, they showed that specifying the baseline rate of detection of an individual as a function of trap- and location-specific covariates (as in our SCR model) implies a particular resource selection function. Therefore, if resource selection occurs it can be estimated from the effects of trap- and location-specific covariates on each individual’s baseline rate of detection.

Our SCR model also shares many of the features of the continuous-time SCR model of [6]. For example, in that model the distribution of detection times was parameterized using a hazard function h_k(t, x_i; θ), wherein x_i denotes the activity center of individual i and k indexes the trap. Similarly, our continuous-time SCR model uses the intensity function ϕ(t, s_i, x_k) (defined in Eq (3)) of an inhomogenous Poisson process to specify the distribution of detection times for an individual with activity center s_i. Both models allow an individual’s baseline rate of detection to be modeled using trap-specific covariates (specified using log-linear formulations of [6] or of ψ_k (our SCR model)). The primary difference between our continuous-time SCR model and that of [6] lies in the treatment of individual activity centers. The model of [6] was used primarily to estimate model parameters, in particular the parameters that specify population density. Parameters were estimated by maximizing the marginal likelihood function obtained by integrating individual activity centers from the joint probability density of the data and the activity centers. In contrast, our SCR model was fitted using a Bayesian approach that retains the individual activity centers as formal parameters; therefore, our model can be used to estimate the spatial distribution of individual activity centers, including valid estimates of their uncertainties. As such, our SCR model essentially provides a species distribution model that allows a map of the distribution of abundance (or density) of individual activity centers to be estimated. Moreover, if spatial covariates of individual density are unknown or unavailable, our SCR model can still be used to reveal spatial variation in the distribution of individual activity centers, as was evident in our analysis of the tiger data.

Given the cost-effectiveness and widespread use of camera-trap surveys, we anticipate increasing use of continuous-time SCR models. Unlike SCR models based on discrete periods of observation, continuous-time models are more accurate representations of the actual sampling process when camera traps or other types of continous-time recorders are used. Continuous-time SCR models use all of the data, which may be advantageous in surveys of species at high density, where individuals are observed many times. An additional benefit of these models is their ability to specify aspects of the detection process more accurately when there are temporal variations in the activity of individual animals. In our analysis of the tiger data, we illustrated use of temporal covariates to model nocturnal behavior; however, semiparametric models of temporal patterns (e.g., splines or other basis functions of time) also can be used in circumstances where covariates of detection are unknown or unavailable. A challenging, but important, area of further research is the extension of continuous-time SCR models for estimating population dynamics. Such “open” models have been developed for discrete-time SCR models ([2, 23], chapter 16). Similar approaches may provide useful extensions of continuous-time SCR models. We plan to address this issue in ongoing work.

Supporting information

[Figure omitted. See PDF.]

S1 Appendix. Markov chain Monte Carlo algorithms.

https://doi.org/10.1371/journal.pone.0176966.s001

(PDF)

S2 Appendix. R source code used to fit models of tiger data observed during camera-trap survey.

This software has been approved for release by the U.S. Geological Survey (USGS). Although the software has been subjected to rigorous review, the USGS reserves the right to update the software as needed pursuant to further analysis and review. No warranty, expressed or implied, is made by the USGS or the U.S. Government as to the functionality of the software and related material nor shall the fact of release constitute any such warranty. Furthermore, the software is released on condition that neither the USGS nor the U.S. Government shall be held liable for any damages resulting from its authorized or unauthorized use.

https://doi.org/10.1371/journal.pone.0176966.s002

(ZIP)

S3 Appendix. Description of spatial covariates used in simulation study.

https://doi.org/10.1371/journal.pone.0176966.s003

(PDF)

S4 Appendix. R source code used to create simulated data and to fit models to these data.

This software has been approved for release by the U.S. Geological Survey (USGS). Although the software has been subjected to rigorous review, the USGS reserves the right to update the software as needed pursuant to further analysis and review. No warranty, expressed or implied, is made by the USGS or the U.S. Government as to the functionality of the software and related material nor shall the fact of release constitute any such warranty. Furthermore, the software is released on condition that neither the USGS nor the U.S. Government shall be held liable for any damages resulting from its authorized or unauthorized use.

https://doi.org/10.1371/journal.pone.0176966.s004

(ZIP)

Acknowledgments

We thank Ravishankar Parameshwaran, Devcharan Jathanna and Jitendra Shankaraiah for assistance with analyses. Andy Royle and an anonymous referee provided review comments that improved an earlier version of this article. We are grateful for institutional support provided by the USGS Wetland and Aquatic Research Center, the Wildlife Conservation Society—India Program, and the Centre for Wildlife Studies in Bangalore, India.

Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Author Contributions

1. Conceptualization: RMD KUK.

2. Data curation: KUK.

3. Formal analysis: RMD.

4. Funding acquisition: RMD KUK.

5. Investigation: KUK.

6. Methodology: RMD.

7. Project administration: RMD.

8. Resources: RMD KUK.

9. Software: RMD.

10. Supervision: RMD KUK.

11. Validation: RMD.

12. Visualization: RMD.

13. Writing – original draft: RMD.

14. Writing – review & editing: RMD KUK.

Citation: Dorazio RM, Karanth KU (2017) A hierarchical model for estimating the spatial distribution and abundance of animals detected by continuous-time recorders. PLoS ONE 12(5): e0176966. https://doi.org/10.1371/journal.pone.0176966

References

1. Borchers D. A non-technical overview of spatially explicit capture-recapture models. Journal of Ornithology. 2012;152(Supplement 2):S435–S444.

2. Royle JA, Chandler RB, Sollmann R, Gardner B. Spatial capture-recapture. Amsterdam: Academic Press; 2014.

3. Efford MG, Dawson DK, Borchers DL. Population density estimated from locations of individuals on a passive detector array. Ecology. 2009;90:2676–2682. pmid:19886477

4. Royle JA, Karanth KU, Gopalaswamy AM, Kumar NS. Bayesian inference in camera trapping studies for a class of spatial capture-recapture models. Ecology. 2009;90:3233–3244. pmid:19967878

5. Dawson DK, Efford MG. Bird population density estimated from acoustic signals. Journal of Applied Ecology. 2009;46:1201–1209.

6. Borchers D, Distiller G, Foster R, Harmsen B, Milazzo L. Continuous-time spatially explicit capture-recapture models, with an application to a jaguar camera-trap survey. Methods in Ecology and Evolution. 2014;5:656–665.

7. Foster RJ, Harmsen BJ. A critique of density estimation from camera-trap data. Journal of Wildlife Management. 2012;76:224–236.

8. Efford M. Density estimation in live-trapping studies. Oikos. 2004;106:598–610.

9. Borchers DL, Efford MG. Spatially explicit maximum likelihood methods for capture-recapture studies. Biometrics. 2008;64:377–385. pmid:17970815

10. Dorazio RM. Bayes and empirical Bayes estimators of abundance and density from spatial capture-recapture data. PLoS ONE. 2013;8:e84017. pmid:24386325

11. Cressie N, Wikle CK. Statistics for spatio-temporal data. Hoboken, New Jersey: John Wiley & Sons; 2011.

12. Møller J, Waagepetersen RP. Statistical inference and simulation for spatial point processes. Boca Raton: Chapman & Hall; 2004.

13. Illian J, Penttinen A, Stoyan H, Stoyan D. Statistical analysis and modelling of spatial point patterns. West Sussex, England: John Wiley & Sons; 2008.

14. Gelfand AE, Smith AFM. Sampling-based approaches to calculating marginal densities. Journal of the American Statistical Association. 1990;85:398–409.

15. Brooks S, Gelman A, Jones GL, Meng XL, editors. Handbook of Markov chain Monte Carlo. Boca Raton, Florida: Chapman & Hall / CRC; 2011. https://doi.org/10.1201/b10905

16. Karanth KU. Estimating tiger Panthera tigris populations from camera-trap data using capture-recapture models. Biological Conservation. 1995;71:333–338.

17. Karanth KU, Nichols JD. Estimation of tiger densities in India using photographic captures and recaptures. Ecology. 1998;79:2852–2862.

18. Karanth KU, Sunquist ME. Behavoural correlates of predation by tiger (Panthera tigris), leopard (Panthera pardus) and dhole (Cuon alpinus) in Nagarahole, India. Journal of Zoology. 2000;250:255–265.

19. Karanth KU, Nichols JD, Kumar NS. Estimating tiger abundance from camera trap data: field surveys and analytical issues. In: O’Connell AF, Nichols JD, Karanth KU, editors. Camera traps in animal ecology. Tokyo: Springer; 2011. p. 97–117.

20. R Core Team. R: A Language and Environment for Statistical Computing; 2016. Available from: http://www.R-project.org/.

21. Fuller AK, Sutherland CS, Royle JA, Hare MP. Estimating population density and connectivity of American mink using spatial capture-recapture. Ecological Applications. 2016;26:1125–1135. pmid:27509753

22. Royle JA, Chandler RB, Sun CC, Fuller AK. Integrating resource selection information with spatial capture-recapture. Methods in Ecology and Evolution. 2013;4:520–530.

23. Gardner B, Reppucci J, Lucherini M, Royle JA. Spatially explicit inference for open populations: estimating demographic parameters from camera-trap studies. Ecology. 2010;91:3376–3383. pmid:21141198