We assessed patterns in timing, location, and age of grizzly bears involved in human-wildlife conflict, represented by conflict kills of bears and attacks on humans in British Columbia (BC), Canada. We modeled the association between annual variation in ecological variables and human-bear conflict frequency. Whereas we were primarily interested in inter-annual patterns for assessing the relative support for our three hypotheses, we included spatial variables in our model to account for ecological differences among populations.

All analyses were performed at the scale of Grizzly Bear Population Unit (hereafter ‘population’), which is designated by the provincial government of BC and is thought to correspond with geographically and genetically relevant sub-populations (50; Supplementary Methods). The ‘habitable area’ (area excluding glaciers and water bodies) of these populations ranges from 2,698 km2 to 49,268 km2 (mean = 13,316 km2; Supplementary Fig. 4).

Conflict kills

We compiled the annual number of conflict-killed grizzly bears from the ‘Compulsory Inspection Database’ (hereafter ‘CID’), which contains the date, location, and cause (e.g. hunt, conflict kill, road accident) of all known human-caused grizzly bear deaths in BC, from 1977–201450; Supplementary Methods; Supplementary Fig. 5). We also extracted age estimates from the database, which were available for 75% of human-caused kills. We used entries only from 1980-onwards because data quality improved considerably after this point (T. Hamilton pers. comm.). We excluded conflict kills from areas where grizzly bears are considered extirpated or threatened (n = 37; Supplementary Methods) because such kills are anomalies, whereas we were interested in generalizable patterns. We included only late-season conflict (occurring from July 1st onwards; 79% of recorded conflict kills; Fig. 2) as a response in our models, because abundance of spawning salmon (a predictor in our model) and bear predation on salmon peak from summer through autumn51. We used the ArcGIS52 spatial analyst kernel density estimator, which uses a quadratic kernel function as described in53 to visualize the spatial density (number per km2) of conflict kills.

Attacks

We combined a database of all known dates and locations of attacks in BC from 1960 to 1997 provided by Stephen Herrero54 with a database of attacks from 1998 to 2014 provided by the BC Ministry of Environment. We included only ‘severe’ attacks for both time periods (for which details on hospitalization differed: pre-1998 severe attacks included fatalities, and injuries requiring >24 hours of hospitalization, whereas from 1998 onwards, when data on hospitalization duration were absent, severe attacks included fatalities, and injuries of a severity requiring hospitalization [e.g. dismemberment and broken bones]) when examining trends through time. We did this because all severe attacks are recorded by the provincial government, whereas the proportion of ‘minor’ attacks (those not requiring medical attention) recorded is unknown (but see Supplementary Fig. 6 for timing of all recorded attacks).

Spatial correlates

We accounted for geographic variation in climate (temperature and precipitation), grizzly bear and human population densities, and presence or absence of spawning salmon. To assess climatic differences among populations we used the program ClimateBC55, which downscales climatic variables obtained from weather stations across the province to an 800 m × 800 m resolution with high accuracy (R2 >> 0.9 between most predicted and weather-station measurements55,56). We created a 4 km × 4 km grid of points across each population’s habitable area (Supplementary Fig. 4) and calculated the log-transformed mean of climate normals (mean spring and summer temperatures and log total spring and summer precipitation from 1981–2010) extracted at each point. To assess differences in grizzly bear densities among populations, we divided 2012 population estimates50 by habitable area in each population and log-transformed the quotient. We used the 2011 Canadian census57 for spatial assessments of human densities. We attributed human population counts from the finest spatial scale available (census subdivisions) to each bear population unit based on percent overlap of the two spatial scales, divided by habitable area of the bear population unit. For a bear population unit with k subdivisions the calculation was:

We used a province-wide database of spawning salmon enumerations58 to attribute presence/absence of spawning salmon to each population.

Annual correlates

Within each population, we assessed inter-annual variation in number of recent conflict kills, number of recent hunting kills, mean spring and summer temperature and precipitation, and annual salmon availability to evaluate the relative support for our three hypotheses (problem individuals, regional population saturation, food supply). We used climate as a coarse but broadly applicable proxy for terrestrial bear food availability because estimates of terrestrial food across BC do not exist, though climate has been linked to food availability and human-wildlife conflict elsewhere59. Specifically, given their broad association with net productivity60, we used measures of temperature and precipitation, during the growing season (spring and summer) of vegetative grizzly bear foods, including shoots, sedges, and berries31,40,61,62. We calculated annual values of mean spring and summer temperature and log-transformed total precipitation from ClimateBC values extracted from a 4 km × 4 km grid of points across each population. We calculated the number of hunting and conflict kills in recent years using a 3-year rolling window, e.g.:

Within each population, recent hunt and conflict kills were scaled by 2 standard deviations with the mean subtracted, providing a measure of inter-annual variation scaled to the magnitude and variability of these measures in each population. We did not include human population as an annual predictor because such data do not exist annually, and there was little temporal variation among censuses in this period. Similarly, we did not include inter-annual variation in grizzly bear densities because such data do not exist in BC26. We assessed inter-annual variation in spawning salmon biomass across the province from 1980 to 2013, using the Fisheries and Oceans Canada nuSEDS database (Supplementary Methods58; see63,64 for caveats). At each salmon count location (Supplementary Fig. 7), we calculated the annual biomass of each species individually, omitting species-stream time series with data missing for more than eight years total, or for three or more consecutive years. We estimated missing counts in the remaining time series by multiple imputation with a Ricker-logistic model fitted to each stream and species (Supplementary Methods; Supplementary Fig 8). We attributed each stream salmon count location to a bear population unit and calculated the total annual salmon biomass for each population as the geometric mean of annual stream biomasses of all salmon species combined (Supplementary Methods). Within each bear population, annual biomass was scaled by 2 standard deviations with the mean subtracted, providing a measure of inter-annual variation scaled to the abundance and variability of salmon in each population.

Analyses

We assessed associations between ecological correlates and conflict patterns using the R65 package glmmADMB, which estimates parameters by maximizing likelihood66. We used a hierarchical modeling approach combining variables that vary spatially among populations with those that vary temporally within each population (Equation 3). We centred and scaled all predictors (subtracted the mean from each observation and divided by 2 standard deviations) to facilitate meaningful comparisons of effect sizes among predictors67. We ran ‘full region’ models that included all populations, and ‘salmon areas’ models restricted to populations with estimated salmon availability. We visually inspected residuals plotted against each predictor and the fitted values and did not detect any remaining strong patterns. Similarly, we visually assessed autocorrelation in residuals and observed little spatial autocorrelation, and substantial temporal autocorrelation in only one population (excluding this population had no qualitative effect on our results), so we did not include autocorrelation terms for model simplicity.

We modelled the number of conflict-killed bears (y) per unit area (km2) in year i and population unit j as

where μ i[j] and ϕ represent the mean and size parameters of the ‘NB2’ parameterization of the negative binomial distribution68 in glmmADBM, which was used because our data were over-dispersed (ϕ in fitted model estimated as 0.72 for salmon areas and 0.62 for full region model); X i[j] represents a vector of annual predictors (salmon biomass, recent conflict kills, recent hunt kills, mean spring and summer temperature, total spring and summer precipitation) with associated β annual coefficients; and Z j represents a vector of spatial predictors (grizzly population density, human population density, mean spring and summer temperature, mean total spring and summer precipitation, salmon present [yes/no]) with associated β spatial coefficients. The α j term represents random deviations from the overall intercept α that vary with population unit and have variance . The term is an offset term representing the habitable area of each population.

We used an information theoretic approach69, whereby we assessed relative variable importance to weigh support for our hypotheses, and conducted model averaging across two sets of candidate models (‘full region’ and ‘salmon areas’; Supplementary Table 2). Spatial predictors were used to account for spatial differences in conflict patterns and allow for inference about our temporal hypotheses. We assessed support for our three hypotheses using coefficient estimates and variable importance of temporal predictors in the two model-averaged models. For the problem individuals hypothesis we assessed previous 3 years of conflict kills, for the regional population saturation hypothesis we assessed previous 3 years of hunting kills, and for the food supply hypothesis we assessed marine-derived food using salmon availability (salmon areas model only), and terrestrial food using climatic variables (annual precipitation and mean temperature; both models) combined. All analyses were performed using R 3.1.265 (for code and source data see github.com/kartelle/ecology-of-conflict/).

In describing overall patterns of conflict, we used all data available (spanning 1960 to 2014), whereas we were limited to data from 1980 to 2013 for modeling purposes because the salmon database and climate data we used did not extend beyond 2013, and reliable conflict-killed bear data were only available from 1980-onwards.