Estimating harvest when hunting bag data are reported by area rather than individual hunters: A Bayesian autoregressive approach

Harvest estimation is a central part of adaptive management of wildlife. In the absence of complete reporting, statishods are required to extrapolate from the available data. We developed a Hierarchical Bayesian framework tailored for partial reporting where hunting areas covered by reporting hunting teams are available. The framework accounts for autocorrelation at the national, county, and hunting management precinct levels. We derived and evaluated an approximation for the probability of harvest on non-reported areas under a non-linear relationship between harvest area per team and harvest rate. We applied the framework to reports of red fox ( Vulpes vulpes ), wild boar ( Sus scrofa ), common eider ( Somateria mollissima ), and grey partridge ( Perdix perdix ) harvest in Sweden from the hunting years 1997/1998 – 2019/2020. The approximation was evaluated and determined to be sufficiently accurate. We showed that accounting for autocorrelation in harvest reduced uncertainty and increased predictive accuracy, particularly for game hunted in low numbers and variably between teams. The analysis also revealed that hunting rate has a sub-linear relationship with a team ’ s area for all considered species. Further, the framework revealed substantial differences across regions in terms of parameters modeling the distribution of huntable land across teams as well as parameters modeling harvest rates. We conclude that the framework is useful to detect shifts in hunting rates and/or practices.


Introduction
Estimation of game harvest is essential for informed management (Elmberg et al., 2006;Apollonio et al., 2017;Aebischer, 2019).Game management must be adaptive, which requires adjustments to changes in ecosystems and society as well as development and application of new knowledge (Walters, 1986).To sustain adaptive management systems, access to relevant data of sufficient quality and spatial/temporal resolution is vital.Harvest statistics are a basic and important component in most game management programs, and, for some species, it may be the only data available (Cretois et al., 2020).Systems for data collection as well as data availability may vary among countries and game species, such as voluntary versus mandatory reporting (Aubry et al., 2020).
In the absence of complete reporting, statistical methods are required to extrapolate from available data to obtain estimates of the total harvest.The approach must be tailored to the data structure and data type, and several statistical methods have been developed to address issues associated with different types of surveys (Aubry and Guillemain, 2019;Aubry et al., 2020).These typically target a subset of a known population of hunters from which the harvest of the entire population is estimated.
The Swedish system, however, faces a different issue.Hunting teams rather than individual hunters report the team's harvest, and there is no information available about the total number of hunting teams or how many hunters belong to each team (Fig. 1).Thus, there is neither information about the number of reporting hunters nor the non-reporting population for which the bag must be estimated.However, the reports include, in ha, the area covered by each team, and the total area of huntable land has been estimated from land-cover maps (Jonsson et al., 2020).Consequently, estimation of the unreported bag must be performed by estimating the bag of the unreported area based on the reported area.
Analysis of longitudinal reporting data faces additional challenges.Ecological time-series typically exhibit temporal autocorrelation (Pimm and Redfearn, 1988;Lindström et al., 2012), which can arise through environmental drivers that are themselves typically positively autocorrelated (Vasseur and Yodzis, 2004;García-Carreras and Reuman, 2011) or through intrinsic factors (Blarer and Doebeli, 1996;White et al., 1996); rarely is the subsequent population size independent of the current.The population of hunters and their interest in various game species are also unlikely to change rapidly over time.Autocorrelated time series can pose a problem for statistical analyses, and tests that disregard dependence over time can increase the risk of both type I and II errors and underestimate parameter uncertainty.
However, autocorrelation also offers statistical opportunities.Quantifying the autocorrelation pattern of an ecological process offers insight into the studied system (Hamel et al., 2012).Also, statistical models that acknowledge dependency between adjacent time steps facilitate borrowing of strength over time.That is, by acknowledging that the behavior of a process at time t +1 exhibits similarities with the behavior at time t, time-specific parameters at time t +1 can indirectly be informed by data for time t and vice versa.This is particularly useful for analyses of weak data.
In this study, we developed a Bayesian framework for analyses of longitudinal hunting report data.The model accounts for autocorrelation at the national, county, and HMP levels, and we compared the framework to equivalent models applied to independent years.The aim was to define a statistical model for the available data, solve computational issues, estimate unreported harvest at all levels of interest, and evaluate the potential benefit of including dependence between years in the estimation.

Study species
This study focuses on four game species representing different taxonomic groups, hunting practises, and spatiotemporal dynamics.
Red fox (Vulpes vulpes, henceforth "fox") is a small predator that is distributed and hunted throughout Sweden.It is hunted by specialized trappers for its pelt, but also to reduce its predation on other game species.
Wild boar (Sus scrofa) was present in Sweden historically but was eradicated.The present population stems from escapes from enclosures during the 1980s and population densities have increased rapidly.Wild boar is present in the southern and middle parts of Sweden and is harvested in large numbers.
Common eider (Somateria mollissima, henceforth "eider") is a marineliving diving duck, meaning that only hunters with access to the coastline have an opportunity to hunt this species.It was historically an important source of meat and down for people living in the Swedish archipelago, but interest from hunters and harvest have diminished over time.
Grey partridge (Perdix perdix, henceforth "partridge") is only found in agriculture areas.It is one of the few bird species that can be released without permission from authorities.The release of captive-reared birds may result in an extremely patchy harvest where some hunting teams with small hunting areas may harvest large quantities whereas similarsized neighboring teams have no harvest.

Data
Starting in 1938, the NGO Swedish Association for Hunting and Wildlife Management (SAHWM) has a public commission from the Swedish Government that includes estimation of annual harvest for all game species with voluntary reporting (presently 32 bird species and 14 mammal species).Annual harvest estimates are thus available since 1939 for many game species.Earlier harvest estimates were calculated per county, but from 1995 and onward, estimates are instead calculated per Hunting Management Precinct (HMP, in Swedish Jaktvårdskrets, a geographic division based on the SAHWM organization, the total number of HMPs varies slightly over time, but there are ≈ 320) (Fig. 1).The main reasons for this change were to increase the homogeneity of conditions within each area of estimation and to provide local management with data at finer spatial scales.
In this study, we included data from the hunting years (July 1 -June 30) from 1997/1998 to 2019/2020.We will henceforth refer to the hunting years by their start year, i.e. 1997-2019.Hunting teams report their annual total harvest to SAHWM on a voluntary basis.Reporting on paper forms occurs, but most commonly harvest is reported directly into the SAHWM-owned database Viltdata (www.viltdata.se).Hunting teams can choose to report their harvest continuously during the hunting year or on one occasion at the end of the hunting year.In the Swedish legislation system, the hunting rights belong to the landowner.As a result, the same families, groups, or organizations may hunt on the ground for generations, and no other hunting teams are allowed to hunt there.
All reports are checked by local personnel of SAHWM, who, if needed, contact the reporting person for clarifications.Examples of inconsistencies that are scrutinized are unusual harvest numbers, nullharvest (report of zero felled game) in areas where the species in question is common, or reported harvest of a species outside its normal distribution area.The reported information includes hunting area (here and henceforth meaning the size of the area for which each team holds the hunting rights), number of individuals harvested for each species, and the HMP in which the hunting ground is located.The geographic position of teams within the HMP is not known, and HMP is the finest spatial resolution that can be considered from the data.All individual reports are confidential, and data is only presented at the HMP and higher levels.For the hunting years 1997-2019, data consists of 131,420 reports.

Fig. 1.
Panel A shows Sweden with county borders indicated in blue and Hunting Management Precinct (HMP) borders of 2019 indicated in grey.Vårgårda HMP, highlighted in red, is used in panels B and C, to illustrate conceptually that the available data includes only the area and harvest of, in 2019, nine (though smaller ones are difficult to make out) reporting teams with areas from 0.01 to 16 km 2 and harvest of red fox ranging from 0 to 18.The shape and position (beyond the level of HMP) of reporting teams are not known.The number and area of non-reporting teams must be inferred from the reporting teams, and panels B and C illustrate two potential realizations with either overall similar sized areas (B) or with large differences between teams (C).Legend indicates harvest (K) for red fox, Vulpes vulpes.

Statistical models
We developed a statistical framework based on the trade-off between sufficient detail and computational pragmatism.In this section, we first introduce models for hunting area per team and hunting rate (i.e.number of animals harvested per year, which depends on the hunting area), focusing on a single HMP in a single year.For equational neatness, we exclude indexing for HMP and year until we, subsequently, extend these models to a hierarchical framework.

Hunting area per team
Previous analysis has demonstrated a nonlinear relationship between a team's area and hunting rate (Lindström and Bergqvist, 2020).Thus, accounting for how an HMPs total huntable area, here denoted H, is distributed across hunting teams is needed to estimate total harvest.Defining A i as the area reported by team i and as the proportion of huntable area covered by team i and where N is the total number of teams in the HMP, the distribution of area across teams was modelled as where a of dimension N is the shape vector, in which all elements are identical and denoted a ∈ R >0 .The reason for implementing the identical concentration parameter a is that we lack any team-level covariates that could explain hunting areas of specific teams, and we can only account for random variability.For a small a, there is a large difference in the relative area of teams within the HMP.For large a, hunting teams hunt on similar sized areas.However, the total number of teams is typically not known, necessitating estimation of the Q number of non-reporting teams. with Thus, Eqs. 3 and 4 sidesteps any concerns with individual xi and models the observable data conditional on two parameters, a and Q.The latter is inherently discrete, yet, for computational reasons (Stan's HMC sampler is not suitable for discrete parameters) and to promote a straightforward hierarchical model, we relax this constraint and define Q ∈ R >0 and N ∈ R >M .That is, we allow the non-reported area to be occupied by, e.g., 10.2 teams rather than 10.Yet, we always make predictions of unreported harvest for the entire non-reported area, using an approximation (see Section 2.5) that side-steps dependence on individual teams' areas.
Eq. 2 can be conditionalized into ĥ ∼ Beta(Qa, Ma) ( and Because HMPs vary in size and non-reported area, it would not be precarious to assume that e.g.HMPs within the same county are necessarily similar in terms of N (total number of teams) or Q (the number of non-reporting teams).Thus, to borrow strength between HMPs, we found it more sensible to construct a hierarchical model around the average area per team, defined as m = H/N = H/(Q + M).

Analysis of hunting rate
The statistical model for hunting rate was defined by three parameters: μ, modeling the average hunting rate per area unit of a team with average-sized area; ϕ, modeling the potentially nonlinear effect of area on hunting rate per team; and β, modeling how hunting rate variability (between teams within an HMP) scales with the average.To derive the model, we started with the assumption that the probability of the reported bag K i by team i can be modeled as where team specific rate parameter r i ∈ R >0 was modelled as Here, the shape parameter α i depends on We modeled the non-linear relationship between a team's area and its hunting rate with the parameter ϕ as where the division by average proportion, x, and multiplication with average team area m facilitates the interpretation of μ ∈ R >0 as the average hunting rate per area unit (km 2 ) of a team with average area.We excluded the possibility of a negative relationship between a team's area and hunting rate by defining ϕ ∈ R >0 .
Defining the Negative Binomial (also denoted Gamma-Poisson mixture) by its mean and shape parameter, Eqs. ( 7)-( 9) can be combined into (10) by integrating out individual r i .
For each species, we excluded from the analysis of hunting rate counties outside of the range where the game is hunted.The reason was primarily pragmatic; including a large number HMPs with null-harvest led to weak identifiability of the model.Though this could be solved with, e.g., informative priors, we deemed it justifiable to exclude HMPs in counties with exclusively null-harvest for all years 1997-2019.However, parameters modeling hunting area per team were still informed by all data.

Spatiotemporal hierarchical models
We implemented a hierarchical model structure for the parameters defining the model at the HMP level, m, a, μ, ϕ, and β.All these parameters were defined as strictly positive, and for simplified notation, we use where ω θ,t , λ θ,l(k),t , and χ θ,k,t are year, t, specific nationwide, county, and HMP level effects, respectively, with l(k) indicating the county in which HMP k is located.We limited a,ϕ, and β to county-level parameterization, i.e. assuming that these higher-order parameters are identical across all HMPs within a county and year, and defined Thus, we assumed that variability in hunting area per team, nonlinear effect of hunting area on hunting rate, and within-HMP variability in hunting rate, as modelled by a,ϕ, and β, respectively, were equal for all HMPs within a county and year.The reason was primarily T. Lindström and G. Bergqvist pragmatic (computational), aiming at reducing the dimensionality of the model.We deemed it more important to prioritize HMP-level parameters for the most central parameters for harvest estimation, i.e., averages m and μ.We used posterior predictive p-values to investigate potential mismatch between observations and predictions, which would indicate that the model lacks important features.
This model structure of Eqs.11 and 12 facilitates borrowing of strength between counties and, in Eq. 11, between HMPs within counties.To incorporate temporal dependence at different spatial levels, we modeled first-degree autoregressive (AR) effects, such that the effects of one year can be written as conditional on the previous.For nationwide effects, this was modeled as where σ ω,θ models the similarity between adjacent time steps.The parent node for ω θ,1 is the prior, which is defined in Appendix B. For county and HMP effects, we obtained an identifiable model by defining and, when applicable, Here, ρ λ,θ , and ρ χ,θ are the county and HMP level correlation parameters and σ λ,θ , and σ χ,θ are the corresponding standard deviation of effects around zero.The interpretation of Eqs. ( 13)-( 15) is that the yearly effect of the focal parameter θ differ from the previous year with standard deviation σ ω,θ .For small values, the focal parameter changes little over time at the national scale.The county (λ θ,l(k),t ) and HMP (χ θ,k,t ) effects are distributed around the yearly effects with standard deviation σ λ,θ and σ χ,θ , respectively, and small values means that there is little difference among counties or among HMPs within counties.To account for the tendency of counties and/or HMPs to remain at a similar difference from the national average over time, ρ λ,θ , and ρ χ,θ models the autocorrelation in these effects.For instance, for ρ λ,μ close to one, there is a high consistency over time in terms of which counties hunt the focal game at high or low rates.For ρ λ,μ close to zero, hunting intensity would change randomly among counties between years.
In most cases, χ ⌣ θ,k,t− 1 = χ θ,k,t− 1 , but HMPs are occasionally changed between years, necessitating special considerations.There are three types of changes in the data: merging of two or more HMPs into one larger HMP (e.g. two becomes one), splitting of one HMP and merging the parts with other HMPs (e.g. three HMPs become two), and complete redrawing of HMPs within a county.In case of the two former, we defined where u is the number of obsolete HMPs merged into the new HMP k and p j,k is the proportion of HMP j incorporated into HMP k.
Complete redrawing of HMPs within a county only occurred twice in the data: the county of Södermanland between the years 2010 and 2011 and Gotland between 2007 and 2008.In these instances, we treated the HMP level effects as unrelated to previous HMPs by modeling χ θ,k,t according to the t = 1 option in Eq. 15.However, temporal autocorrelation at the county-level was still included.
To provide overall estimates of parameters we defined i.e. the weighted (by area) mean estimate.For m and a, the summation included all HMPs and years.For β and ϕ, the summation excluded HMPs located outside of the hunting range of the species.To provide comparable measures of nationwide hunting rate, we included HMPs outside of the counties where hunting occurs as μ k,t = 0 when defining μ.
We define prior distributions in Appendix B. For θ parameters, we obtained prior distributions for plotting by kernel-smoothed samples from the prior.
To evaluate the benefit of accounting for autocorrelation when estimating harvest, we analyzed the reporting data for each species and data configuration (full data or split into training and validation sets, see Section 2.6) in two ways: jointly with all 23 considered years with autoregression (AR) or as 23 independent year (IY) analyses, each focusing on a separate year.The latter obviates all ρ parameters and implements models strictly according to the t = 1 condition in Eqs. ( 13)-(15).

Computation
We used Stan (Carpenter et al., 2017;Stan Development Team, 2021) and the R-package rstan (Stan Development Team, 2020) to sample from posteriors.Stan's Hamiltonian Monte Carlo algorithm efficiently samples from continuous target distributions and often outperforms other MCMC software (Monnahan et al., 2017).However, efficient sampling may require reparameterization of the model.Specifically, we found that applying non-centered parameterization and avoiding sharp gradients were necessary to facilitate computation.Also, we avoided repeatedly evaluating the same calculations whenever possible.
Divergent transition, i.e., when the numerical integration breaks down due to the simulated trajectory of the HMC deviating from its true trajectory (see chapter 15 in Stan's refernce manual for details, (Stan Development Team, 2021)), posed a problem for some analyses, and we addressed this by gradually increasing the targeted average acceptance probability in Stan's sampler.For diagnostics, we used scale reduction factor (SRF) to investigate if chains had converged and effective sample size (ESS) to estimate how many independent samples the (typically autocorrelated) posterior samples represent (Gelman et al., 2004).Computational details and results are elaborated on in Appendix C.

Estimating unobserved harvest
We focus on prediction at the HMP level and again drop indices for HMP and year.Assuming gamma-distributed variables with equal β, the bagging rate of all non-reporting teams are where accent ∧ indicate parameters of non-reporting teams.Defining the average hunting rate of the non-reporting teams as ν = ∑ νi /Q, the likelihood of the total unreported bag K is Parameters Q and β were directly inferred by the statistical model, but νrequires further deduction, which was done in two steps.First, we derived an approximation for x, the unreported area proportions, based on the assumption that Q is large, in which case the distribution's properties could be approximated by a continuous distribution, f(x).
Secondly, we derived an expression g(ν) as a function of f(x), from which νis given as the first moment as which is finite for a and ϕ ∈ R >0 .This is derived in Appendix A.

Evaluating the approximation of K
To investigate if the approximation was "good enough" for the purpose, we performed a simulation study for selected HMPs.For each species, we selected three 2019 HMPs based on estimates of Q being low (median estimate of Q at the 2.5th percentile of all 2019 median estimates of Q), intermediate (median estimate at the 50th percentile), or high (median estimate at the 97.5th percentile).For each exemplifying HMP, we performed both exact and approximate simulations of K based on posterior samples.Exact simulations require discrete Q, and we rounded sampled Q to the nearest integer when evaluating the approximation.For the approximate simulations, K was directly simulated according to Eq. 19 for each posterior sample.
For exact simulations, the unreported harvest was simulated explicitly for each non-reporting team.This required the proportion of unreported area covered by each non-reporting team to be sampled according to Eq. A.4. Subsequently, each team was assigned a unique rate of harvest according to Eq. 8, and the harvest was simulated according to Eq. 7. Finally, K was calculated as the sum of all simulated teams' harvests.
To minimize Monte Carlo errors in the comparison, which could obscure small differences, we performed 100 simulations for every posterior sample.

Model comparison and validation
To investigate if AR modeling reduced uncertainty in estimated harvest, we calculated for each level of interest (nation, county, and HMP) the average ratio of the width of the 95% credible intervals estimated with AR over the corresponding width for IY estimates.We refer to this statistic as the 'range ratio', with values below one indicating less uncertainty for AR estimates compared to IY estimates.
To evaluate predictive performance, we split the data into training and validation sets by randomly assigning (with equal probabilities) each report to either set.In analyses of these data, only the training data were included for parameter estimation.To evaluate predictive performance, we subsequently performed predictions of harvest for the area covered by the validation data, imitating the issue of an unknown number of teams covering the unreported area.Using Eqs.19 and 20, we evaluated for each HMP and year with at least one report assigned to the validation data the probability of observing harvest K in the validation data.We define MPM as the mean predictive mass for each HMP, which integrates over the posterior by taking the mean value of posterior predictive probability over all post-warm up samples from the MCMC.
The same reports were assigned to either the training or the validation set, making it possible to compare predictive performance of AR and IY modeling.Akin to the log pointwise predictive density commonly used in Leave One Out Cross Validation (LOO-CV) (Vehtari et al., 2017), we considered HMP-wise difference in log MPM between AR and IY models and defined ALMPMD (average log mean predictive mass difference) as the average (over validation HMPs) difference between models.We also calculated the associated standard error to provide an estimate over how variable the difference in predictive performance is across HMPs.Further, to provide an estimate of total support of one model over the other, we calculated the product (over HMPs and years) of MPM ratio (AR over IY), abbreviated MPR.
We also utilized the validation data to evaluate posterior predictive p-values.This is elaborated on in the supplementary information, Appendix D.

Results
In this section, we first present the results of the parameter estimation stage, focusing on estimates from the AR model applied to the full data.Subsequently, we present results from posterior predictions based on parameter estimates, including for comparison also the IY analyses and the training and validation data application.Finally, we evaluate the approximation used for posterior prediction and to calculate the probability of the validation data.Analyses of posterior predictive pvalues indicated a good fit between model predictions and data.See Appendix D for details.
Computational diagnostics revealed that all chains converged and gave acceptable effective sample size (ESS).Computational performance was consistently better for AR than for IY analyses, both in terms of scale reduction factor and ESS.For IY analyses, there were also MCMCs where divergent transitions could not be circumvented by increasing the targeted average acceptance probability.These years were excluded from all further analyses.Computational aspects are elaborated on in supplementary information, Appendix C.

Parameter estimates
Figs. 2 and 3 present posterior estimates of high-level parameters modeling hunting area and hunting rate, respectively.For the considered period 1997-2019, the average area per team, m, was estimated at 33.2 [33.1, 33.8] km 2 (presented here and henceforth by median and 95% central credible interval in brackets).The average concentration parameter, a, modeling variability in area per team, was estimated at 1. 42 [1.41, 1.43].As a reference, a = 1 indicates an expected coefficient of variation of 1, and a > 1 indicates that, at the HMP level, the standard deviation of hunting areas per team is smaller than the mean hunting area.
The standard deviation σ ω,θ (θ again indicating either μ, ϕ, or β) models differences between subsequent years in terms of national averages, whereas σ λ,θ and σ χ,θ parameters models variability in space as the difference among counties and among HMPs within counties, respectively.Across all parameters and species, we found support for a smaller σ ω,θ than the corresponding σ λ,θ and (when applicable) σ χ,θ with a probability > 0.99.This indicates that between-year changes in system are smaller than the differences that are observed across space.
Estimates of yearly variability (Fig. 3 D-F) were largely overlapping across species, but σ ω,μ estimates for wild boar stood out as higher than for other species (Fig. 3 D), indicating comparably large differences between years in terms of average hunting rate.Wild boar also obtained the highest estimates of σ λ,μ .These posterior estimates however largely overlapped with other species, except for fox.This indicates that, compared to other species, fox is hunted at more similar rates across counties.Fox also exhibited the lowest estimate of σ χ,μ (Fig. 3 G), indicating greater similarity also at the HMP scale.
To investigate if variability was smaller at the HMP scale compared  T. Lindström and G. Bergqvist to the county scale, we calculated (based on posterior samples) the probability P ( σ χ,μ < σ λ,μ ) , which estimated at 1.00, 0.914, 0.999 for wild boar, partridge, and fox, respectively, providing support (if however weak for partridge) for larger variability in hunting rate across counties than at the scale of HMPs within counties.The corresponding statistic for eider was 0.074, thus showing a weak trend in the opposite direction.Similarly, the probability P ( ρ χ,μ < ρ λ,μ ) was estimated at 1.00, 0.990, and 0.995 for wild boar, partridge, and fox, respectively, providing support for higher temporal consistency at the scale of counties compared to HMPs.For eider (P ( ρ χ,μ < ρ λ,μ ) = 0.421), no such support was found.
However, both ρ λ,μ and ρ χ,μ indicated high spatiotemporal autocorrelation for all species (Fig. 3 H and N), with the lowest estimate found for ρ χ,μ for fox at 0.976 [0.971, 0.980].Estimates of ρ χ,ϕ and ρ χ,β (Fig. 3 J and L) indicated high temporal consistency in terms of spatial differences in higher-order parameters ϕ and β.Spatial differences were also largely consistent over time for parameters modeling the distribution of hunting area per team, as is captured by estimates of ρ λ,m ρ χ,m , and ρ λ,a in Fig. 2 E.

Harvest estimates
Fig. 4 shows the yearly harvest as estimated by either AR or IY models.The AR model produced tighter predictive envelopes, particularly for eider, partridge and early years of wild boar harvest.The range ratio analyses (Table 1) showed that 95% predictive credibility intervals were substantially reduced at the national, county and, HMP levels.Notably, AR modeling reduced intervals to less than 20% of the corresponding IY estimate for eider at both national and county levels and to less than 30% for partridge at the county level.The training and validation data analyses showed that AR modelling also increased predictive accuracy (Table 1, ALMPMD and MPR).

Evaluating the approximation used for prediction
The parameter estimates and unreported area of the exemplifying HMPs used to evaluate the approximation of unreported harvest, K, are listed in Table A.2.The corresponding estimates of K are shown in Fig. 5.The bars indicating the approximate methods (red) are largely concealed by the bars indicating exact sampling (dashed black) due to almost indistinguishable estimates.

Discussion
In this study, we developed a framework for analysis of harvest data for systems where harvest is reported by area.The primary aim was to extrapolate from available information to predict annual unreported harvest at the national, county, and HMP scale.We found strong support for positive temporal autocorrelation across all considered parameters and species (Fig. 3).This is expected-areas with, for instance, high hunting rates one year are likely to remain high the next year, due to similar game availability or cultural differences in game preferences.Beyond the insight that the system exhibits spatiotemporal consistency, the presence of autocorrelation also permits borrowing of strengths between years.With limited information, it is beneficial to inform what happens the subsequent year by what happened the previous and vice versa.Incorporating autocorrelation (AR) into the model resulted in less uncertainty at all considered spatial scales as well as higher predictive performance (Table 1, Fig. 4).Thus, our study concludes that AR modeling substantially improves estimation of harvest rates.
From a wildlife management standpoint, the proposed method has obvious benefits.Estimates with narrow uncertainty envelopes will increase the possibilities to detect changes in harvest rate, potentially in response to, e.g., altered climate, habitat, population structure or hunting practises.Identifying shifts in harvest rate may warrant further investigations and/or changes in hunting regulations for that species.Accounting for autocorrelation particularly reduced uncertainty in harvest estimates for eider and partridge, species that are hunted at Fig. 4. Annual harvest, indicated by median and 95% credibility interval estimated with autoregressive (AR) and independent years (IY) models.

Table 1
Comparison of autoregressive (AR) and independent years (IY) models.Average log mean predictive mass difference (ALMPMD, parenthesis indicate standard error) and mean predictive ratio (MPR) evaluated over validation data after fitting to training data.Values above 0 (ALMPMD) and 1 (MPR) indicate higher predictive accuracy for AR compared to IY.The range ratio is based on analysis of all available data and is defined as the average ratio between the width of 95% credible intervals at different levels of interest (nation, county, and HMP) comparably low rates and with high variability.Substantial differences between AR and IY models are also seen for wild boar in the first considered years (Fig. 4 C), i.e., when wild boar were rare (Bergqvist and Elmhagen, 2016).This indicated that, for the purpose of estimating harvest, the benefit of AR modeling is most beneficial for game that is hunted variably and in low numbers.
The understanding that borrowing strength (also labeled Bayesian shrinkage) through hierarchical Bayesian models reduces uncertainty is well established in ecology (Webb and King, 2009;Iknayan et al., 2014;Lindström et al., 2015;Huang et al., 2020).For harvest data, He and Sun (2000) and Lindström and Bergqvist (2020) have presented methods to borrow strength in space.In this study, we demonstrated the potential for improving harvest estimation by also borrowing strength in time, particularly when data are weak.We believe this insight could be useful for hunting harvest estimation also in other systems, and our study demonstrates the value of systematic long-term data collection.
The framework also addresses aspects beyond average harvest rate, which promotes quantitative insight into hunting patterns.For instance, we found the most pronounced sublinear relationship between a team's area and its harvest rate, as indicated by ϕ, for eider.This may reflect how and where eiders are hunted; there is a large discrepancy between those who do and those who do not hunt by the coastline.Though the fractal properties of coastlines (Kappraff, 1986;Husain et al., 2021) may have non-trivial effects, a reasonable null-expectation is that access to a one-dimensional coastline is approximately proportional to the square root of the area, thus giving an expected value of ϕ = 0.5.Yet, most of posterior density of ϕ is located < 0.5 for eider and < 1 for all species (Fig. 3 B), necessitating other explanations.It is plausible that teams with small areas simply keep their areas small because it is located where game availability is sufficiently large.Alternatively, to compensate for lower game availability, hunters with smaller areas may increase their effort, which has been shown to be an important predictor of harvest rate where this information is available (MMcDonald and Har-riscDonald and Harris, 1999;Tomeček et al., 2015;Vajas et al., 2020).
The estimated harvest rate for eider varied substantially among HMPs (Fig. 3 M) within counties and among teams within HMPs (Fig. 3 C).Unlike other species, the estimated autocorrelation was not lower at the HMP level than at the county level (Fig. 3 H, N).These findings also agree with the expectation that eider hunting is driven largely by access to specific habitats that are stable over time.
The variability of hunting rate among teams was also high for partridge (Fig. 3 C).Partridge can be released for immediate harvest and this is done to such an extent that the hunting bag may increase at the same time as the wild population is decreasing (Ottvall et al., 2009).Similar practices occur also in other countries (Ewald et al., 2010;Delibes-Mateos et al., 2013).Thus, specialized teams release a large number of partridge that are subsequently harvested.Other teams do not hunt this game at all or only rarely.This leads to a large discrepancy at the local scale.The high posterior estimate of σ λ,β (Fig. 3 K) suggests that the release and harvest practice vary among counties and, as inferred by large ρ λ,β (Fig. 3 L), that the practice is consistent over time.
Conversely, fox was the most homogeneously harvested game at all considered scales (Fig. 3 C,G,M).Fox is a generalist that can adapt to a wide range of habitats (Larivière and Pasitschniak-Arts, 1996).It is an available game to most hunters.Across the period, fox was also hunted at the highest rate of the considered species(Fig.3 A), but wild boar harvest has surpassed that of fox in recent years (Fig. 4 C, D).
These examples demonstrate that higher-order parameters-beyond average harvest rates-gives additional insight for managers into harvest dynamics and hunting practices.We argue that this can be valuable when assessing the status of species.For instance, game could face greater risk of stochastic over-hunting if the harvest is driven by a few teams that harvest in high numbers (low β).While we here have focused on the average parameters, the framework explicitly estimates unique non-linear effects of hunting area on harvest rate (ϕ) and random within-HMP variability (β) for each year and county.Studying spatiotemporal patterns of these higher-order parameters can pinpoint shifts in hunting practices.These aspects of the model also have important implication when estimating average harvest rates, particularly when reporting is low; the harvest observed per area unit for a handful of reporting hunting teams may not be representative for the entire population of hunters within an HMP.The framework also provides insight into the distribution of habitats across teams.Some between-year variability was observed, both in terms of average area per team and concentration parameters, but this variability was small in comparison to the substantial differences among counties as well as among HMPs within counties-differences that remain constant over time (Fig. 2 E).This indicates that, in terms of the size of hunting area per team, the Swedish system is highly variable across the nation but mostly consistent over time.
Autoregressive modeling presented fewer computational challenges than IY (Table C.2) due to more well-behaved distributions.Certainly, it would be possible to address the challenges found in IY modeling by other means than increasing the targeted average acceptance probability in Stan's sampler.This could involve either parameter transformations or more restrictive prior distributions.However, special solutions for individual species or years would render the methodology cumbersome to implement for the 46 game species for which SAHWM are commissioned to estimate the harvest.Though running the analysis for multiple years takes longer than for single years, it reduces man-hours.
Through the approximation of unreported harvest, we sidestepped sampling of individual non-reporting teams.The closed-form solution also permitted us to evaluate the likelihood of total harvest of the validation area in the training and validation data analysis.Even though the approximation was derived under the assumption of many nonreporting teams, the simulation analysis revealed that the approximation produced indistinguishable predictions from those of the exact sampling even at low number of non-reporting teams (Fig. 5).We believe the methodology-inferring a likelihood for count data of an unknown number of unobserved sections based on observed sections given non-linear effects of section size-could be useful outside the scope of analyzing harvest reports.
Reporting and sampling biases are important to consider when analyzing hunting data (Aubry et al., 2020), and the Swedish system is not immune to these challenges.Non-response may be higher for null harvest (Aubry and Guillemain, 2019), which could overestimate the total harvest.There are two aspects of the Swedish system that may to some extent mitigate such bias.First, because harvest is reported at the team level, it reduces the impact of null-harvest of individual hunters.Second, reports are entered jointly for all species, meaning that species for which the bag is zero are reported whenever a team enters data into the system.Nevertheless, the underlying assumption that the nonreporting population does not differ from the reporting is an important caveat that should be acknowledged.Future studies should aim to address this.Though it is beyond the scope of this study, we believe carefully designed surveys could gather additional information about hunting rates.A benefit of the Bayesian paradigm is its flexibility, and the presented framework could be expanded to incorporate such additional data.We believe the closed-form approximation is useful for such endeavors since it can be used to probabilistically define the hunting rate over Sweden's entire huntable land.Thus, it would be possible to incorporate in the framework a non-reporting factor that scales the rate of non-reported areas to best fit with additional, (hopefully) unbiased information.
Also, while reporting bias poses a problem for estimating the absolute harvest numbers, the framework is arguably less sensitive in terms of revealing trends, at least if reporting-bias does not change dramatically between years.The reduced uncertainty and increased precision obtained through AR modeling will increase the reliability of harvest estimates when they are used for wildlife monitoring.Previous harvest estimates have been used to analyze, e.g., climate effects on mammal populations (Elmhagen et al., 2015), population trends in goose species (Liljebäck et al., 2021) and effects of deer community composition on forest damage (Pfeffer et al., 2021).
We conclude that the presented framework provides important insight into hunting rates and practices.While the analysis was tailored for the focal data, we believe the benefit of accounting for dependence between years is applicable across systems.The study demonstrates the value of long-term harvest data, and by borrowing strength in time, information is transferred between years.We hope the framework will inspire new statistical approaches for valuable yet often imperfect citizen science data.

Declaration of Competing Interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Tom Lindstrom reports financial support was provided by Swedish Association for Hunting and Wildlife Management.

Fig. 2 .
Fig. 2. Parameter estimates and priors for parameters modeling the distribution of area among hunting teams.Parameter m (panel A) quantifies average hunting area per team and concentration parameter a (panel B) quantifies variability in hunting area per team within hunting management precincts (HMPs).For θ indicating either m or a, standard deviation σ ω,θ(panel C and D)  quantifies between-year variability in the focal parameter at the national level, σ λ,θ and (when applicable) σ χ,θ quantifies variability between counties and between HMPs within counties, respectively, and ρ λ,θ and ρ χ,θ (panel E) are the corresponding autocorrelation parameters, quantifying continuity in county and HMP differences over time.

Fig. 3 .
Fig. 3. Parameter estimates and priors for parameters modeling hunting rates.Parameter μ (panel A) quantifies average hunting rates per km 2 of hunting teams with average-sized hunting area, ϕ (panel B) quantifies how hunting rate per team depend on team area (with ϕ=1 indicating a linear response), and β (panel C)quantifying variability in hunting rates among teams within hunting management precincts (HMPs) (lower β indicating more variability).For θ indicating either μ,ϕ, or β, standard deviation σ ω,θ (panel D-E) quantifies between-year variability in the focal parameter at the national level, σ λ,θ (panel G, I, K) and (when applicable) σ χ,θ (panel M) quantifies variability between counties and between HMPs within counties, respectively, and ρ λ,θ and ρ χ,θ (panel H, J, L, N) are the corresponding autocorrelation parameters, quantifying continuity in county and HMP differences over time.

Fig. 5 .
Fig. 5.Posterior predicted distributions of unreported harvest, K, for three HMPs with low (top-row panels), medium (second-row panels), and high (third-row panels) estimated number of non-reporting teams, using either exact sampling or approximate method of simulation.
. Values below 1 indicate narrower predictive intervals by AR compared to IY models.