Estimating the effect of hunting on wild animals and their population status is a key component of wildlife management. Such estimates will help determine the appropriate level of hunting pressure for the management of distinct issues, over-exploitation (Brashares et al. 2004) and overabundance (Goodrich and Buskirk 1995). Recent advances in statistical modeling for time-series analysis and the increasing availability of long-term monitoring data can offer more precise and effective measures of population status.

Because it is notoriously difficult to directly obtain the exact abundance in the field, indirect measures of abundance (‘population indices' hereafter) are essential for evaluating temporal and spatial trends in animal populations under hunting pressure and hence for evaluating the impact of hunting on the population stock. Population indices are easily observable values that are, on average, proportional to population size. There are various population indices, such as harvest size (Cattadori et al. 2003, Ueno et al. 2014), catch per unit effort (CPUE, Roseberry and Woolf 1991), sight per unit effort (SPUE, Ericsson and Wallin 1999, Ueno et al. 2014), spotlight count (Uno et al. 2006), camera trap index (Rovero and Marshall 2009) and dung count (Goda et al. 2008), and the validity of each index has been studied extensively by examining correlations with independent population size estimates or other population indices. Harvest size (i.e. the number of individuals captured) is the most widely used statistic for game animals; however, its validity as a population index is controversial. Some studies have shown high correlations between harvest size and population size (Ranta et al. 2008, Ueno et al. 2014) and others have found biases (Cattadori et al. 2003, Imperio et al. 2010). These differences can be explained by different patterns in capture effort (i.e. the amount of effort devoted to capture animals, such as hunter-days and trap-days) among regions and species. When capture effort varies substantially both in time and in space depending on human social conditions (e.g. population dynamics of hunters and cash rewards for nuisance control) or changes in wildlife behavior in response to hunting activity (Iijima 2017), the harvest size can exhibit bias, resulting in inaccurate evaluations of population dynamics and the effect of hunting pressure, which are crucial for sustainable game management and successful population control. Nevertheless, the performance of population indices for the estimation of population parameters under different scenarios of capture effort has rarely been studied.

Recently, harvest-based estimation (HBE; Roseberry and Woolf 1991, Yamamura et al. 2008, Chee and Wintle 2010, Ueno et al. 2014, Iijima and Ueno 2016), an estimation method for population size and hunting mortality from longitudinal observations of population index (or indexes) and harvest size, has been applied in wildlife population management. Especially, HBE is a common approach to estimate sika deer Cervus nippon in Japan and used to develop management plans by the local governments (Iijima 2018). In recent studies, HBE is often formulated by generalized state-space models, which include state processes (i.e. population size dynamics determined by population growth and harvest) and observation processes of population index(es) (Iijima 2020). In particular, HBE yields successful estimates when hunting has a non-negligible impact on population size. Longitudinal data for the population index and harvest size are used to estimate model parameters.

In some applications of HBE in Japan, harvest size has been used as a population index in the observation process for HBE (Iijima 2018). However, behavior of the estimator is not well studied. Although temporal heterogeneity in capture effort is expected to result in biases in parameter estimates even if catchability (i.e. capture probability per unit capture effort) is constant over time, the robustness of estimation with respect to variation in capture effort is not clear. When capture effort is constant over time, harvest size is truly proportional to population size and there is no difference between data-generating and estimation models. However, the behavior of estimates in such an ideal situation is not clear. Studies on the behaviour of estimators based on harvest size will offer valuable insights into the limitations of harvest size as a population index to evaluate the management of exploited animal populations.

Table 1.

Capture effort scenarios. In the decreasing and increasing scenarios, capture effort decreases and increases linearly with time, respectively. Means and 50% range of Pearson's correlation coefficients between harvest size and population size are also shown.

In the present study, we conducted simulation experiments to demonstrate how HBE behaves when harvest size is used as a population index in the observation model (hereafter, proportional harvest size estimator (PHSE)). We generated a simulated dataset iteratively under different capture effort scenarios (i.e. constant, increasing and decreasing) and evaluated the identifiability and precision of estimates. To show how temporal heterogeneity in capture effort introduces estimation bias for PHSE, we also compared estimates of PHSE with those from a model that reflects the true data-generating process in which harvest size depends on capture effort.

Material and methods

Simulation settings and dataset generation

To show the effect of violations of the assumptions of PHSE, such as temporal heterogeneity in capture effort, we simulated harvest data of a hunted population under different capture effort scenarios and applied PHSE. Temporal heterogeneity in capture effort results in the overdispersion of harvest size around the average trend, which can result in biased estimates of population size and other parameters. We predicted that estimation bias would be more severe for scenarios with larger variation in capture effort.

The following seven scenarios with different temporal patterns in capture effort were used to generate datasets: a constant effort scenario (scenario 1), three decreasing effort scenarios (scenario 2–4) and three increasing effort scenarios (scenario 5–7) (Table 1, Supplementary material Appendix 1 Fig. A1). According to each scenario, the population dynamics were simulated for a 20-year period. In scenario 1, capture effort in every year was 100 hunter-days and invariant with time. In decreasing and increasing effort scenarios, capture effort changed linearly with year. Three levels of change (low, moderate and high) were established for decreasing and increasing effort scenarios. The average capture effort across years was 100 hunter-days for all scenarios.

Although stochastic structures of the underlying model for HBE vary among studies (Yamamura et al. 2008, Chee and Wintle 2010, Fukasawa et al. 2013, Iijima et al. 2013, Osada et al. 2015, Iijima and Ueno 2016), they reflect population growth and artificial removal from the population. In this study, the Poisson-binomial surplus-production model (Chee and Wintle 2010, Fukasawa et al. 2013) was applied because it explicitly incorporates a process of stochastic population growth and the removal of individuals from a finite and countable population. For simplicity, we assumed exponential population growth. Let N_t denote the population size at the beginning of year t, the population size after population growth, S_t, is assumed to follow a Poisson distribution as follows:

where exp(r) is the annual population growth rate. Harvest size, C_t, is determined by a binomial process with capture probability p_t and population size at year t:

and the population size in year t+1 is determined by removing harvest size from S_t:

To model the relationship between capture probability p_t and capture effort, the Poisson catchability model (Seber 1982, Skalski et al. 2005), which assumes random encounters of individuals and unit capture effort, was applied as follows:

where exp(φ) is the catchability coefficient (i.e. capture probability per unit effort) and E_t is the capture effort (i.e. hunter-days). Equation 4 corresponds to an inverse complementary log–log link function that constrains p to between 0 and 1 for any given real number φ.

For the seven scenarios, virtual harvest sizes were generated by Monte Carlo simulations following the Poisson-binomial surplus-production model. True r, φ and N₁ were ln(1.21), –6 and 1000, respectively. The population growth rate, exp(r), was determined with reference to the sika deer population growth rate in Shiretoko peninsula, Japan (Kaji et al. 2004). The other parameters, exp(φ) and N₁, were determined so that 1) the population never goes extinct, 2) population dynamics are clearly affected by hunting pressure and 3) computational time for parameter estimation is realistic. The order of computational time is cubic of the maximal possible population size estimate (for details, Supplementary material Appendix 2), and setting initial population size too large can result in unrealistic computational time in parameter estimation. Monte Carlo simulations were run for 100 iterations. The number of iterations was determined to finish our calculation within a realistic computational time.

Specifications of estimation models

Prior to defining the estimation model for PHSE, we considered the statistical model for HBE which has the same structure as the data-generating model. In terms of statistical modelling, the Poisson-binomial surplus-production model can be regarded as a discrete-valued generalized state–space model. The generalized state–space model is a statistical modeling framework for time-series analyses for the estimation of system parameters and latent states, considering two sources of variability: observation error and process variability. It can be written as a set of three probability distributional functions, the initial state distribution, state process distribution and observation process distribution (Buckland et al. 2004). The state–space model which has the same structure as the data-generating model is described as follows:

Initial state distribution:

State process distribution:

Observation process distribution:

where the vector of parameters θ=(r, φ).

PHSE is based on the assumption that the expectation of harvest size is proportional to population size, [C_t] ∝ N_t; it is equivalent to assume that the capture probability p_t is constant. The underlying PHSE model can be derived by a slight modification of the state–space Poisson-binomial surplus-production model by replacing Eq. 7 with the following:

Because Eq. 8 does not include E_t in contrast to Eq. 7, exp(φ) is no longer interpreted as a catchability coefficient but as a transformed capture probability. Note that the Poisson-binomial surplus-production model for scenario 1 has an identical structure to the model for PHSE because true capture effort is invariant over time, and they are expected to return the same population size estimates. Under this model, [C_t] = exp(r)pN_t, satisfying the condition that the harvest size is, on average, proportional to population size.

Results

Population size generated by our simulation model showed different nonlinear temporal trajectories among scenarios. Under constant effort, both harvest size and population size showed exponential decay with time, on average (Supplementary material Appendix 1 Fig. A2a, A3a). When effort decreased with time and the rate of decrease was high, the harvest size showed temporal decay in which the harvest size growth rate was temporally heterogeneous but population size exhibited a U-shaped trajectory (Supplementary material Appendix 1 Fig. A2b–d, A3b–d). Harvest size and population size were unimodal when effort increased with time (Supplementary material Appendix 1 Fig. A2b–d, A3b–d). Pearson's correlation coefficient for the relationship between harvest size and population size was high (> 0.9), except in scenario 6 and 7 (Table 1).

The sensitivity of PHSE to a change in capture effort was high. Initial population size and annual growth rate estimates of PHSE exhibited downward and upward biases from true values when the true capture effort changed with time, respectively (Table 2). Estimation bias was larger when the rate of change in effort was larger and was more severe in increasing effort scenarios than in decreasing effort scenarios. Estimated initial population sizes for scenario 6 and 7 were 51.0 and 7.45, on average, respectively. Pearson's correlation coefficients of estimated mean population size and true population size were also low for these scenarios (Table 2). The 50% range of population size estimates did not cover the true value, except in scenario 1 and 2. Figure 1 shows the ratio of estimated and true population sizes under PHSE against Pearson's correlation coefficients of harvest size and true population size for each scenario. Even when correlation coefficients exceeded 0.9, there were scenarios in which initial population size estimates by PHSE were largely biased. The original Poisson-binomial surplus-production model did not show such a downward bias and correlation between estimated and true population size was substantially high (Table 3). Although a weak upward bias was found, on average, the 50% range of initial population size estimates captured the true value for all scenarios. The original Poisson-binomial surplus-production model yielded more precise estimates with lower variance when the change in capture effort was larger.

Table 2.

Summary of estimates of PHSE applied to 100 simulated datasets. Means and 50% ranges were derived from estimates of identifiable iterations. Means and 50% range of Pearson's correlation coefficients between harvest size and population size are also shown.

The ratio of identifiable iterations was low when there was little or no temporal variability in capture effort for both the PHSE and original Poisson-binomial surplus-production model (Table 2, 3). In particular, it was 10% when capture effort was constant, the case in which the data-generating model and estimation model for PHSE are identical.

Discussion

Although harvest size is often highly correlated with population size and is thought to be a viable population index (Ranta et al. 2008, Ueno et al. 2014), our results indicate that caution is needed when it is used to estimate population size and population parameters by HBE. With low variation in capture effort, population size and population growth rate estimates exhibited serious bias. Even when the correlation between population size and harvest size was high (e.g. Pearson's correlation coefficient > 0.9), there were cases in which population size was underestimated by about half, on average. Another drawback of PHSE is the low identifiability in the ideal situation in which capture effort is constant and the estimation model is identical to the data-generating model. In the constant effort scenario, a model of PHSE is identical to an original Poisson-binomial surplus-production model. Therefore, the low identifiability under constant effort is an intrinsic constraint of the Poisson-binomial surplus-production model. These results imply that longitudinal data for harvest size only would not offer sufficient information to evaluate the effect of harvest on wildlife populations and PHSE could not compensate for the lack of information of data itself.

Figure 1.

Divergence between the estimated initial population size and the true value when harvest size is treated as the population index in relation to the Pearson's correlation coefficient for the relationship between harvest size and true population size under different data-generating scenarios (1–7), as described in Table 1. Dots and bars indicate means and 50% ranges of identifiable iterations, respectively. (A) full-scale and (B) extended figures when the correlation coefficient is larger than 0.9 are shown.

The underestimation of population size by PHSE can be explained by the relationship between model parameters, the resulting harvest size trend and dispersion around the average trend. Yamamura (2016) showed that population growth rate and population size are not identifiable from the growth rate of the population index. In a similar way, we can derive the following equations to explain the relationship between parameters and harvest size trends under the assumption of the Poisson-binomial surplus-production model from Eq. 1 to 4:

If we set up a system of nonlinear simultaneous Eq. 9 over different ts, there is a unique solution for the log population growth rate r and log capture efficiency φ, but special combinations of capture efforts E_t at different times, making the system singular. However, the PHSE assumption of constant capture effort is a special case in which simultaneous equations do not have a solution: and . Therefore, the identifiability of model parameters of PHSE depends on the dispersion component of harvest size, rather than the average trend. The Poisson-binomial structure of the surplus-production model is the same as that of the N-mixture model, which is well-studied owing to its capacity to estimate abundance for unmarked populations using information from dispersion among repeated counts (Royle 2004, Kéry 2018). Population size estimates by the N-mixture model are highly sensitive to the distributional assumption of the model (Duarte et al. 2018, Knape et al. 2018), and the Poisson-binomial surplus-production model without information for different capture efforts would share the same feature. In the case of PHSE, overdispersion of observed harvest size would increase the estimate of p because removing a large proportion of individuals induces large fluctuations in population size and hence fluctuations in harvest size. Under a ‘false’ assumption of constant capture effort, the overdispersion of harvest size induced by temporal heterogeneity in capture effort can be misinterpreted as a high capture probability, large population growth rate and small initial population size. Then, the overestimation of p and r and underestimation of N₁ would occur.

Table 3.

Summary of estimation of the Poisson-binomial surplus-production model applied to 100 simulated datasets. Means and 50% ranges were derived from estimates of identifiable iterations.

Low identifiability for a constant actual capture effort is also due to the dispersion dependence of parameter estimates. As mentioned above, the system of Eq. 9 is singular when capture effort is constant and it can be difficult to determine the model parameters from the dispersion of harvest size. Considering that estimates of the original Poisson-binomial surplus-production model were more precise when the change in capture effort was large, the change in capture effort and the harvest size response would be important information to ensure parameter identifiability. Our results further indicate that harvest-estimation should be applied when capture effort varies with time and it can be explicitly modelled. If capture effort has low variation, we recommend the inclusion of auxiliary information, such as independent estimates of population size (Fukasawa et al. 2013, Iijima et al. 2013) and a reliable estimate of intrinsic growth rate (Yamamura et al. 2008, Yamamura 2016), to improve model identifiability and stability.

Although we considered single time-series and assessed the effects of ignorance of temporal heterogeneity in capture effort on estimation bias in our study, heterogeneity in space would have similar effect on estimation. HBE can be extended to accommodate multiple spatial units with different population sizes, capture efforts and harvest sizes (Iijima et al. 2013, Osada et al. 2018). In such cases, ignoring spatial heterogeneity in capture effort as well as temporal heterogeneity can induce overdispersion of harvest size hence estimation bias of model parameters. However, spatial heterogeneity in capture effort would improve the model identifiability if capture effort were explicitly modelled.

In this study, we applied marginal maximum likelihood estimation via recursive filtering to HBE and easily and clearly determined model identifiability. Although Gibbs sampling is often used for HBE (Chee and Wintle 2010, Fukasawa et al. 2013) and data cloning (Lele et al. 2007) is an available option for the evaluation of model identifiability within a framework of Gibbs sampling, the computational cost is high, which is inappropriate for our study in which state–space models were applied to large quantities of simulated data. For the application to actual datasets, marginal maximum likelihood estimation presented in this study is useful to check model identifiability.

In conclusion, we recommend wildlife managers monitor both harvest size and capture effort to make sound estimation of population size and harvest rate. The use of harvest size as a population index in HBE can lead to seriously biased estimates of population size and growth rate when capture effort shows a temporal trend. Additionally, independent survey of population size or population growth rate would be needed to overcome low identifiability of parameters if temporal change in capture effort is small. Obviously, such problems would be shared by other population indices that are closely correlated with harvest size. Generally, the management of large mammals requires reliable estimates of the population status because maintaining their density at a socially acceptable level requires fine-tuned control. In particular, wildlife managers in many regions must address the overabundance of large mammals and the urgent need for population control (Fagerstone and Clay 1997, Côté et al. 2004, Uno et al. 2009). The underestimation of population size leads to the failure of population control projects via the establishment of inadequate quantitative targets for removal. We hope our research encourages managers to collect unbiased population indices, not only harvest size.