Hierarchical Bayesian Model for Estimating Migratory Bird Harvest in Canada

1 Adam C. Smith 2 Canadian Wildlife Service, Environment and Climate Change Canada 3 National Wildlife Research Centre 4 1125 Colonel By Dr. 5 Ottawa, ON, K1A 0H3 6 adam.smith2@canada.ca 7 8 Smith et al. • Bayesian Harvest Model 9 Hierarchical Bayesian Model for Estimating Migratory Bird Harvest in Canada 10 Adam C. Smith, 1 Canadian Wildlife Service, Environment and Climate Change Canada, 11 National Wildlife Research Centre, 1125 Colonel By Drive, Ottawa, ON K1A 0H3, Canada 12 Thomas Villeneuve, Canadian Wildlife Service, Environment and Climate Change Canada, 13 National Wildlife Research Centre, 1125 Colonel By Drive, Ottawa, ON K1A 0H3, Canada 14 Michel Gendron, Canadian Wildlife Service, Environment and Climate Change Canada, 15 National Wildlife Research Centre, 1125 Colonel By Drive, Ottawa, ON K1A 0H3, Canada 16 ABSTRACT The Canadian Wildlife Service (CWS) requires reliable estimates of the harvest of 17 migratory game birds, including waterfowl and murres, to effectively manage populations of 18 these hunted species. The National Harvest Survey is an annual survey of hunters who purchase 19 Canada’s mandatory migratory game bird hunting permit. We use these survey data to estimate 20 the number of birds harvested for each species, as well as hunting activity metrics such as the 21 number of active hunters and days spent hunting. The analytical methods used to generate these 22 estimates have not changed since the survey was first designed in the early 1970s. Here we 23 describe a new hierarchical Bayesian model, which replaces the series of ratio estimators that 24 comprised the old model. We are now using this new model to generate estimates for migratory 25 bird harvests as of the 2019-2020 hunting season, and to generate updated estimates for all 26

migratory game birds, including waterfowl and murres, to effectively manage populations of 1 8 these hunted species. The National Harvest Survey is an annual survey of hunters who purchase 1 9 Canada's mandatory migratory game bird hunting permit. We use these survey data to estimate 2 0 the number of birds harvested for each species, as well as hunting activity metrics such as the 2 1 number of active hunters and days spent hunting. The analytical methods used to generate these 2 2 estimates have not changed since the survey was first designed in the early 1970s. Here we 2 3 describe a new hierarchical Bayesian model, which replaces the series of ratio estimators that 2 4 comprised the old model. We are now using this new model to generate estimates for migratory 2 5 bird harvests as of the 2019-2020 hunting season, and to generate updated estimates for all We modeled the year-effects for both harvest and days using a random-walk, first-difference for the variance of the year-effects as weakly informative priors on the standard deviations, density at values < 1.0, but includes a relatively long tail that allows for much larger values, if prior is only very weakly informative (Gelman et al. 2006  We estimated the observation-level hunter-effect parameters for days and harvest as zero-mean, 2 1 2 t-distributed random-effects with caste-specific variances and degrees of freedom ). Using the t-distribution to model these over-2 1 4 dispersion effects allows the modeled hunter-level variation to fit heavier tails than a normal 2 1 5 distribution. We used a normal distribution in early versions of the model, but in most zones and 2 1 6 years the empirical distributions of these hunter-level effects showed much heavier tails than a 2 1 7 normal distribution. These heavy tails capture the influence of particularly active and successful gamma distributions with shape and scale set to 2 and 0.2 respectively.
We estimated the parameters of the zero-inflation for the harvest (ߩ ௭ ) using a logistic was a function of the value in the previous year, plus random variation. We estimated the binomial probability that a hunter in caste-c and year-y actively hunted ) using the known numbers of HQS respondents who purchased a permit in that year, 2 2 9 and the number of respondents who indicated they had > 0 days of hunting zero-inflation component.
In the same way, we estimated the binomial probability that a hunter in caste-c and year-y who 2 3 5 actively hunted was successful (i.e., harvested > 0 birds, ) using the known numbers of HQS respondents who indicated they had > 0 days of hunting, and the number of HQS 2 3 7 respondents who indicated they harvested > 0 birds (݊ sub-model identical to the one used for the proportion that were active.
Similarly, we used similar time-series, logistic regression models to correct for inter-zone hunting. We modeled the annual probability that hunters sampled in a given zone would hunt 2 4 2 mostly outside that zone (ߩ ௩ ) and the annual probability that a hunter hunting mostly inside 2 4 3 the zone was sampled outside that zone (ߩ ௩ ). We derived the data for these two sub-models given zone and year, we modeled the probability that a hunter sampled in that zone would hunt sampled in the zone who hunted outside the zone (݊ ).

4 9
We also modeled the annual values of ߩ ௩ and ߩ ௩ using a time-series logistic model, with a structure identical to the one used for the other binomial probabilities.

5 1
We combined the correction factors for inter-zone hunting, the proportion of active hunters, and ).
We re-scaled all mean values for hunters, castes, and species to population-level totals using the 2 5 5 estimated number of active hunters for a given caste and year ‫ܣ(‬ , ௬ ).

5 7
We calculated the estimated mean number of days hunting waterfowl (݀ , ௬ ) for caste-c and year-2 5 8 y as a derived statistic using the exponentiated sums of the relevant parameters, plus some added 2 5 9 variance components to account for the asymmetries in the retransformation.
is an approximation of the half-variance retransformation from the 2 6 1 mean of a log-normal distribution to the mean of the normal, accounting for the t-distributed . This term is an this re-transformation is an area of ongoing research, we have found that it generates estimates of 2 6 7 total harvest and activity that are on the same scale as estimates from the previous model. ). In this case, we also included the estimated probability of a 2 7 1 non-zero harvest to account for the zero-inflation. We divide the estimates of mean total harvest by year and caste into species-specific estimates of 2 7 5 harvest using data on the seasonal harvest patterns from the HQS calendars (daily estimates of

Harvest by Species and by Age and Sex
the number of waterfowl harvested for each respondent) and data on the species-composition 2 7 7 information collected from the SCS. On the SCS envelopes in which wing and tail fans are received, hunters also provide information on their current year permit numbers as well as the 2 7 9 date and location of harvest. These parts are then identified to species, and in most cases,  For calculating the species composition, we divided the harvest season into periods to account 2 8 5 for the declining response rates to the SCS as the season progresses (e.g., due to response- entire season (Cooch et al. 1978). We based the periods for a given zone and harvest group (e.g. sharing of information on the seasonal patterns in harvest across years, we kept the periods 2 9 2 consistent across all years. For ducks, we divided the season into 6 to 13 periods, depending on season. In most cases, the earlier periods in a zone each represented a single week and the later 2 9 5 periods may include more than 1 week. The final period often included many weeks, to capture 2 9 6 the low level of hunting that continues through the winter. The model includes three similarly structured sub-models that rely on multinomial distributions 2 9 8 to estimate three key proportional distributions: 1) the proportions of annual harvest that occurred in each period using the calendar data; 2) the proportion of the harvest in each period 3 0 0 that can be attributed to each species; and 3) the proportion of each species harvest that can be 3 0 1 attributed to the age and sex categories (e.g., adult-female, immature male, etc.). We then 3 0 2 combined these proportions with the total harvest estimates to estimate the number of birds category for a given species. For the sub-model for the proportional distribution of the harvest across W-periods, we used data ).
We gave Dirichlet-priors to the individual probabilities for each period-w and year-y   fixed-effect for all other periods using a normal prior with mean of 0 and a variance of 10 For the sub-model for the proportional distribution of all S-species in period-w and year-y, we  ).We estimated the parameters of the Dirichlet using a hierarchical, log-link, random- ) as well as a random effect for the mean abundance across periods for a given species (߶ ௦ , ௪ ).
In the first year, we estimated the value for the year-effect component in each period and species 3 2 4 as a random effect using a prior with mean = 0 and a variance specific to the period for the random species effect by period, we fixed it at 0 for the first species, and we estimated it 3 2 7 as a fixed effect for all other species, using a normal prior with mean of 0 and a variance of 10 Finally, we gave Dirichlet-priors to the individual probabilities for each demographic group-d 3 3 0 (i.e., each combination of the age and sex categories, such as adult-female, immature-male, etc.) harvest across the periods. Finally, for the derived estimates, we re-scaled the relevant mean values (e.g., mean harvest of ). We summed these derived estimates of harvest, days, number of successful hunters, etc. in the annual analysis.

4 4
Uncertainties for each of the estimates represent summaries of the full posterior distributions.

5 6
Other Species 3 5 7 As part of the HQS, we also collect information on the harvest of other species of migratory coot, rails, and band-tailed pigeon. We capture this harvest information using a simplified year-to-year for data-sparse species, including species of conservation concern (Fig. 5).

0 1
With the new model, we also generate formal estimates of age ratios in the harvest, which have 4 0 2 estimates of uncertainty and are similarly less sensitive to annual fluctuations (Fig. 6). We all lower than ratios from the old model (Fig. 7). This bias in the old estimates was due to the times as many parts in relation to the total harvest, than were received in western provinces such 4 1 0 as Saskatchewan (Fig. 7). concern" by the Committee analysis would estimate no harvest for less common species if no parts were submitted in a given 4 2 7 year, such as all national estimates for harlequin duck since 2012 (Fig. 5). Harvest estimates for from the traditional design-based approach, although generally more precise. This similarity suggests that there is no need to revisit past management decisions, and that these new estimates 4 3 1 can integrate seamlessly with existing decision processes for harvest management. We see this model as an initial step in an ongoing evolution of our estimation and understanding conditions during peak hunting periods. Additionally, the spatial stratification provided by the increase response rates for the survey and to generate more data to help fill this and other gaps.

6 7
With these additional data and by sharing information in a hierarchical framework and assuming Wildfowl 69: 206-220. Evolution 4: 132-143.  The good, the bad and the not great. Methods in Ecology and Evolution 11: 882-889. in Canada. Ottawa. modeling. Ecological Applications 19: 553-570.   for North American sea ducks. PLOS ONE 12: e0175411. Bayesian hierarchical spatial models: Implementing the Besag York Mollié model in stan.