As a globally widespread apex predator, humans have unprecedented lethal and non-lethal effects on prey populations and ecosystems. Yet compared to non-human predators, little is known about the drivers and consequences of human hunting behavior. Here, we characterized the hunting modes, habitat selection, and harvest success of 483 rifle hunters in California using high-resolution GPS data. We used Hidden Markov Models to characterize fine-scale behavior, and k-means clustering to group hunters by hunting mode, on the basis of their time spent in each behavioral state. Hunters exhibited three distinct and successful hunting modes (“coursing”, “stalking”, and “sit-and-wait”), with stalking as the most successful strategy. Across hunting modes, there was variation in patterns of selection for roads, topography, and habitat cover, with important differences in habitat use of successful and unsuccessful hunters across modes. Our study indicates that hunters can successfully employ a diversity of harvest strategies, and that hunting success is mediated by the interacting effects of hunting mode and landscape features. Such results highlight the breadth of human hunting modes, even within a single hunting technique, and lend insight into the varied ways that humans exert predation pressure on wildlife.

Methods

To understand patterns of hunter movement behavior, we collected GPS tracks of all hunters in the study area (2015-2022). We first classified the behavioral state of each location for each hunter, using Hidden Markov Models (HMMs; location-level classification of behavior). Next, we used k-means clustering to group hunters into distinct hunting modes based on the relative time that they spent in each behavioral state (hunter-level classification of behavioral strategy). Finally, we used Resource Selection Functions (RSFs) to evaluate patterns of habitat selection for each hunting mode, comparing habitat selection between successful and unsuccessful hunters.

Study area

We conducted primary data collection at the 2,168-hectare Hopland Research and Extension Center (HREC) in Mendocino County, California (Latitude: 39.002, Longitude: -123.084; Figure 1). The site features habitat types including grassland, oak woodland, and chaparral, with a network of dirt roads and fences. The site hosts an annual public hunt, in which twenty hunters per day are selected by lottery from a pool of applicants, for 4-6 days each year. In 2020, a restricted multi-day hunt was introduced for a small number of hunters.

Data collection

Our study took place each August-September from 2015-2022, excluding 2018 due to wildfire. We invited all hunters at the study site to participate in our study. We had a 100% rate of participation (n = 483 hunters representing 648 hunter-days).

We provided each hunter with a GPS unit (i-gotU GT-600) that was programmed to take a GPS fix every 5 seconds from 5am to 10pm to encompass legal hunting hours at the study site. We asked hunters to keep the GPS unit in a pocket that would remain on their person, even when they were moving on foot. All harvested deer were brought back to headquarters, and we confirmed with the hunters whether each logger was associated with a successful or unsuccessful hunt. Upon data retrieval, we resampled all tracks to a fix rate of 3 minutes to accommodate GPS error and computational limitations. We followed data cleaning procedures described in detail in the supplementary methods.

Spatial data

We identified environmental features that we a priori hypothesized to influence hunter behavior and habitat use: distance to nearest road, ruggedness, viewshed, and density of each of the three habitat types (woodland, grassland, and chaparral). These hypotheses were drawn from existing literature on human hunter movement and behavior. We generated raster layers for each feature in the study area. Additional details on the development of spatial variables are provided in the supplementary methods section. We extracted spatial covariates at each point, and we calculated the elapsed time since sunrise for each point using the suncalc package. We standardized all covariates prior to modeling.

Behavioral state classification with Hidden Markov Models

To identify fine-scale behaviors of hunters, we used the moveHMM package to fit a hidden Markov model to the hunter movement data. We ran a global model with all predictors (distance to road, viewshed, ruggedness, woodland density, chaparral density, and time since sunrise). We assigned movement points to one of three behavioral states, as initial modeling indicated that three-state models performed better than two-state models (based on AIC), and best corresponded to self-described hunter behavior. We interpreted State 1 as corresponding to a stationary state (searching, resting, or processing deer), State 2 to walking on foot, and State 3 to driving in a vehicle.

We followed best practices when choosing initial parameter values. We included a zero-mass parameter for step length given the high proportion of step lengths equal to 0 (17% of all steps). To determine whether our models were sensitive to initial parameter choice, we ran 100 iterations of the model with randomly-chosen starting parameters for step length mean, step length standard deviation, step length zero mass, and turning angle concentration. Our model converged on the same parameters for 82 of 100 of the iterations, and this model had the maximum likelihood, indicating numerical stability. We then used the parameter values from the best model as our starting values for all subsequent modeling. Based on the global model, we determined the most probable behavioral state at each step for each hunter, and determined the percentage of time that each hunter spent in each of the behavioral states. Following identification of the three states from movement parameters, we further distinguished between resting behavior on road (<10 meters from road) and off road (>10 meters from road), as these behaviors are associated with different hunting strategies (resting on the road to visually scan for deer on the landscape vs. resting off the road in a sit-and-wait hunting strategy).

Identification of distinct hunting modes

To identify the dominant hunting mode of each hunter, we used k-means clustering to group hunters on the basis of their time spent in each fine-scale behavioral state. We determined the optimal value of k using the elbow method heuristic. Specifically, we plotted the total within-cluster sum of squares as a function of k, and determined the value of k at which this sum of squares began declining linearly.

We then ran logistic regressions to evaluate the effect of hunting mode on harvest success. Additional model covariates included year (as we were interested in whether hunting success changed over time) and whether the track came from a single-day or multi-day hunter. We tested all possible covariate combinations and we also explored interactions among hunting mode and the other covariates, to examine whether the effectiveness of different hunting modes changed over time, or varied between single- and multi-day hunts. We compared models using AIC. We also determined relative variable importance (RVI), as calculated by summing the Akaike weights of all models in which the variables appeared.

We also evaluated whether the time of day at which deer were harvested varied across hunting modes, for hunters for which we had known harvest times (n = 37 of 39 successful tracked hunters). We compared harvest time (elapsed time since sunrise) for each of the three clusters using an Anderson-Darling test, a non-parametric rank test of whether samples from different groups came from the same distribution.

Evaluating habitat selection

To evaluate patterns of habitat selection by hunters using different hunting modes, and to evaluate connections between habitat selection and harvest success, we used Resource Selection Functions (RSFs). RSFs compare environmental features of used versus available locations in a logistic regression. We compared locations recorded by hunter GPS trackers (used locations) to locations that we systematically sampled throughout the huntable area at a 30 x 30 m resolution. To evaluate potential links between habitat selection and hunting success, we ran separate models for successful and unsuccessful hunters in each of the hunting modes, and we used the same predictors in all models to facilitate comparison of model coefficients. Model covariates included the same spatial covariates used in the HMM: ruggedness, viewshed, chaparral density, and woodland density. We assigned a weight of 5,000 to the available points, and 1 to the used points. We first ran mixed models with a random intercept for track ID, but among-individual variance was 0 for all models, resulting in a singular fit. We therefore removed the random intercept to ensure estimate stability. In addition, because RSFs assume spatiotemporal independence between points, we checked the effect of fix interval. We thinned the “used” points to a 30 min interval to reduce spatial autocorrelation between points, while retaining sufficient data for each individual hunter (mean of 20.3 used points per hunter). Conclusions remained unchanged despite the 10-fold reduction in fix rate, and the results of this model are presented in the supplementary material.

To rule out any potential issues of circularity when using some of the same spatial covariates to classify behavior (which was then used to identify hunting mode) and to compare differences in habitat selection across hunting modes, we also re-ran the HMMs without any spatial covariates, classifying behavior based only on step length and turn angle. We then re-ran the k-means clustering analysis and RSFs with the updated behavior, and conclusions again remained unchanged. We have chosen to retain the spatial covariates in the HMM for all analyses presented here, to improve the accuracy of behavioral classification (particularly with regard to driving). The HMM uses information about the spatial covariates to improve predictions of behavioral states at a given location, while the RSF compares all used locations (regardless of behavioral state) to a random set of available points within the study area to evaluate patterns of habitat selection.