## Abstract

Predation risk is an important driver of ecosystems, and local spatial variation in risk can have population-level consequences by affecting multiple components of the predation process. I use resource selection and proportional hazard time-to-event modelling to assess the spatial drivers of two key components of risk—the search rate (i.e. aggregative response) and predation efficiency rate (i.e. functional response)—imposed by wolves (*Canis lupus*) in a multi-prey system. In my study area, both components of risk increased according to topographic variation, but anthropogenic features affected only the search rate. Predicted models of the cumulative hazard, or risk of a kill, underlying wolf search paths validated well with broad-scale variation in kill rates, suggesting that spatial hazard models provide a means of scaling up from local heterogeneity in predation risk to population-level dynamics in predator–prey systems. Additionally, I estimated an integrated model of relative spatial predation risk as the product of the search and efficiency rates, combining the distinct contributions of spatial heterogeneity to each component of risk.

## 1. Introduction

Risk of predation alters prey distribution and demography, the effects of which can cascade across trophic levels [1]. As mediated by changes to prey behaviour, spatial heterogeneity in predation risk is thus an important driver of ecosystems [2]. Characterizations of spatial predation risk are often limited to predator density models, where predation risk is indexed with estimates of predator presence, density or habitat suitability [3]. However, spatial heterogeneity may alter broad patterns of predator density (i.e. predator numerical response), redistribute local predator density through effects on predator spatial allocation of search effort (i.e. aggregative response) or alter the rate at which predators successfully kill prey while searching (i.e. functional response). Spatial models of each of these components are needed to explicitly link space and predation to population dynamics, as well as to estimate complete and integrated spatial depictions of predation risk.

At small temporal or spatial scales, the predator numerical response varies relatively little, and the product of predator spatial aggregative and functional responses estimates spatial patterns in prey mortality [4]. Lima & Dill [5] accordingly represented predation risk as a function of the rate of encounter between predator and prey, *α*, and the probability of death given an encounter, *d*. Hebblewhite *et al*. [6] subsequently estimated spatial predation risk, using distinct spatial models of *α* and *d*. In multi-prey systems, spatial drivers of both predator encounter [3] and kill rates [7] can vary by prey species, and thus spatial variation in risk becomes a complex multi-dimensional landscape of species-specific rates of encounter and vulnerability. However, an overall multi-prey spatial predation risk model might also be estimated from the perspective of the predator alone as a function of *a* and *d*, where *a* is the spatial search rate of predators and *d* is the predator's efficiency in making kills of any prey species while searching. With this perspective, the spatial drivers of predator search and efficiency rates include underlying but unmeasured patterns of relative prey density and selection by predators. Scaling predation risk back up to population dynamics, spatial models of *a* and *d*, respectively, characterize the predator aggregative response in units of predator search effort per unit of space and the predator functional response in kills per search time.

Spatial decomposition of predation risk has been achieved previously, using resource selection function (RSF) approaches [6]. While purely spatial RSF approaches can be used to estimate a local density-based search rate, they may fall short for consideration of predation efficiency by lacking a temporal scale. Merrill *et al*. [8] recently suggested time-to-event (or here time-to-kill) proportional hazard (PH) analysis to incorporate the effect of time on predation efficiency by accounting for a non-static underlying probability of a kill occurring. Modelling how the baseline probability of a kill occurring (or the baseline hazard) varies with time is a necessary precursor to assessing the effect of spatial covariates upon efficiency. Furthermore, the time-to-kill-dependent variable is akin to the inverse of the functional response of predators [9].

Here I demonstrate a novel combination of spatial RSF and temporal PH analyses to decompose and subsequently integrate the spatial variation of predation risk in a predator–prey system of conservation relevance. Studies of wolf search rates have previously suggested an increased search rate of wolves with proximity to anthropogenic linear features such as roads, seismic exploration lines and trails [10,11]. Wolf selection for linear features as travel routes may also theoretically increase wolf predation efficiency [12], offering a potential mechanism for studies relating linear features to functional habitat loss [11,13] and decreased population growth rates [14] for a threatened prey species, woodland caribou (*Rangifer tarandus caribou*). However, empirical findings regarding search rates do not necessarily imply increased efficiency in terms of cumulative kill rates. I used this system as a platform for decomposing wolf predation risk with the goals of explicitly testing whether anthropogenic changes to spatial predation risk are manifested in spatial variation in the wolf search rates, efficiency rates or both. I demonstrate how RSF and PH methods can be used to assess spatial drivers of each component of risk, with conservation implications for woodland caribou.

## 2. Methods

### (a) Predation monitoring

I studied wolf predation within the Rocky Mountains of western Alberta and eastern British Columbia, Canada. The study area included both mountain and foothill ecoregions, which differed in topography, vegetation and jurisdictional management practices. For more study area description, see [13]. I compiled data from global positioning system (GPS) telemetry collars (Lotek GPS3300 and GPS4400 models) deployed on 46 wolves from 23 packs between 2000 and 2011. For packs with two or more GPS-collared wolves, I included data from only the single wolf with the most complete GPS data in analyses. I restricted analyses to GPS data collected at fix intervals of up to 2 h [15], which included intervals of 15, 30, 60 and 120 min. I conducted analyses separately by season for summer (16 May–16 October) and winter (17 October–15 May) seasons as defined by DeCesare *et al*. [13].

I modified a rule-based algorithm previously programmed in Python language (Python Software Foundation, Hampton, NH) by Knopff *et al*. [16] to identify GPS location clusters that represented putative wolf kill sites (see the electronic supplementary material, appendix A). For a subset of 12 packs, I used field validation to detect 258 (23%) kill sites out of 1107 clusters visited. I restricted analyses to kills of large mammal prey species, including moose (*Alces alces*), elk (*Cervus elaphus*), woodland caribou, mule deer (*Odocoileus hemionus*), white-tailed deer (*Odocoileus virginianus*), bighorn sheep (*Ovis canadensis*), mountain goat (*Oreamnos americanus*) and beaver (*Castor canadensis*), and I treated GPS clusters at kill sites of ungulate neonates (less than 30 days old; *n* = 4) as non-kills for this analysis. For each cluster, I estimated a set of seven parameters found to be predictive in previous work delineating kill sites with wolf GPS location clusters (see the electronic supplementary material, appendix A). I then built seasonal predictive models of the probability of a given cluster being a kill site, using generalized linear mixed-effects models (GLMMs) with the logit link and a random intercept for each pack [15]. The area under receiver operating characteristic curves (AUC; winter = 0.906, summer = 0.911) indicated excellent discrimination by model, and I applied models to the remaining clusters for all wolves monitored during the study to estimate a total of 1072 kill sites across 23 packs (see the electronic supplementary material, appendix A).

### (b) Spatial drivers of predator searching

I used RSFs to characterize the spatial variability in the wolf search rate by modelling how underlying landscape attributes affect where wolves concentrate their search effort. I estimated distinct aggregative response RSFs for each of summer and winter seasons, and mountain and foothill ecoregions, pooling data across years. I subsampled GPS data to standardize analyses to data collected at 2 h fix attempt intervals, and I limited analyses to searching (or hunting) wolves by removing repeated locations at kill or rest site location clusters, as defined earlier. This restricted RSFs to characterize selection by wolves specifically while they were moving and presumably searching for prey. I treated GPS locations as samples of wolf resource use while searching and drew an equal number of random locations within each pack's 95 per cent fixed kernel home range to sample spatial availability. I quantified the spatial attributes underlying each wolf and random location with a suite of topographic, vegetative and anthropogenic attributes shown to be important in previous wolf studies. These included the slope, topographic position index (an index of curvature where drainages are represented with negative values and ridges with positive values), distance to streams, land cover type (forest, grass/shrubs, alpine vegetation and rock/ice), distances to all roads, other human linear features (seismic lines in the foothills ecoregion and maintained hiking trails in the mountain ecoregion) and forestry cut-blocks (see the electronic supplementary material, appendix B). I estimated RSFs using GLMMs with the logit link and a random intercept for each pack [17]. I built global models, including all variables with weak univariate significance (*p* < 0.25) and low multi-collinearity (|*r|* < 0.5). I then reduced these to final models using a manual backward-stepping procedure until all variables were significant (*p* < 0.05). I used this model selection routine instead of an AIC-based model averaging procedure as it allowed careful observation of coefficient signs, magnitudes and variance inflation in multivariable space, and resulted in single, parsimonious models that were easily compared among seasons and ecoregions. All GLMM analyses were conducted using R v. 2.13.1 [18] with the lme4 package.

### (c) Spatial drivers of predation efficiency

I used parametric proportional hazards (PPHs) time-to-event models [8,19] to assess the spatial factors that increased the instantaneous hazard, or probability of making a kill, along wolf search paths. This quantified predation efficiency as a response variable in units of kills per search time, and modelled spatial variation in efficiency as a spatial hazard function of spatial factors affecting the rate of making kills while searching. I divided telemetry data into hunting sequences where a given hunt began with the first location away from the previous kill site location cluster and continued until the first location at the next kill site cluster. I censored the time spent revisiting previous kill sites during hunts, thus restricting analysis to inter-kill searching time. I measured time in days and set the start of each hunt at time = 0. I included only those packs with 10 or more kills for a given season and censored data for hunts that extended more than 15 days since the previous kill, as they appeared to be the result of non-hunting movement behaviour such as extra-home range exploratory forays [9].

I included annual estimates of pack size and an ordinal variable characterizing the ordered month per season as covariates potentially predictive of time to kill [20]. For winter models, I also included annual winter-averaged anomaly estimates of the North Pacific Oscillation (NPO) index, which is a reliable index of winter severity for this study area [21]. To characterize spatial attributes underlying each step of wolf search paths, I intersected telemetry locations with the same topographic, vegetative and anthropogenic covariates as considered for aggregative response analysis. In addition to considering the distance to each of three anthropogenic disturbance types (roads, seismic lines and trails, and cut-blocks), I also estimated whether wolf locations were on or off these features (see the electronic supplementary material, appendix B). I built separate PPH models for each season and ecoregion, and used shared frailty models to account for autocorrelation within packs [19]. I began each model selection process by using AIC to select among exponential, Weibull and Gompertz parametrizations of the baseline hazard [19]. I used univariate significance (*p* < 0.25) to compile global multivariable models, and I reduced these to final models with a backward-stepping procedure until all variables were weakly statistically significant (*p* < 0.1).

As a form of validation, I tested whether predicted spatial instantaneous hazard models could be successfully scaled up to predict among-pack variation in the observed functional response of kills per unit search time. I predicted the cumulative hazard of each hunting path to combine the effects of space and time into a metric of the total predicted hazard encountered by each hunt from start to finish. I estimated the cumulative hazard at the time of completion of each hunt, *H*(*t*), as the integral of instantaneous hazard across the entire search path,
2.1where *u* is the hazard function across variation in spatial attribute data *x*, model coefficients are represented by *β* and the integral of *h*(*u*)d*u* represents the cumulative baseline hazard, *H*_{0}(*t*). I estimated the cumulative baseline hazard specific to each model's parametric parametrization, using *H*_{0}(*t*) = exp(*β*_{0})*t* [exponential], *H*_{0}(*t*) = exp(*β*_{0})*t ^{p}* [Weibull] or

*H*

_{0}(

*t*) = exp(

*β*

_{0}){exp(

*γt*) – 1}

*γ*

^{−1}[Gompertz], where

*p*and

*γ*are shape parameters for Weibull and Gompertz distributions, respectively [19]. I then estimated a daily rate of the spatial hazard searched per pack-season by summing per-hunt cumulative hazard values within packs and dividing by the total search time during which a kill was possible. These steps allowed me to scale the local spatial hazard model predictions up to a total rate of hazard encountered by each pack while searching. For comparison, I calculated the observed predation efficiency per pack-season in terms of total kills per unit search time. I counted the number of GPS location clusters that were kill sites for each pack-season and divided this by the total search time monitored.

I used linear regression with a clustered sandwich variance estimator (to account for two seasons recorded for some packs) to test whether the increased predicted cumulative hazards underlying wolf hunting paths were indeed associated with increased kill rates per pack. Cumulative hazards for each pack-season were estimated according to season- and ecoregion-specific models, which caused variation among pack seasons in this comparison to be a function of variation in both the space use of the pack and the underlying baseline hazard for a given season and ecoregion. Thus, I also conducted four separate linear regressions for subsets of data within each season and ecoregion combination, wherein differences among packs could be solely attributed to variation in spatial factors. All efficiency analyses were conducted in Stata v. 10 [22].

### (d) Integrated modelling of spatial predation risk

RSFs based on use-availability data allow the estimation of relative probabilities of selection rather than that of true probabilities [23]. Thus, I estimated the search rate, *a*′, as the relative probability of selection using
2.2

These predicted values are akin to the odds ratio formulation of individual model coefficients, wherein the relative differences among predicted values represent the relative probabilities of a wolf searching event, given underlying spatial covariates. PPH models offer a parallel hazard ratio formulation of relative predator efficiency rate as
2.3where *β _{k}* is the PH coefficient for the

*k*th spatial covariate [19]. Hazard ratios are interpreted, similar to odds ratios, as the relative probability of a wolf making a kill while searching, and are assumed to remain constant over time. I used ratio formulations (exp[

*β*]) of search and efficiency models to estimate spatial models of integrated predation risk according to 2.4

_{j}x_{j}## 3. Results

RSFs revealed consistent spatial drivers of the wolf search rate across seasons and ecoregions, including topographic, vegetative and anthropogenic conditions (table 1). Wolves selected gentle slopes, drainage topography close to streams, grassland and shrubland vegetation types, and areas close to roads, seismic lines and trails, and cut-blocks when searching (table 1). However, the primary spatial variable driving predation efficiency in areas searched by wolves was the topographic position index, wherein the relative hazard, or probability of making a kill, increased in drainage topography across all seasons and ecoregions (table 2). Additionally, the hazard increased with winter severity, as indexed by the NPO during the winter season, and varied inconsistently with changes in slope and forest land cover (table 2). Anthropogenic features did not affect predation efficiency, with the exception of the winter season in the mountain ecoregion, where the relative hazard was lower (*β* = −0.227, *p* = 0.043) when wolves were searching within 250 m of seismic lines and trails.

The predicted cumulative hazard rate underlying the GPS paths of wolves correlated with wolf predation efficiency, as measured with kill rates per pack-season (*β* = 1.167, *p* < 0.001, *R*^{2} = 0.693; figure 1*a*). However, when separate regressions were conducted within each combination of the winter (w) and summer (s) seasons and the mountain (m) and foothill (f) ecoregions (wherein variation was solely due to wolf space use by ensuring all packs had the same predicted baseline hazard), this relationship was consistently positive, but statistically significant only for two of the four models (*β*_{wm} = 1.440, *p* = 0.002; *β*_{wf} = 1.176, *p* = 0.222; *β*_{sm} = 0.390, *p* = 0.690; *β*_{sf} = 1.833, *p* = 0.033; figure 1*b–d*). When estimating the overall predation risk as the product of search and efficiency rates, the increase in both wolf search effort and killing efficiency in topographic drainages led to high predation risk in these areas, while the contribution of anthropogenic disturbance to predation risk was largely driven by increasing search effort alone (figure 2).

## 4. Discussion

PH models revealed a previously undocumented link between local spatial variation in the instantaneous hazard, or probability of a kill, and the functional response of predators at broad scales (figure 1). These results highlight the potential of PH models as a tool for scaling local changes in spatial heterogeneity of predation risk up to landscape-level variation in predator kill rates, and for linking behavioural studies of risk to population-level consequences [1–3]. I also used RSFs to model the predator aggregative response, or variation in local predator density due to the space use of a fixed number of predators. Assessment of spatial drivers of both the search and efficiency rates allowed improved understanding of how different spatial factors affect different components of predation risk. It also allowed the estimation of an integrated spatial predation risk model accounting for the effects of spatial heterogeneity on each component of risk (figure 2). Ideally, broad-scale spatial models of predator numerical responses could also be integrated with functional and aggregative response models to convert relative *per capita* rates of predation into population-level modelling of both predator and prey population dynamics.

Predator–prey interactions are complex games, wherein predators search for prey and prey avoid predators [24,25]. While predator-based studies such as this one have revealed spatial variation in multiple stages of predation such as search and kill rates [26], prey-based studies have also revealed differential avoidance and escape tactics as prey-mediated components of spatial risk [7]. Spatial models of predation risk become increasingly difficult to generalize as one acknowledges the potential for behavioural variation in both predators and prey among species [7,27], among individuals within species [28–30] and among behavioural states within an individual [31]. Rather than a model of risk for any given prey species or individual, I estimated a multi-species composite model of risk for the average prey individual killed by wolves in my study area. The net spatial pattern of risk with this treatment inherently includes the interacting behavioural mechanisms of predator selection of prey species [32] and the spatial avoidance behaviour and vulnerabilities of prey [27]. This approach tests for links between spatial heterogeneity and predator responses, though such links are undoubtedly mediated by complex and multi-species behavioural games between both predators and prey [7]. While spatial risk may in fact vary among species, the overall predator–prey dynamics in systems of shared predation, and thus apparent competition, are typically driven by preferred and abundant primary prey species [33]. As such, predation risk measured from the perspective of the predator alone may best characterize patterns of risk most relevant to the greater ecological community.

*Per capita* predation risk may also vary with spatial or temporal changes in local densities of predators and prey [34]. Thus, my treatment of risk also includes an unmeasured effect of multi-species prey densities over space [32]. McPhee *et al*. [9] showed that prey density was less predictive of predator efficiency than other spatial factors, suggesting that the effect of prey density may manifest more in predator searching behaviour than efficiency. In other systems, predator searching behaviour has been linked to prey density [35], prey vulnerability [36,37] and an interaction of both [38]. Thus, prey vulnerability and prey density are likely to mediate one or both of the relationships between spatial heterogeneity and components of predation found here.

The integration of RSF and PH models revealed both consistencies and discrepancies in the relationships between spatial heterogeneity and components of predation. In the case of wolves, topographic drainages proved important in driving both search and efficiency rates. Proximity to anthropogenic linear features was consistently associated with higher wolf search effort across seasons and ecoregions, but did not positively affect the kill rate of wolves while searching. Thus, the previously hypothesized positive relationship between anthropogenic disturbance and wolf predation risk seems to be supported as mediated by the aggregative response alone, and not through a mechanism of predation efficiency. This result offers mechanistic insight to explain previous results showing that despite consistent selection of wolves for linear features as travel routes [10,11], neither wolf kill sites [10] nor caribou mortalities [11] were found significantly closer to such features. Furthermore, this example highlights the strength of considering both predator search and efficiency rates when assessing the effects of risk upon prey behaviour and demography. Rates of overlap or encounter alone do not necessarily forecast demographically meaningful variation in prey survival [7].

Baseline hazards and related kill rates increased during the winter season (figure 1) and with increased winter severity among years (table 2). This corroborates other findings of increased predator kill rates [39] and prey vulnerability [40] in predator–prey systems with snow, yet seemingly contradicts recent studies indicating higher summer kill rates of wolves, specifically [20,41]. I used different seasonal models and cut-off values for predicting kills from GPS location clusters (see the electronic supplementary material, appendix A) and my kill rates are in units of kills per search time rather than kills per total time. Thus, the potential for seasonal differences in the relative allocation of searching and handling time—which should be expected with seasonal differences in the foraging costs for predators or vulnerabilities of prey [42]—may also explain seasonal kill rate differences in this study. Additionally, Sand *et al*. [20] and Metz *et al*. [41] both showed that despite higher overall summer kill rates, wolf predation was lower on older and larger ungulate age classes, which corresponded to my results, wherein I modelled these age classes exclusively by excluding neonates.

Overall, integration of predator search and efficiency models allowed a more complete depiction of local predation risk than estimates of predator density alone (figure 2), and the modelling process revealed empirical relationships between spatial factors and each component of risk. Future efforts in the research of predation risk effects upon prey would serve well to assess how the relative contributions of each component and the integrated product of predation risk ultimately drive prey behaviour, demographic responses and the cascading of risk effects to ecosystem processes.

## Acknowledgments

Wolves were captured during winter using helicopter net-gunning, and during summer with rubber-padded foothold traps (no. 7 EZ Grip Trap; Livestock Protection Company, Alpine, TX), according to approved animal care and use protocols (University of Montana Animal Use Protocol no. 059-09MHWB-122209).

Much of the conceptual foundation and management behind this project was the work of M. Hebblewhite, and data collection was led especially by S. Hazenberg, but also by M. Bradley, G. Kuzyk, L. Neufeld, F. Schmiegelow and J. Whittington. My sincere and humble gratitude goes to these generous collaborators. Please submit data-sharing requests to the author of this paper, and they will be distributed to multiple primary holders of these data. Support for this research programme was provided by the Alberta Conservation Association, Alberta Sport, Recreation, Parks and Wildlife Foundation, Government of Alberta, British Columbia Ministry of the Environment, Canadian Association of Petroleum Producers, Canadian Forest Products, Foothills Research Institute, Montana Institute on Ecosystems, NSERC, Petroleum Technology Alliance of Canada, Parks Canada, Royal Dutch Shell Canada, University of Alberta, University of Calgary, University of Montana, Weyerhaeuser Company and World Wildlife Fund. I thank J. Berger, L. S. Mills, M. Musiani and D. Pletscher for helpful reviews, and C. Carli, M. Colombo, D. Dupont, C. Mamo, W. Peters, B. Weckworth and others for their contributions.

- Received July 20, 2012.
- Accepted August 22, 2012.

- This journal is © 2012 The Royal Society