A Bayesian Hierarchical Model for Estimating American Woodcock Harvest

Estimates of total harvest help inform harvest management decisions, but such data are also useful for estimating population size and composition in demographic models. Historical estimates for U.S. harvest of American woodcock ( Scolopax minor ; hereafter woodcock) are available from 2 separate surveys: the 1964−2001 duck stamp survey (DSS) that sampled woodcock hunters who also hunted waterfowl, and the 1999−2016 Harvest Information Program (HIP) that sampled all licensed woodcock hunters. During overlap years (1999−2001), HIP estimates of total woodcock harvest were approximately twice as large as DSS estimates, but with only 3 years of overlap there was little potential to develop robust correction factors for historical DSS data. I developed a model of historical woodcock harvest that posited 3 groups of woodcock hunters, including those who always, sometimes, or never hunted waterfowl. During the HIP survey all 3 groups were included in harvest surveys; during the DSS years, however, only woodcock hunters who always hunted waterfowl were reliably sampled during all years, but I used annual duck stamp sales as a covariate to help predict harvest by woodcock hunters who hunted waterfowl irregularly. Using a reverse-time (2016 to 1964) model that assumed these 3 proportions of harvest remained constant through time, I was able to estimate total harvest in all years by estimating the latent component of harvest by non-waterfowl hunters. Averaged over all harvest jurisdictions, this model estimated that hunters who always, sometimes, or never hunted waterfowl contributed 43%, 32%, and 25% of the total woodcock harvest. Using these relationships, I estimated total harvest during all years (1964−2016) using data from both harvest surveys, although estimates based only on DSS data had greater uncertainty. In combination with band recovery data and harvest composition from the Parts Collection Survey, analysts could use estimates of historical harvest to estimate population size, composition, fecundity, and survival dating back to the initiation of harvest surveys in 1964.

Reliable estimates of total harvest are necessary for harvest management, but are also important for demographers because they can be combined with band recovery data to estimate fecundity, population size, and age-sex composition at time of banding (Zimmerman et al. 2010, Alisauskas et al. 2014, Hagen et al. 2018. For species such as American woodcock (Scolopax minor; hereafter woodcock), which are largely monitored through surveys of displaying males and for which polygynous mating systems may give a distorted view of effective population size of breeding females (Ziel et al. 2010), determining population composition is an important additional component of monitoring programs (Hagen et al. 2018). From 1964From -2001, the U.S. Fish and Wildlife Service (USFWS) estimated woodcock harvest in the United States using a survey of waterfowl hunters who had purchased a federal Migratory Bird Hunting and Conservation Stamp (Padding et al. 2010; hereafter "duck stamp survey" or DSS). But because woodcock hunters are not necessarily waterfowl hunters, the sampling frame for the federal DSS was incomplete with respect to woodcock harvest. More importantly, proportional coverage of active woodcock hunters by the DSS likely varied through time because more duck stamps were sold during years when waterfowl populations were high (Vrtiska et al. 2013), and the proportions of woodcock hunters who were inciden-tally captured by the sampling frame for waterfowl hunters likely varied on an annual basis. Moreover, these temporal relationships likely varied among states due to different opportunities and traditions for hunting waterfowl and woodcock over time.
In 1999, harvest surveys for waterfowl, woodcock, and other webless migratory game birds were modified to include samples of all licensed migratory bird hunters within each state (Padding et al. 2010). The joint federal-state Harvest Information Program (HIP) was specifically designed to identify an appropriate sampling frame for woodcock hunters and other groups of migratory gamebirds (Raftovich et al. 2015). Both surveys were conducted concurrently during 1999-2001, and estimates of woodcock harvest under the HIP survey averaged 1.6-fold higher in the Eastern Management Region and 2.1-fold higher in the Central Management Region (Padding et al. 2010); these increases were undoubtedly attributable to inclusion in HIP of woodcock hunters who rarely or never hunted waterfowl.
My objective in this paper was to develop a Bayesian hierarchical model to estimate total harvest of American woodcock from 1964-2016 by achieving a mechanistically appropriate synthesis of incomplete harvest estimates based on the duck stamp survey and with more recent harvest estimates based on the HIP. By modeling harvest of woodcock hunters who never hunt waterfowl as a latent (i.e., unobserved) parameter, and by controlling for annual

Methods
Data sources and conceptual model For each state, I obtained annual estimates of total woodcock harvest under the HIP (1999−2016) and federal duck stamp survey (1964−2001), and annual duck stamp sales (1964−2011). There were strong positive correlations between duck stamp sales and DSS harvest estimates in many states (Table 1), even though a federal duck stamp is not required for woodcock hunting. This could occur because numbers of waterfowl and woodcock hunters have both declined, and declines in annual duck stamp sales (Vrtiska et al. 2013) are also indicative of declining numbers of woodcock hunters (Luukkonen and Frawley 2010). Alternatively, many woodcock hunters might hunt waterfowl irregularly, with much of the annual variation in estimated harvest from the DSS arising from incomplete and variable coverage of the sampling frame of woodcock hunters due to their variable participation in waterfowl hunting. From 1999−2016, the sampling frame for the HIP survey of woodcock harvest included all individuals purchasing a state hunting license who answered "Yes" to the question: "Will you hunt migratory birds this season?" The HIP survey thus potentially included nearly all licensed woodcock hunters, excluding only small numbers of juniors, seniors, and landowners who are exempt from license purchase (Padding et al. 2010). However, concerns about non-compliance with HIP registration by hunters, and the failure of some license vendors to ask all screening questions needed for stratification, raise concerns that HIP might produce a slightly biased or less efficient sampling frame (Ver Steeg and Elden 2002).
Reconciling DSS and HIP surveys requires an understanding of how the sampling frames differed between the 2 surveys: the DSS potentially included all woodcock hunters who also hunted waterfowl each year, and partially sampled woodcock hunters who occasionally hunted waterfowl ( Fig. 1). By contrast, the HIP survey had a high sampling intensity for woodcock hunters who successfully harvested woodcock during the previous season, and a smaller sampling intensity for hunters who did not hunt woodcock or were unsuccessful during the previous season. Reconciling differences between these 2 surveys might be achieved by recognizing 3 components to the annual woodcock harvest, including harvest by woodcock hunters who: 1) always hunted waterfowl, and were therefore present in both the DSS and HIP sampling fames during all years; 2) irregularly hunted waterfowl, and were consistently sampled by HIP but more likely to be absent from the DSS sampling frame during years with lower duck stamp sales; or 3) never hunted waterfowl, and were therefore only sampled during the HIP survey (Fig. 1) An annual woodcock harvest should be closely related to the previous year's harvest, due to year-to-year similarities in the numbers of active woodcock hunters and the regional abundance of woodcock. I therefore modeled annual woodcock harvest using an autoregressive model (Link and Barker 2010) that presumed total harvest would be related to harvest in the previous year. For 1999−2016, I assumed that the HIP survey measured all 3 components of total woodcock harvest, including harvest by hunters who always, sometimes, or never hunted waterfowl. For 1964−2001, I assumed that total woodcock harvest was partially measured by the DSS, including the full complement of woodcock harvest by hunters who always hunted waterfowl, but a variable fraction of woodcock harvest by hunters who sometimes hunted waterfowl. I indexed this fraction of occasional waterfowl hunters by using annual deviations from the long-term trends  in duck stamp sales within each state. I presumed that these long-term trends represented relative changes in total hunter numbers (i.e., declines in numbers of both woodcock and waterfowl hunters), whereas annual residuals represented relative participation by active woodcock hunters in harvesting waterfowl (and hence, probability of inclusion in the DSS).
Formal hierarchical model Because the most complete sampling frames of woodcock hunters occurred under the HIP survey, I developed models in reverse time, treating the 2016 harvest as year t = 1 and 1964 as year t = 53. For each state (s) and year (t), I treated true harvest (H) as an unobserved latent variable: where logH s,t is natural logarithm of true woodcock harvest in state s during year t and ε s,t+1 is the annual rate of change in true harvest from the previously modeled year. I used state-specific HIP estimates during 2012-2016 to develop priors for year 1 (2016) log harvest (logH s,1 ) and I modeled ε s,t+1 ) as random normal variables with potential positive or negative state-specific trends in annual harvest (∆r s ) and state-specific annual variation(σ s 2 ) from this trend: where ∆r s~N ormal(μ r ,σ r ) using vague uniform priors for μ r (-1,1) and σ r (0.01,1) and σ s 2~l ogNormal(μ σ , σ σ ) using vague uniform priors for both μ σ and σ σ (0.01,10).
For observation error, I used state-specific estimates of annual variance in harvest from HIP surveys (Padding et al. 2010, Seamans andRau 2016). If annual harvest and variance of annual harvest were estimated at 0 for a particular state due to a small sample of hunters, none of whom reported harvesting woodcock (e.g., FL 1999; P. Padding, US Fish and Wildlife Service, pers. commun.), I replaced these estimates with the minimum observed estimates of means and variances for that state to avoid taking the logarithm of 0. For states that were combined for analysis (e.g., Vermont + New Hampshire), component variances were added together to obtain total variance. Variance estimates are not calculated for the DSS, but Geissler (1990) estimated variance for a variety of waterfowl species and sampling units (e.g., flyways, states) using bootstrap methods (Efron and Tibshirani 1986). I summarized data from Geissler (1990: Tables 1-2) and estimated the variance to mean relationship for 86 harvest estimates using Taylor's (1961)  I modeled duck stamp sales as a continuous time series using the smooth.spline function in program R with 4 degrees of freedom, which captured long-term trends in total duck stamp sales within each state without removing shorter-term variation potentially caused by year-to-year variation in participation in waterfowl hunting (Fig. 2a). I used residuals from state-specific splines as an annual index of potential short-term participation in waterfowl hunting by woodcock hunters.
For 1999 where π s1 represents the proportion of harvest by woodcock hunters who always hunted waterfowl, π s2 represents the proportion of harvest by woodcock hunters who sometimes hunted waterfowl, frac st represents the fraction of occasional waterfowl hunters who were hunting waterfowl that year (as indexed by residual duck stamp sales), and Ŝ E DSS,s,t is estimated SE of the DSS harvest estimate, derived from Taylor's power law. I modeled harvest and duck stamp data in JAGS 3.3.0 (Plummer 2012) using the jagsUI package (Kellner 2015) as implemented in R, using an adaptation phase of 1,000 iterations, followed by 3 Markov-chain Monte Carlo (MCMC) chains of 110,000 iterations each, with the first 10,000 iterations discarded as burn-in, and retaining every tenth remaining iteration, giving 30,000 observations for each posterior distribution. I verified model fit by examining trace plots and verifying that all R values were <1.03 (Gelman and Rubin 1992). Data sets and R code for running all models are provided online .

Results
During 1964−2001, annual harvest estimates of woodcock for most states were strongly correlated with annual duck stamp sales (mean r = 0.55). This was most prominent in the Eastern Management Region (Table 1) and was typically driven by strong declines in both woodcock harvest and duck stamp sales (Figs. 2-3). Correlations between DSS and residual duck stamp sales were weaker (mean r = 0.36), but positive (r > 0) in all states (Table 1).
During 1999−2001, when the duck stamp survey (DSS) and Harvest Information Program (HIP) overlapped, estimated harvest components averaged 57, 26 and 16% for hunters who always, sometimes, or never hunted waterfowl, respectively, but with tremendous variation among states (Table 1). Four northern states, including Michigan, Minnesota, and Maine, had large (≥40%) components of woodcock harvest by hunters who never hunted waterfowl (Table 1), and these states therefore exhibited large discrepancies between HIP and DSS estimates during the 1999−2001 overlap years (Fig. 3). Seven states, including New Jersey, had >70% of estimated woodcock harvest by hunters who always hunted waterfowl, with <10% harvested by hunters who never hunted waterfowl (Table 1). These states demonstrated little difference between HIP and DSS estimates during overlap years. For Louisiana, more than 40% of the woodcock harvest was estimated to have come from hunters who irregularly hunted waterfowl, but given that 95% of irregular waterfowl hunters were estimated to have hunted during the 1999-2001 overlap years, the DSS and HIP surveys were similar (Fig. 3). Michigan and Minnesota also had large proportions of woodcock hunters who hunted waterfowl irregularly, and accounting for annual duck stamp sales changed the shape of harvest trajectories considerably from the raw DSS data (Fig. 3).
Harvest estimates for Michigan based on an independent survey conducted by the Michigan Department of Natural  Minnesota, Louisiana, New York, Maine, andNew Jersey, 1964−2005 (y-axis is on loge scale: 9, 10, and 11 represent ~8, 22, and 60 thousand stamps sold, respectively). Long-term trends were characterized using smoothing splines with 4 degrees of freedom (red lines). Fitted lines are presumed to represent long-term trends in numbers of small game hunters, whereas residuals are presumed to represent relative participation of small game hunters in waterfowl hunting (i.e., purchasers of duck stamps). For Michigan (top left panel), additional lines demonstrate fitted splines for 3, 6, and 10 degrees of freedom (black, blue, and green lines, respectively).

Harvest Estimation in Woodcock · Arnold
Resources during 1964-2015 were strongly correlated with model-based estimates from this study (r = 0.78, P < 0.001). For Minnesota, the correlation between independent harvest estimates by the Minnesota DNR and model-based estimates from this study was also strongly positive (r = 0.73, P < 0.0001).
Most state-specific harvest estimates followed a common pattern of increasing total harvest from 1964 through the mid-or late-1970s, followed by strongly declining harvest until 2016 ( Fig. 3; Supplemental Materials). Aggregated harvest, summed across all states within each management unit, exhibited similar patterns (Fig. 4). The Central Management Unit exhibited tremendous variation in annual harvest estimates during the late 1970s and late 1980s, and this variation was largely driven by variation in harvest estimates from Louisiana (Fig. 3).

Discussion
By hypothesizing that woodcock harvest consisted of birds shot by hunters who: 1) always, 2) sometimes, or 3) never hunted waterfowl, I was able to model latent components of total harvest corresponding to these 3 groups. During years when total woodcock harvest was measured using the federal duck stamp survey (DSS), only the first component and a variable but unknown fraction of the second component were measured each year. However, by model-  ing harvest by occasional waterfowl hunters as a function of annual variation in duck stamp sales within each state, I was able to obtain estimates of total woodcock harvest by regular and occasional waterfowl hunters. By estimating these 2 harvest components during 1999−2001, when the duck stamp survey and Harvest Information Program (HIP) were conducted concurrently, harvest of woodcock by hunters who never hunted waterfowl also became estimable. Furthermore, by using models where annual harvest was presumed to be correlated to harvest in adjacent years, I was able to discriminate between annual process variation (i.e., real annual variation in woodcock harvest) and measurement error (i.e., estimation uncertainty due to incomplete sampling of woodcock hunters). This approach allowed me to obtain seamless estimates of total woodcock harvest across both survey frameworks, 1964−2016, including measures of uncertainty during early years (Figs. 3−4) that had formerly included only point estimates (Padding et al. 2010).
In addition to the myriad assumptions necessary to avoid bias in harvest estimates (Sheriff et al. 2002), this model required several additional assumptions. Because separate components of harvest were only measured during the 1999−2001 overlap years, I had to assume that the proportion of woodcock harvest from hunters who never purchased federal duck stamps remained relatively constant during 1964−2001. I also assumed that short-term annual variation in duck stamp sales was positively correlated with the surveyed fraction of woodcock hunters who sometimes hunted waterfowl, and positive correlations between DSS-based harvest estimates and residual duck stamp sales suggested this assumption was true (Table 1). I further assumed that observation error from the DSS could be approximated using variance-mean relationships derived from waterfowl harvest estimates (Geissler 1990). An alternative model that considered state-specific variance inflation factors that allowed DSS harvest estimates to have greater than predicted sampling variances resulted in similar harvest estimates, proportions, trends, and measures of annual process variation (T. Arnold, unpublished data), suggesting that my assumption about observation error had little impact on harvest estimates.
The primary goal of this model was to obtain seamless estimates of total woodcock harvest that could be used for Lincoln population estimates of total and cohort-specific population sizes (Alisauskas et al. 2014;S. P. Saunders, et al., 2019). I believe these estimates are more reliable than simple survey-specific estimates from either the DSS or the HIP, but they could be refined based on better knowledge of hunter behavior (e.g., contemporary or historical estimates of proportions of woodcock hunters who always, sometimes, or never hunt waterfowl; and the relative harvest rates of each group). Although my primary moti-vation was to obtain long-term total harvest estimates for use in developing integrated population models (Schaub and Abadi 2011), further developments of this approach might be useful for understanding changing dynamics of woodcock hunters. Moreover, some state-specific harvest estimates from the HIP survey appear to be highly variable (e.g., Maine, Fig. 3c), and state-space approaches helped achieve more consistent annual harvest estimates even without the added goal of trying to synthesize estimates from the DSS and HIP surveys. Bayesian hierarchical models are highly flexible, and the preliminary models developed here could be readily modified to include other measured variables known or suspected to affect annual woodcock harvest (Luukkonen and Frawley 2010).