Capture-recapture (CR) model for the closed population

First, I treat the estimation way of abundance in the closed population by the CR survey because the CR model for the closed population is considered to be the most basic of all abundance estimation models. “Closed” means that individuals in the population do not increase or decrease by birth and death, and there is no migration during the study period. If individuals in the population can be identified and detected more than once during several sampling occasions, their abundance can be estimated as below (Williams et al. 2002):

is the estimator of abundance, C is the number of counted wildlife, and is the estimator of the detection probability. Although the name of the method is “capture-recapture”, researchers do not need to capture wildlife physically if they can distinguish individuals from any patterns or characteristics of their body. Then, the CR method is also known as the “mark-resight” (McClintock et al. 2009) and “mark-recapture” (Greenwood et al. 1985). At the same time, “detection” is also the same meaning as “capture.”

Let's imagine the case of two sampling occasions in estimating and under the constant detection probability. n₁ is the number of animals caught and marked on the first occasion, n₂ is the number of animals caught on the second occasion, m₂ is the number of animals caught on both occasions, and N is the true population size. If all animals have equal capture probabilities, the proportion of marked animals in a population after the first sample should approximate the proportion of marked animals in the second sample, i.e.,

Rearrangement of terms in equation (2) leads to the estimator of N () as below, and is known as the Lincoln-Petersen estimator (Petersen 1896; Lincoln 1930):

In this case, is derived as m₂/n₂.

It is useful to generalize the procedure to obtain the estimators by using probability distribution. The multi-nomial distribution is as appropriate as the model in estimating abundance and detection probability. In the case of the above example (i.e., two sampling occasions with a constant detection probability), the probability for the data can be expressed below:

The above probability can be extended to K times the sampling occasions as below:

x_ω is the number of caught wildlife in each capture history (ω) and M_K+1 is the total number of marked individuals caught during the study. Specifically, in the case of three sampling occasions, ω = {(0,0,1), (0,1,0), (0,0,1), …, (1,1,1)}, where the “0” and “1” indicate non-capture and capture in each sampling occasion, respectively.

Recently, the formalization of the CR model for the closed population as a hierarchical model was conducted by Royle et al. (2007). The hierarchical model along with the ecological model describes the ecological and latent process of the targeted wildlife, and the observation model describes the process to observe the latent state of the targeted wildlife (Royle and Dorazio 2008). In the hierarchical CR model for the closed population, the true but latent state about the presence/absence of each individual is estimated. Then, the detection/non-detection of each individual in the population is treated in the hierarchical CR model for the closed population. Implementation of the hierarchical CR model for the closed population is accomplished by the “parameter-expanded data augmentation (PX-DA)” technique (Royle et al. 2007). PX-DA consisted of two things by Royle et al. (2007): (1) adding an arbitrary number of zeros to the data set and (2) analyzing a reparametrized version of the original model. Let's assume that n is the number of marked individuals, N is the latent population abundance, and M is the augmented data size. The relationship among them is n < N < M and N = ΨM: Ψ is the inclusion probability. Then, the number of zeros that are added to the original data should be larger than N. The zeros can be divided into two components, i.e., “real” but unmarked individuals and individuals that do not exist in practice. The sum of n and the “real” but unmarked individuals is N. Hierarchical CR model for closed population by using data augmentation is expressed below:

Ψ is the inclusion probability as stated above, p is the detection probability, z_i is the latent state of ith individual whether ith individual truly exists in the target population (“1”) or not (“0”), and y_i,k is the augmented data of detection (1) or the non-detection (0) of the ith individual in the kth sampling occasion.

The detection probability can differ during the survey if you obtain the appropriate covariate(s) of the detection probability. Otis et al. (1978) organized various models for detection probability. M₀ is the model with the constant detection probability (as above example), M_t is the model with the time-dependent detection probability, M_b is the model with the behavior-dependent detection probability, and M_h is the model with the individual-dependent detection probability. M_t is suitable when the detection probability changes with time. For example, if the detection probability increases in winter because of the food shortage, M_t should be used. M_b is suitable when the detection probability changes by the behavioral response to trapping. For example, if the detection probability differs before and after the first capture of wildlife, M_b should be used. The change of the detection probability by trapping is a well-known phenomenon that is known as “trap happy” and “trap shy” (Tuyttens et al. 1999; Zwolak and Foresman 2008; Augustine et al. 2014). If there is a behavioral response to trapping, the duration time of the behavioral response should be considered.

Because of too many parameters of M_h, the parameters of individual effect cannot be identifiable unless each of the parameters is regarded as a random effect that follows any probability distributions (Otis et al. 1978). Then, Huggins (1989) and Alho (1990) reorganized the models as fully-observable covariate model (M_t, M_b) and individual covariate models. Although the difference between fully-observable covariate models and individual covariate model seems difficult to understand, the former treats the effect of factor (i.e., time or the catch history) as a fixed effect and the latter treats the effect of individual heterogeneity as a random effect. As stated above, we cannot essentially observe “individual heterogeneity”. Then, it is reasonable to consider the effect of individual heterogeneity as a random variable. Because all individuals including unobserved individuals are treated in the hierarchical CR model as we see in equations (7) to (11), the effect of individual heterogeneity on detection probability can be easily treated in the hierarchical CR model.

Spatially explicit capture-recapture (SECR) model

CR models basically cannot estimate the density of wildlife because the effective sampling area for estimation cannot be defined in the model. To estimate the density of wildlife, the location of individuals in an arbitrary space should be defined. Although there are many studies about SECR models (Hartstack et al. 1971; Efford 2004; Borchers and Efford 2008), I introduce the idea of Royle and Young (2008). The treatment of the activity center of each individual enables the estimation of the effective sampling area. Royle and Young (2008) modeled the detection probability as the relationship between the activity center of each individual and the location of observation.

Here I explain the SECR model to estimate the wildlife abundance and their activity centers for the closed population from the camera-trapping survey data. J camera traps are set in a square space and the photographed wildlife can be identified. The locations of the J camera traps with coordinates are c_j = (c_1j, c_2j); j = 1, 2, ..., J. n individuals are totally identified, and the number of photos of the ith individual of the jth camera in the kth sampling occasions (y_{i, j,k}) is recorded. However, some individuals may not be detected during the study period. In order to estimate true abundance N in an arbitrary space region, data augmentation is useful in estimating N as explained in the CR model for the closed population. An arbitrary large number of individuals with all zero records are added to the data, and N is modeled as equations (7), (8), and (9).

Next, I introduce another square space S that encloses the region where J camera traps are set. The activity center of each N individual in S such as s_i = (s_1i, s_2i); i = 1, 2, ..., N, is modeled as below:

In this case, N is the estimate of abundance, the area of S is the effective sampling area, and N/S is the estimate of density.

With the location of activity centers of all individuals, the detection probability is modeled as the distance between the activity center of each individual (i.e., s_i) and the location of observation (i.e., c_j). Although many functions can be used, the half-normal function is a popular function. By using the half-normal function, the detection probability can be modeled as below:

p_i,j is the detection probability of the ith individual by the jth observation point, is the squared distance between the ith individual's activity center and the jth observation point (i.e., (s_i – c_j)²), and σ² is the scale parameter of the half-normal function. Therefore, the number of photos of each individual in each sampling occasion can be expressed as below:

λ₀ is the baseline detection intensity.

Here we assume that the population is closed. However, the SECR model can be extended to the open population (Gardner et al. 2010). Furthermore, Chandler and Royle (2013) extended the SECR to be able to be applied to the population without the individual identification or partially identified population (unmark SECR).

Capture-recapture models for open population

The estimation of the population abundance by the CR model for the closed population is robust because the relationship between the true abundance and observed counts is obviously expressed by the model. However, the assumption of the population closure is sometimes difficult to accomplish in ecological studies. Furthermore, the temporal change of abundance itself can be the interest of research. Therefore, the CR model for the closed population has been extended to the open population. “Open” means that the individuals in the population can increase or decrease by their natural mortality, birth, and migrations during the study period. To estimate the abundance of open population, the survival and entering of individuals for the target population, in addition to the detection of wildlife, should be treated.

The Jolly-Seber (JS) model (Jolly 1965; Seber 1965) estimates the entry probability in addition to the survival and detection probability. Although there are some variations of the JS model, I explain the super-population approach (Williams et al. 2002) here. For simplicity, I explain the JS model of super-population approach as the hierarchical model. In the super-population approach, we assume that the individuals enter the target population from a “super-population” by each sampling occasion. The entry probability of the individuals from the super-population is defined as below:

b_k is the entry probability of the kth survey, a = (α₁, α₂, …, α_k) is the parameter of Dirichlet distribution, and η_k is the conditional entry probability of the kth survey. Usually, each component of a is set as 1 (a non-informative prior distribution).

By using the entry probability in addition to the survival probability, we can model the survival and entry of each individual in each sampling occasion below. In this case, we assume the constant survival probability:

φ is the survival probability and z_i,k is the latent variable of the presence or absence of the ith individual in the kth sampling occasion in the target population. Equation (19) indicates the presence or absence of the ith individual in the k+1th sampling occasion (z_i,k+1) is determined by the survival of the ith individual (φ) that existed in the target population in the previous sampling occasion (z_i,k) and the products of absence probability until the kth sampling occasion and the entry probability in the k+1th sampling occasion (η₊₁).

Next, I explain the observation model. Since the detection of the individual is generally imperfect in the field survey, the detection probability should be estimated. In this case, I assume that there are no covariates among the detection probability as below:

Because the detection is imperfect, the unmarked but existing individuals in the target population should be considered to estimate the population abundance. As explained in the closed population, the data augmentation technique is again useful. The estimation of the latent state of M individuals can be modeled as below:

Ψ is the inclusion probability as stated above, and w_i is the latent variable among the presence (1) or absence (0) of the ith individual in the target population during the whole survey.

By using the above parameters, the observed data including the added zeros can be modeled as below.

y_i,k is the detection (1) or non-detection (0) of the ith individual in the kth sampling occasion. Furthermore, abundance of the kth sampling occasion (N_k = (N₁, N₂, ..., N_K)) and total abundance of the targeted population (N) can be derived as below:

Nind_i is the number how many times the ith individual is detected, and Nalive_i denotes its state (alive or dead, or included or not-included). These two derived parameters are necessary to calculate N.

It is possible that the survival and detection probabilities fluctuate with covariates like time, individual, and capture history, as is case of the detection probability of the CR model for the closed population. However, in this case, there are two parameters. It is necessary to measure covariates in each sampling occasion and to model the survival or detection probability with the measured covariates. However, if there are correlations among the factors of fluctuation of survival or detection probability, parameter identifiability may be lost. Therefore, I recommend researchers to adopt the “robust design” (Pollock et al. 2002). The central philosophy of the robust design is to conduct multiple surveys (i.e., replication) under a constant condition. Under the robust design, the sampling occasions are divided into two hierarchies. As an example of the robust design, let's assume the situation in which you would like to estimate the effects of the climate on the survival and detection probabilities in the count survey of some wildlife species. In the first hierarchy, factors like climate can change in-between the sampling occasions. In the second hierarchy, more than two samplings should be conducted within a short period when the factors can be regarded as stable. Because it is reasonable to assume that there is no mortality among samplings within a second hierarchy, the variation of the number of counts in the second hierarchy reflects the detection probability. Therefore, the usage of the robust design is important in guaranteeing the identifiability of parameters in a model.

Background

All models that are explained above coerce the researchers to identify individuals in the population. However, the identification of wildlife is sometimes impossible or practically difficult. The difficulty to identify wildlife has motivated the development of population abundance estimation way without identification of the wildlife. From here, I review such estimation models.

Distance sampling

Distance sampling is to count the number of individuals from fixed points or fixed transects. There are many ways to conduct distance sampling, but here I introduce the distance sampling by using the line transect sampling and the abundance estimation method by following Buckland et al. (2001). The basic model for the estimation of the wildlife density by distance sampling is expressed below:

is the estimator of wildlife density, n is the total number of counted wildlife from a line of L length within w distance from the line, and is the probability of detection for an object within an area a. can be expressed as below:

g(x) is detection function in relation to x length which is the distance between the line and the object. For simplicity, is sometimes expressed as µ and equation (28) becomes as below:

is the estimator of an effective strip half-width (ESW) that is the distance from the line for which as many objects are detected beyond µ as are missed within µ of the line. Various probability functions are available for g(x), and the selection of probability function can be done by the likelihood ratio test or AIC.

The assumptions of distance sampling by the line transect sampling are categorized into experimental design and observation. The assumptions of the experimental design are that 1) objects are spatially distributed in the area to be sampled according to some stochastic process with rate parameter D (i.e., number per unit area) and 2) randomly placed lines or points are surveyed, and a sample of n objects are detected, measured, and recorded. The assumptions of observation are that 1) objects directly on the line or point are always detected (i.e., they are detected with probability 1); 2) objects are detected at their initial location, prior to any movement in response to the observer; and 3) distance (and angles where relevant) are measured accurately (ungrouped data) or objects are correctly counted in the proper distance interval (grouped data).

Recently, the formulation of the distance sampling model as a hierarchical model (hierarchical distance sampling model) was accomplished by Kéry and Royle (2015). Although the details of the hierarchical distance sampling model are not introduced here, the hierarchical distance sampling model is especially useful in modeling the effects of the individual covariates and the spatial heterogeneity of abundance.

N-mixture model

The N-mixture model was proposed by Royle (2004) to estimate the local and latent abundance of wildlife and detection probability of the latent abundance by repeated counting in space and time. The N-mixture model is expressed below:

λ is the mean abundance of the target population, N_i is the local abundance in the ith location, p is the detection probability, and y_i,k is the observed abundance in the ith location and the kth sampling occasion.

To apply the N-mixture model, researchers need to conduct repeated observations in space and time under closed conditions. The assumptions of the N-mixture model are as follows:

Link et al. (2018) examined the biases of the N-mixture model under the three violations of model assumptions and recommended to collect reliable data about detection probability by CR methods because the violation of constant detection probability especially caused the biased estimates.

While, some family models of N-mixture model can estimate the abundance even if the duplicated count of the same individual occurs in a single sampling occasion. Nakashima (2020) showed that the Royle-Nichols model and Poisson/Poisson N-mixture model provided reasonable estimates of the number of animals in such situation. The Royle-Nichols model (Royle and Nichols 2003) treats the duplicated counts by reducing the information of data from abundance to presence/absence. The formulation of the Royle-Nichols model is as below:

In this case, y_i,k is the presence/absence of target wildlife in the ith location and the kth sampling occasion. However, it should be noted that the Royle-Nichols model can estimate realistic abundance only when abundance and detection are low (Kéry and Royle 2015). The Poisson/Poisson N-mixture model (Kéry and Royle 2015) treats the duplicated counts by using Poisson distribution as the alternative of binomial distribution of the observation model. The formulation of the Poisson/Poisson N-mixture model is as below:

In this case, y_i,k is the counts of target wildlife in the ith location and the kth sampling occasion and ϕ_ik is a rate parameter in the ith location and the kth sampling occasion.

Random encounter (RE) model

Rowcliffe et al. (2008) applied the model about rates of the collision between gas molecules into density estimation of wildlife by using camera-trapping data and named as the RE model. The RE model for wildlife density estimation is expressed as below:

D is the density of wildlife, y is the observed counts of wildlife, t is the time of sampling, π is the ratio of the circumference of a circle to its diameter, v is the animal speed of movement, r is the radius of detection, and θ is the detective angle from the observer.

Assumptions of the RE model are as follows: (i) animals conform adequately to the model used to describe the detection process; (ii) photographs represent independent contacts between animal and camera; and (iii) the population is closed. Furthermore, as stated above, researchers need to clarify the movement speed of the target wildlife. The movement speed of target wildlife can be obtained by radio-tracking of target wildlife, but it is time-consuming to attach many radio-collars or GPS-collars to wildlife.

Random encounter and staying time (REST) model

Nakashima et al. (2018) extended the RE model as researchers can estimate wildlife density by only the data from camera trap. They changed the RE model as below:

D is the population density of target animal, E(Y) is the expected number of encounters, E(T) is the expected staying time, s is the area of the detection zone, and H is a research period. The important point is to use video mode instead of camera mode to evaluate the animal speed of movement. For evaluation of the animal speed, researchers should set signs in front of the camera trap to clarify when an animal will enter and how long the animal will stay in the detective zone in front of the camera.

There are seven assumptions to apply the REST model. (1) Camera traps must be placed randomly with respect to the spatial distribution of animals. It is the issue of experimental design. (2) Cameras must certainly detect the animals entering the focal area throughout the research period. It is the issue of performance of cameras. (3) Animal density must not vary during the research period; (4) animal movement and behavior are not affected by cameras; (5) observations are independent events. These are the issue of animal behavior. (6) The observed distribution of staying time in the focal area must represent a good fit for the distribution that animal movements actually follow; (7) the observed staying time must follow a given parametric distribution. These are the issue of model structure. Nakashima et al. (2018) examined the violation of some assumptions and showed that the REST model provided unbiased estimates of animal density even when animal movement speeds varied among individuals, and when animals traveled in pairs. However, the REST model is vulnerable to unsynchro-nized activity patterns among individuals.

Removal sampling

Removal sampling is to remove individuals from the target population repeatedly. Although the name is “removal” sampling, the sampled individuals are not necessarily removed eternally from the target population. For example, the removed individuals can be returned to the target population after the survey. Otherwise, an individual can be regarded as “removed” if an observer observes the individual and intentionally never counts the individual again (a mental removal protocol, Kéry and Royle 2015). Under the assumption of closed population and constant detection (catch) probability, the relationship among the number of removed individuals in K times sampling occasions (X_k = (X₁, X₂, …, X_K)), abundance (N_k = (N₁, N₂, …, N_K)), and catch probability (p) can be modeled as below (Dorazio et al. 2005; Bord et al. 2014):

N_k is the abundance in kth sampling occasion, M is the arbitrary large number, and p is the catch probability. In addition to the above assumptions, the removal should reduce the true abundance marginally. In other words, removal sampling cannot be applied when the catch probability is very low or very high. For example, if the removal is not effective to reduce wildlife and then the catch probability is very low, the number of removed individuals may be similar regardless of sampling occasions. The abundance and catch probability can be estimated even in such a situation if researchers conduct many, many samplings and the population closure is certainly satisfied during the survey period, but it is practically difficult to conduct such a sampling especially in the field.

The catch probability per unit effort may decrease with the increase of sampling occasions (Iijima 2017) because individuals with low vigilance (bold individuals) tend to be sampled earlier than those with high vigilance (Honda and Iijima 2016). Mäntyniemi et al. (2005) developed the heterogeneous catchability model, which assumes that the catch probability differs among individuals in the target population and individuals with high catch probability are caught at earlier sampling occasions. The model is expressed as below:

X_k is the number of caught individuals in the kth sampling occasion, N_k is the abundance in kth sampling occasion, qk is the catch probability in kth sampling occasion, and µ is the mean catch probability of individuals in the target population. Because q_k decrease with the sampling occasions from equation (44), the mean catchability declines between consecutive removals but the decrease rate depends on the variation of the catchability in the initial population.

Harvest-based model

In recent years, the increase of some wildlife (e.g., sika deer Cervus nippon; Iijima et al. 2013) or the expansion of exotic species (e.g., the small Indian mongoose Herpestes auropunctatus; Fukasawa et al. 2013) causes the necessity of population management. In the population management, it is necessary to estimate population abundance and the effect of harvest on population abundance. In such situation, the population management rarely finishes within a short period (i.e., closed population). Then, the model to estimate wildlife abundance of open population under harvesting pressure is necessary. Because the sampled individuals are not necessarily removed eternally from the target population and the target population should be closed in the removal sampling model, another model should be defined. Here, I define “the harvest-based model” as the model that estimates wildlife abundance of open population under harvesting pressure and has the ecological process model that explicitly describes the decrease of abundance by the number of harvested wildlife.

Especially in targeting open population, sampling effort tends to differ spatially and temporally. Then, constant catch probability cannot be satisfied. Furthermore, demographic rates should be estimated because the population size can change between any two sampling occasions by birth, natural mortality, or migration. The basic model structure of the harvest-based model is as below:

r is the population growth rate and p_k is the catch probability in kth sampling occasion. As shown in equations (49) and (50), we need to estimate r and pk.

The Weibull catch-effort model (Barron et al. 2011; Fukasawa et al. 2013) treats the effect of sampling effort and its change with the increase of sampling occasions on catch probability as below:

p_k is the detection probability in kth sampling occasion, a is the coefficient of catchability, N_k is the abundance in kth sampling occasion, E_k is the effort in kth sampling occasion, and β is the shape parameter of Weibull distribution. The idea of the Weibull catch-effort model is that the detection probability increases with sampling effort but the degree of detection probability increase with sampling effort differs as the difference of parameter β.

To estimate demographic rates, it is useful to collect an abundance index in each sampling occasion. The abundance index is any data that is expected to correlate with the true abundance. However, because the harvest-based model explicitly treats the number of harvested individuals to reduce abundance in latent ecological process and the number of harvested individuals is affected by hunting effort and hunting efficiency (Iijima 2017), the number of harvest individuals must not be used as abundance index in the harvest-based model although the number of harvested individuals are known to correlate with abundance (e.g., Ueno et al. 2014).

Let's imagine the case of the three-year survey in estimating the wildlife abundance and demographic rate. The number of harvested wildlife of each year (sampling occasion) is known as X_k = (X₁, X₂), and the wildlife abundance index in each sampling occasion is obtained as I_k = (I₁, I₂, I₃). As the underlying process of wildlife, it is assumed that the wildlife abundance decreases by hunting (X_k) and increases by a constant annual population growth rate (r). If the true abundance of wildlife is assumed as N_k = (N₁, N₂, N₃) and the initial abundance index is assumed as I₁ = 1, population dynamics can be expressed as below:

The abundance index of the second and third sampling occasions can be written as below:

By assigning equation (54) to equation (52) and equation (55) to equation (53), r can be obtained as below:

Because we can obtain X₁, X₂, I₂, and I₃ as data, we can estimate r. With r, we can also estimate N_k by assigning equation (56) into equations (52) and (53). The parameter identifiability will greatly improve if the abundance index is the absolute abundance or density. Then, the combination of abundance estimation at certain period by models for the closed population with the harvest-based model is effective to increase parameter identifiability. Such a model can be regarded as an integrated population model (IPM).

The monitoring of covariates that affect the demographic rate is also useful in estimating demographic rates. Sæther et al. (2008) and Iijima et al. (2013) used environmental covariates like rainfall and the type of landscape to model temporal and spatial differences of the annual population growth rate of wildlife. Iijima and Ueno (2016) used the percentages of deciduous forests, evergreen forests, and artificial grasslands to model spatial variation of carrying capacity of sika deer. Use of covariates contributes not only to increase parameter identifiability but also to strengthen the biological meaning of results.