This analysis used open (government website-posted) data of hunt unit-level white-tailed deer harvest and winter weather in northeastern North America from 2003-2024 to evaluate trends in deer harvest with respect to some commonly used deer Winter Severity Indices (WSIs). From these data I calibrated a novel formula for WSI (the LSDNT index). In comparison to other WSI methods, this new WSI demonstrated a modest increase in local deer harvest prediction and substantially greater prediction of deer harvest trends at the scale of the northeastern USA. After fitting harvest-based WSI models, I extrapolated deer harvest trends to a general predictive model of deer annual population impacts from winter weather.
This analysis provides a benefit to deer managers presenting a predictive context of winter severity at the scale of individual deer management units or finer scales. Managers can use these locally-predicted results for insight to questions such as:
How does the winter severity in this area fall within the overall deer-climate spectrum of the northeastern United states?
What are the likely effects and magnitude of a given winter-year’s weather on this area’s deer population?
Within this area, how many years are the effects of a given winter likely to influence hunter harvest or other deer population metrics?
These questions of winter weather effects are important considerations where deer managers confront issues such as disease, predation, habitat management, and harvest policy. Managers often look to guidance from research and outcomes in different jurisdictions (states or provinces), but directly accounting for winter severity in cross-jurisdiction comparisons is hindered because almost every jurisdiction in northeastern North America uses a unique system of quantifying deer winter severity. Therefore, a unified winter severity index would reduce barriers to cooperative deer research and management.
I assessed deer harvest and annual WSI values from 415 hunt units comprising > 7,000 hunting seasons across 7 USA states (MN, WI, MI, NY, NH, VT, ME). Deer harvest tended to decrease following more severe winters in 87% of hunt units (negative correlation), and this effect was statistically significant in 24% of hunt units. Hunt units in areas with more severe winters had greater harvest sensitivity to winter weather, and I used this relationship to classify 3 zones of winter severity (low, medium, and high). I further used the average WSI relationship with WSI effect size to generate annual predictive effect maps of WSI on deer in the northeastern USA.
The lagged effects of winter severity on deer harvest were statistically significant up to 4 years after a given winter. These lag effects of winter grew stronger and persisted for more years as the average winter severity of hunt units increased. Overall, annual winter severity and lag effects explained about 20-30% of annual variation in deer harvest across the study interval.
This approach to measuring deer winter severity may offer several advantages. First, the evidence in this analysis suggests that LSDNT index has greater predictive capacity than other common WSI across northeastern North America. Second, scaling predictions to local climate mitigates a problem of local baselines that comes with “absolute” WSI measures; i.e., a WSI score that is relatively mild in a northern Maine locality may represent a historically severe winter in coastal southern Maine. Finally, classification of local WSI trends across the northeastern USA aids comparison across jurisdictional boundaries when examining management and ecological interactions between climate and deer populations.
This analysis focused on adult male deer harvest, which was consistently offered to hunters as an “over-the-counter” option for all licensed hunters across the study area. As a result, this dataset was a variety of annual sampling without replacement for male deer >1 years-old. Certainly the data has many limitations that come with use of harvest as a population indicator. However, at very large spatio-temporal scales it is a plausible assumption that much, if not most, of hunter success where tag availability is constant is explained by changes in deer abundance. Therefore, WSI values indicating a positive or negative effect on male deer harvest are likely indicative of directional effects on the deer population as a whole.
The Linear Snow Depth, Nonparametric Temperature (LSDNT) Winter Severity Index is calculated as the average daily snow depth from October to May, plus weighted contributions from the number of days where the minimum daily temperature falls below certain thresholds:
\[ \text{WSI} = \overline{S} + 0.33 \cdot T_{-20} + 1 \cdot T_{-25} + 5 \cdot T_{-30} + 10 \cdot T_{-35} \]
Where:
Winter weather has many established effects on white-tailed deer populations, through which more severe winters (those with deeper snow, colder temperatures, and longer duration of winter conditions) suppress deer population growth. As such, deer management agencies widely use Winter Severity Indices (WSIs) to gauge potential consequences of annual weather conditions on localized deer populations.
Current methods for calculating deer WSIs are inconsistent among jurisdictional agencies in eastern North America, but most use threshold values of daily snow depth and temperature to calculate a cumulative annual winter score from late winter to early spring (approximately November-April). These indices have decades of use with proven relationships to deer population effects, but few attempts have been made to explore new methods for greater predictive value. Further, current WSIs, though similar, are subtly different among states with regards to calculation values, start and end dates, weather data sources, and interpretations (e.g., what annual value constitutes a mild or severe winter). These differences make general comparison of winter weather effects on deer populations at a continental scale difficult.
Among the obstacles to generating a more accurate winter severity index is the large spatio-temporal scale of winter severity effects on deer populations. Each year represents only a single winter value, and even long-term intensive studies of deer populations with telemetered samples and robust annual population estimates generally last less than a decade. Additionally, geographic differences in winter weather often are not easily contrasted until compared at large distances.
With the limitations of broad inference into deer and winter weather from planned observational studies, hunter harvest is a potential data source to compare winter weather effects across large spatio-temporal scales. Though hunter harvest data has many caveats as a representation of deer population change, it also has some strengths that no scientifically designed sampling efforts can rival. In the Eastern USA where deer hare harvested at high rates and resident hunter licenses are offered as over-the-counter opportunities, the annual hunting season represents a consistent sampling scheme without replacement spanning decades of data. This data set can be viewed in some respects as akin to the applications of harvest yields within fisheries science.
My objective for this analysis was to measure relationships between winter weather variables and hunter harvest in eastern North America as a test and calibration of deer Winter Severity Indices. Further, I sought to develop a broad framework of winter severity to facilitate comparison of deer populations among different jurisdictions through the lens of winter weather.
The target area of the analysis was the zone of white-tailed deer harvest across Eastern North America encompassing the, northern Great Lakes, Eastern Canada,and New England regions. This area represents a broad range of conditions in deer population density, habitat, winter weather, and management practices. Including a wide variety of deer population contexts was desired to assess WSI performance at a broad scale. However, deer management practices across these regions also share some characteristics useful for assessing deer population trends through harvest. First, all the jurisdictions in this region generally offer at least one antlered deer tag to all resident hunters each year. Though a very small portion of areas such are urban areas and military installations have exceptions to general resident hunting, it is a generally true condition that harvest effort for adult male deer is present across the targeted study zone with minimal influence from changes in tag availability.
My goal was to use adult male harvest rates as a base data to assess broad-scale WSI predictive value to deer populations. This approach has limitations. Despite consistently available antlered deer tags for hunters, many confounding variables exist that potentially disrupt the relationship between antlered deer harvest and overall deer population trends. First, hunter effort is variable among years are the number of hunters and days spent afield changes. Second, factors beyond hunter effort can affect harvest success such as changes to season length and dates, changes to allowed methods of harvest, and even weather conditions within season. Third, the abundance of adult male deer does not vary in lockstep with adult female and juvenile deer as age and sex ratios vary among years, so adult male abundance is an imperfect index of total deer population changes. These caveats are important limitations of interpreting adult male harvest as overall deer population trends.
Why use antlered deer harvest as a metric to assess population responses to winter weather? First, antlered deer harvests in most areas is heavily composed of 1.5 (yearling) to 3.5 year-old deer, with yearlings often comprising a third or more of antlered harvest. Fluctuations in yearling buck abundance are tied to the recruitment of the cohort’s year of birth and overwinter fawn survival of cohort’s first winter. Winter severity generally impacts overwinter fawn survival and fawn recruitment of the following summer, which may be expected to produce a signal in buck harvest that is detectable within a year (via yearling buck cohort) and may fade or vanish within 4 years (via 3.5 year-old buck cohort). This hypothetical relationship provides a tractable basis for time-series analysis of annual harvest in response to preceding winter weather. Although the aforementioned caveats in relating adult male harvest to overall deer population change are important to consider, is is also reasonable to expect that in most cases, the total harvest of adult males is positively correlated to deer abundance.
It is worthwhile to look for insights into deer-weather relationships within deer harvest because there is no other deer population monitoring data-set that compares to the spatial, temporal, and numerical scale of harvest data, especially when combining data across multiple jurisdictions such as states and provinces. Much research on deer and winter weather has leveraged intensive local sampling such as survival monitoring, movement analyses, and biological sampling (e.g. measuring nutritional reserves and reproduction. While a very large telemetry study may provide a high resolution sample for up to several local deer populations and last a decade or more, with direct live deer monitoring it is ultimately very difficult to sample more than perhaps 10–20 local population-level responses to winters. By contrast, hunt-unit level deer harvest can provide a low-resolution sample from hundreds of local deer populations across several decades, comprising thousands of local deer deer population-level responses to winters. To augment knowledge of deer winter ecology gathered in intensive local field studies, there is value in a large-scale perspective that captures extreme severe and mild winters conditions across many local climates, management regimes, and deer habitats.
Annual Deer harvest estimates were collected from summary reports published publicly on jurisdictional agency websites as available for the year range of 2004-2023 within Minnesota (MN), Wisconsin (WI), Michigan (MI), New York (NY), New Hampshire (NH), Vermont (VT), Maine (ME), Ontario (ON), Quebec (QC), and New Brusnwick (NB). Each jurisdiction reported deer harvest for management sub-units within the jurisdiction, which have various names (e.g. wildlife management units, hunt zones, wildlife management districts) but are hereafter collectively referred to as hunt units in this analysis (see Figure 2 for a map of hunt units).
All jurisdictions except for Maine provided annual estimates at the hunt unit scale. Maine provided annual deer harvest at the township level. I converted Maine harvest by assigning harvest to the hunt unit (Maine Wildlife Management Districts) at the geographic center of each township. Adult male-specific harvest rates were available and used for all jurisdictions except for Maine and New Brunswick, which only provided estimates of total deer harvest. Therefore, all deer harvest was used for Maine and New Brunswick, which provided a sub-optimal response co variate in those jurisdictions because agencies often change antlerless tag quotas in response to winter conditions, with could exaggerate the effect of winter severity on deer harvest. By comparison, buck-only harvest is typically a consistent opportunity for all hunters during all years. Nonetheless, I considered inclusion of “all deer” harvest for these units acceptable as the analysis broke down weather effects by jurisdiction, so specific trends to ME and NB could be identified and interpreted through the lens of “all deer” harvest.
Table 1. Sample sizes and date ranges for white-tailed deer hunt unit-level harvest data included in the analysis.| Jurisdiction | Harvest data | Data range | Harvested deer | Number hunt units | Hunt units analyzed | Season-years |
|---|---|---|---|---|---|---|
| NY | adult male | 2004-2023 | 2,149,176 | 89 | 89 | 1,780 |
| VT | adult male | 2006-2023 | 160,715 | 20 | 19 | 380 |
| NH | adult male | 2004-2023 | 140,561 | 20 | 20 | 340 |
| ME | all deer | 2005-2023 | 518,362 | 29 | 28 | 504 |
| MI | adult male | 2004-2023 | 3,616,306 | 85 | 78 | 1,437 |
| WI | adult male | 2007-2023 | 2,464,612 | 70 | 70 | 1,120 |
| MN | adult male | 2010-2023 | 1,368,673 | 113 | 111 | 1,443 |
| USA total (model fitting) | NA | NA | 10,418,405 | 426 | 415 | 7,004 |
| ON | adult male | 2008-2023 | 554,183 | 112 | 68 | 816 |
| QC | adult male | 2004-2024 | 705,208 | 18 | 12 | 228 |
| NB | all deer | 2011-2023 | 89,346 | 27 | 20 | 180 |
| Canada total (extrapolation only) | NA | NA | 1,348,737 | 157 | 100 | 1,224 |
Hunt units were removed from the data-set if they did not have at least consecutive 8 years of harvest data and a minimum of at least 15 deer harvested during any year. This ensured that the included hunt units had a meaningful time series of winter conditions and were less susceptible to extreme proportional changes in annual harvest due to very small harvest estimates.
After screening for short time series and very small harvest estimates before using harvest as a response variable to winter severity measures, it was necessary to account for the large differences in harvest magnitude among hunt units because using the absolute harvest values would allocate greater weight to units with greater harvest. and the potential effects of long-term harvest trends operating on larger timescales than winter weather (e.g. 10-20 year changes in deer harvest due to steady growth/decline in deer populations or hunter participation). To accommodate long-term harvest trends, each hunt unit was fit to a linear regression of harvest ~ year, from which the residual values of harvest were extracted. These residual harvest values were then standardized within each unit to proportional change from the previous year, such that the absolute magnitude of harvest was replaced with a relative change in harvest.
With controls in place for long-term linear harvest trends and proportional scaling to unit harvest magnitude, a pre-screening was made to identify and remove anomalous hunt units that would exert excessive leverage on winter weather estimates. To identify units with outlier relationships between harvest and winter severity that could confound assessment of general winter weather effects, I fit an initial General Additive Mixed Model (GAMM) of residual annual harvest change (the proportional change in harvest from the the preceding year, after accounting for linear trend of time over the study interval) in relation to average October-May snow depth. Average snow depth was selected as a representative measure of WSI for the sake of screening, as snow depth is the primary basis of most established deer WSI methods. I assessed the per-unit GAMM residuals to select a reasonable outlier cutoff, and I excluded from further analysis 4 hunt units which showed major divergence from the rest of the data-set. The number of hunt units included in the analysis are summarized in table 1. For the USA training data, this was 415 hunt units after pre-screening the data.
Figure 1. Hunt unit residual values from a GAMM fitting residual annual harvest change (the proportional change in harvest from the the preceding year, after accounting for linear trend of time over the study interval) in relation to unit-scaled average October-May snow depth. Hunt units with a residual > 0.4 were ecluded from further analyses due to exceptional leverage.
I included a variety of candidate WSI calculation methods for comparison of predictive value in deer harvest change. There were two prevailing WSI forumulas currently used by state agenciesin the eastern USA. Minnesota and New York use a furmula that accumulates 1 point for each day of snow depth > 15 inches and 1 point for each day with a minimum temperature less than 0 F (maximum daily value would be 2 when both conditions are met). Wisconsin, Michigan, Vermont, and New Hampshire use a similar formula that accumulates 1 point for each day of snow depth > 18 inches and 1 point for each day with a minimum temperature less than 0 F. I included both formulas as WSI candidates, named “S15_T0” and “S18_T0”, respective to their threshold values. These states employ variable start and end dates to winter accumulation, but all start in late fall (Oct-Nov) and end in mid-spring (Apr-May). In my analysis, I used Oct-May as the interval for both formulas.
In addition to the two commonly used WSI methods, I included several measures of winter severity that have been used in deer research. these included average snow depth, and days of snow cover, and average minimum temperature. I used Oct-May for the annual interval on these indices as well. For days of snow cover, I set a minimum threshold of >2 mm to define snow coverage, as snow depths less than this would be trace “dusting” conditions with no likely implications to deer.
I made an initial comparison of a priori WSI fit using a linear mixed-effects model (LMM). Proportional annual harvest change was the response variable, consisting of 7,004 hunting season harvest estimates from 415 hunt units within the seven USA jurisdictions. Each model included a random intercept for hunt unit and jurisdiction, and a fixed effect for the the respective WSI formula calculated annually for each hunt unit. From this initial analysis, average daily snow cover had greater statistical support than threshold-based WSI methods.
After the initial fitting and AIC evaluation of WSI methods, there was not a clear best WSI formula. In both single year and lagged year model approaches, certain indices best predicted harvest within states but there was relatively even support for average snow depth, inch18_F0, and inch15_F0. Among these three competing indices, it was noteworthy that average snow depth included only one weather input (snow depth), while the other two metrics included both snow depth and temperature data. Considering that average snow depth was a competitive WSI metric to snow and temperature threshold indices, I hypothesized that adding a weighted temperature component to average snow depth may produce a more predictive WSI.
I incorporated a post-hoc method to fit additional temperature weights to the average snow depth index using an information-theoretic approach and and underlying non-parametric temperature function. The base of this approach was the average snow depth model from the initial WSI comparison analysis. From this base model, I fit structured repeated iterations of permuted minimum temperature thresholds and weights, comparing each iteration using AIC until convergence was achieved at minimum AIV value. In essence, this procedure expanded upon the temperature threshold premise already used by some WSI formulas, but alowed for multiple thresholds of varying weight and tested likelihood ratio for combinations of temperature thresholds and their respective WSI weights against the deer harvest data. I allowed up to 5 temperature thresholds, because adding more thresholds was too demanding for the computational power available. Each iteration applied semi-random changes to temperature threshold values and WSI weights, then selected the model with lowest AIC as the base for the next iteration. Initial iterations used larger step sizes (e.g., +/- 20%) and to escape local likelihood maxima. As convergence approached, random step size decreased (e.g., +/- 1%). A final convergence was reached after ~50,000 iterations.
This process is similar to a visit to the eye doctor where patients are screened though a series of corrective lens combinations and declare each change as “better” or “worse”. Essentially, a very large number of binary comparisons suggested a most AIC-supported combination of temperature thresholds and weights to add to the average snow depth.
The top supported WSI formula from this process included 4 minimum daily temperature thresholds at -20°, -25°, -30°, and -35°C. The overall pattern suggest that temperature weight was unimportant until temperatures fell below -20°C, and then the weight of temperature increased exponentially to -35°C. Truncating temperature weights at <-35°C into one value likely does not suggest that temperature effects on deer do not continue to increase at even lower temperatures. Rather, it is likely a reflection that days with lower temperatures than -35°C were very rare and therefore additional low thresholds contributed little additional data to WSI prediction.
The result of this post-hoc WSI calibration was to treat daily snow depth as a linear negative effect on deer, while treating temperature as a non-parametric variable that roughly imitates exponentially negative effects on deer when temperatures fall below -20°C. Due to the combined functions use to weight the index, I refer to this as the Linear Snow Depth, Nonparametric Temperature (LSDNT) index.
An ecological interpretation of the LSDNT index is that deeper snow has greater negative effects on deer, and these negative effects are better described by the actual snow depth than by binary depth thresholds (e.g., above or below 18 inches). Low temperatures also have negative effects on deer, but direct impacts of temperature may be minimal at >-20°C, then become very strong at extreme lows (e.g., <-35°C). Inherent to this relationship is the understanding that temperatures below freezing contribute to sustaining greater snow depth, but that contribution is accounted for within the snow depth variable.
Figure 2. Effect coefficients of LSDNT, scaled within hunt unit, on annual rate of change in white-tailed deer harvest in northeastern USA, 2003–2023. Vertical dashed lines and labels Z0–Z3 represent the ad-hoc delineations of WSI zones used in this analysis.
Figure 3. (Interactive): Unit map of Effect coefficients of LSDNT, scaled within hunt unit, on annual rate of change in white-tailed deer harvest in northeastern USA, 2003–2023. Select a harvest unit for correlation details.
Deer harvest was log-transformed before fitting to mitigate units with greater harvest carrying disproportionate estimate weight.
Figure 4. Effect estimates of scaled LSDNT WSI lagged 0–4 years on log-transformed white-tailed deer harvest in northeastern USA, 2003–2023. Error bars represent 95% confidence intervals. The current year represents the Oct-May period immediately preceding each hunting season (e.g., Current year WSI for fall 2023 hunt would be Oct 2022 – May 2023).
Figure 5. Histogram of r values (Pearson correlation) between hunt unit annual observed deer harvest and predicted deer harvest from the 4-year lag model. Fill color represents jurisdiction (state) of each hunt unit. Harvest was log-transformed for model fitting and prediction, but back-transformed to absolute harvest numbers before computing r values.
Figure 6. Histogram of r values (Pearson correlation) between hunt unit annual observed deer harvest and predicted deer harvest from the 4-year lag model. Fill color represents WSI zone of the each hunt unit. Harvest was log-transformed for model fitting and prediction, but back-transformed to absolute harvest numbers before computing r values.
| WSI Index | Calculation Description |
|---|---|
| null_model | No WSI parameter |
| snow_cover | Number of days Oct-May with snow depth > 2mm |
| snow_depth | Average Oct-May daily snow depth (mm) |
| temperature | Average Oct-May daily minimum temperature (C) |
| inch15_F0 | Number Oct-May days with snow depth > 15 inches + Number of Oct-May days with minimum temperature < 0 Fahrenheit |
| inch18_F0 | Number Oct-May days with snow depth > 18 inches + Number of Oct-May days with minimum temperature < 0 Fahrenheit |
| LSDNT | Average Oct-May daily snow depth (mm) + daily weighted nonparametric temperature values (see overview) |
| Index | MI | MN | NH | NY | VT | WI | ME | All |
|---|---|---|---|---|---|---|---|---|
| LSDNT | 9.5 | 11.0 | 0.0 | 6.1 | 7.0 | 0.0 | 0.0 | 0.0 |
| inch18_F0 | 0.0 | 79.6 | 13.9 | 1.3 | 2.7 | 17.2 | 19.3 | 26.0 |
| snow_depth | 11.4 | 0.0 | 2.2 | 7.0 | 6.6 | 4.4 | 8.5 | 30.2 |
| inch15_F0 | 6.6 | 55.1 | 13.7 | 4.3 | 5.4 | 15.0 | 23.7 | 33.8 |
| snow_cover | 14.0 | 72.8 | 30.4 | 7.1 | 0.0 | 23.4 | 76.0 | 123.1 |
| temperature | 1.3 | 131.5 | 28.2 | 0.8 | 7.2 | 43.8 | 53.0 | 139.3 |
| null_model | 9.8 | 214.0 | 20.1 | 0.0 | 3.7 | 84.1 | 88.8 | 254.8 |
| Index | MI | MN | NH | NY | VT | WI | ME | All |
|---|---|---|---|---|---|---|---|---|
| LSDNT | 11.6 | 15.5 | 26.0 | 0.0 | 18.7 | 4.2 | 14.6 | 0.0 |
| snow_depth | 28.8 | 53.1 | 30.7 | 10.7 | 14.2 | 0.0 | 19.7 | 48.9 |
| inch15_F0 | 13.5 | 26.9 | 0.0 | 3.0 | 17.9 | 94.7 | 0.0 | 191.6 |
| inch18_F0 | 0.0 | 65.1 | 13.6 | 8.0 | 25.2 | 95.3 | 37.8 | 250.7 |
| temperature | 25.7 | 0.0 | 25.3 | 16.1 | 26.4 | 58.7 | 138.2 | 551.0 |
| snow_cover | 18.3 | 65.9 | 70.4 | 86.2 | 68.2 | 0.2 | 141.8 | 569.7 |
| null_model | 82.8 | 110.3 | 32.5 | 24.5 | 0.0 | 84.8 | 162.0 | 767.4 |
Figure 7. Mean WSI values during 2004–2023 calculated using the LSDNT formula, Northeastern USA.
Figure 8. Ad-hoc classification of WSI zones (see Figure 2) using mean WSI values and fitted effects on deer harvest during 2004–2023, Northeastern USA.
If the premise is accepted that over-the-counter adult male deer harvest is strongly correlated to overall deer population performance (i.e., abundance and recruitment) at large spatio-temporal scales, then the predictive models from this analysis can be used to predict the general impact of winter weather on deer populations.
For the predictive maps below, I used the following parameters to predict annual impacts during 2003-2024 from LSDNT WSI values:
| Deer Population Effect Class | Predicted Harvest Change | Map Color |
|---|---|---|
| Strong Negative | < -20% | |
| Moderate Negative | -10% to -20% | |
| Mild Negative | -5% to -10% | |
| Marginal Negative | -2% to -5% | |
| Neutral/Average | -2% to +2% | |
| Marginal Positive | +2% to +5% | |
| Mild Positive | +5% to +10% | |
| Moderate Positive | +10% to +20% | |
| Strong Positive | > +20% |