Efficacy of the 2004 model in Michigan

Treves et al. (2004) developed the following linear equation using a single-sample discriminant function analysis (Morrison 1990), to calculate relative risk (R) of wolf predation on livestock in Wisconsin and Minnesota:

This model used survey townships (White 1983; 92.16 km²) as the spatial resolution and included landscape variables (% area/township) and deer density (deer/km²) as independent variables. The public land survey system refers to a square unit of land, which is nominally six miles (9.7 km) on a side. Each 36 square mile (93 km²) township is divided into 36 one-square mile (2.6 km²) sections, and each section covers exactly 640 acres (2.6 km²).

We combined deer density estimates and the following spatial data sets from the UP to evaluate the 2004 model's initial predictive capability: 1) wolf population range in Michigan since 1996, 2) locations of all 87 verified wolf predations on livestock during 15 April 1996 - 14 April 2009 and 3) 2001 remotely sensed land-cover data. We considered wolf range in Michigan to consist of the mainland UP (about 42,610 km²) based on MDNRE radio-telemetry data (D. Beyer, MDNRE, unpubl. data) and winter snow track surveys (Michigan Department of Natural Resources 2008). We included all verified wolf predation events involving livestock (N = 87) that occurred during 15 April 1996 - 14 April 2009, provided by MDNRE. We considered a verified predation event to consist of ≥ 1 domestic or captive livestock animal killed or injured in a single occasion and cause of death or attack confirmed and attributed to wolves by MDNRE personnel or their agent (i.e. US Department of Agriculture, Wildlife Services). We included records of cattle, sheep, ducks, white-tailed deer, poultry and rabbits deemed livestock by The Michigan Animal Industry Act (Michigan Public Act 466 of 1988). Of livestock predation events, 77% involved cattle, followed by sheep (10%), chickens (7%), ducks (2%), white-tailed deer (2%), pheasants (< 1%), turkeys (< 1%), geese (< 1%) and rabbits (< 1%). We recorded locations to township, range and section and analyzed spatial information at the township scale.

We estimated deer density (deer/km²) for each deer management unit (DMU) in the UP (422 - 4,604 km² in 2009, N = 22) following Treves et al. (2004) and using sex-age-kill model data or deer pellet count surveys (MDNRE, unpubl. data). We used 2001 National Land Cover Data (NLCD, Landsat TM data, 30-m resolution; Vogelmann et al. 2001) to reduce disparity between time of predations (1996 - 2009) and land cover data collection. We measured the percent area of pasture/hayfield, coniferous forest, crops, emergent wetland and open water for each township. We applied the following reclassification scheme from the 2001 NLCD to fit land covers in the 2004 model: we reclassified non-vegetated farmland and forage crops/non-tilled herbaceous agriculture as pasture/hayfield, row crops as crops and pines, other upland conifers, lowland coniferous forest and mixed upland conifers we reclassified as coniferous forest. We did not reclassify open water and emergent wetlands.

We mapped all 87 verified wolf predations using a geographic information system (ArcMap 9.3.1, Redlands, California). We identified affected townships (N = 36) as those in which at least one verified livestock predation occurred and we identified unaffected townships (N = 578) as those where no verified livestock predations occurred. We used a matched-pair design by randomly assigning each affected township an unaffected, contiguous township (Treves et al. 2004). We compared landscape variables between paired affected and unaffected townships. Potentially, eight unaffected townships could neighbour each affected township. We excluded one affected township because the area was < 50% of a standard township (92.16 km²). Our final sample consisted of 35 township pairs (11.4% of townships in the UP, N = 614; Fig. 1).

We measured variables for all townships and applied them to the 2004 linear equation to obtain R values. We used these values to compare predicted classifications of townships to actual classifications (i.e. affected or unaffected) to determine the model efficacy. We determined predicted classification by the proximity of the resulting R value of each township to the classification cut-off value (i.e. mean of group centroids, affected and unaffected; Morrison 1990). Townships predicted as affected had R values greater than the cut-off value, whereas townships predicted as unaffected had R values less than the cut-off value. We considered townships within our sample (N = 70) as being correctly classified if the actual and predicted classifications were congruent. We determined the percentage of affected and unaffected townships that were correctly classified. We averaged the percentages of affected and unaffected townships correctly classified to determine the overall efficacy of the 2004 model for the UP.

Model adaptation

Our initial analysis revealed that, without adaptation to the UP, the 2004 model had a predictive accuracy of 61% overall in the UP (affected and unaffected townships) compared to 76.5% in Minnesota and Wisconsin (Treves et al. 2004). Based on these results we adapted the 2004 model to the UP to maximize the predictive accuracy. We first computed univariate tests of association using one-sample t-tests (affected minus unaffected) and binomial sign tests to assess the discriminatory ability for each variable in the 2004 model using the Michigan data. We accepted statistical significance for all analyses when α < 0.10. The t-test indicated whether the mean difference between landscape features of affected and unaffected township differed from zero. The binomial sign test indicated whether the variable in question discriminated matched pairs of townships with better than the chance probabilities. We retained any variable considered significant in either test for model development.

We performed a single-sample discriminant function analysis using Predictive Analytic Software (PASW Version 18, IBM SPSS, Chicago, Illinois) following Morrison (1990). This linear combination of landscape variables consisted of the sum of products of the variable coefficients with the observations and yielded the discriminant scores (R), which we used to predict risk of wolf predation on livestock. We used R values to form predictive classifications of townships (i.e. affected or unaffected) from our adapted linear discriminant equation. We computed the percentage of original cases (N = 70) correctly classified to determine the efficacy of our adapted model (i.e. actual classification vs predicted classification). We also calculated effect size of R values between affected and unaffected townships (Cohen's d; Cohen 1988). Finally, we computed a leave-one-out cross-validation (LOOCV; Moore & Lee 1994), where we classified each case using a discriminant function based on all the remaining cases.

Risk maps

We used the adapted linear equation to produce two predictive maps of relative risk of livestock predation by wolves in the UP, following Treves et al. (2004). We estimated the relative risk using our adapted model so that the township risk values > 2 SD above the mean of our sample of 35 affected townships were coded as highest risk, those >1 SD above the sample mean were coded moderately-high risk, those within ±1 SD of the sample mean were coded as moderate risk, those > 1 SD below the mean were coded as moderately-low risk and those > 2 SD below the sample mean were coded as lowest risk. We coded townships with 0% pasture as lowest risk regardless of their other landscape attributes. Our first map assumed a uniform distribution of wolves in the UP. Our second map incorporated road density (km/km²) per township, where any township with > 0.70 km/km² of roads was reclassified as lowest risk (unlikely wolf territory establishment; Potvin et al. 2005), and townships with < 0.70 km/km² retained their original risk assignments. We considered the value 0.70 km/km² as an effective indicator of wolf territory establishment in the UP (Potvin et al. 2005). We present both maps because road density has been reported as an effective predictor of wolf territory establishment in the region; however, the importance of this variable remains equivocal (Thiel 1985, Thurber et al. 1994, Mladenoff et al. 1995, Mech 2006, Mladenoff et al. 2006). To generate maps, we extrapolated from the 35 affected townships to all mainland townships in the UP (N = 614). We assumed that the discriminating variables determined from univariate tests were discriminatory across all UP townships and could be applied to predict risk.

Efficacy of the 2004 model in Michigan

The predictive model developed by Treves et al. (2004) correctly discriminated 61.4% of UP townships (i.e. affected 57.1% and unaffected 65.7%) from the universe of townships tested, which was considerably less than the accuracy reported for Wisconsin and Minnesota in Treves et al. (2004). Therefore, we adapted the model based on spatial data derived from the UP to maximize the model's predictive accuracy in the UP.

Model adaptation

To adapt the 2004 model, we used the subset of variables (Table 1) which passed one or both univariate tests of association (P < 0.10) in the UP. These included pasture/hayfield (sign test: P = 0.01, t₃₄ = 3.89, P = < 0.01), crops (sign test: P = 0.04, t₃₄ = 3.40, P = < 0.01) and coniferous forest (sign test: P = 0.27, t₃₄ = 3.15, P = < 0.01). The adapted linear combination of variables (see Equation 2 below) correctly classified 70% (sign test: P = < 0.01, t₃₄ = 4.11) of affected and unaffected townships overall (63 and 77% of affected and unaffected townships, respectively) meeting our overall acceptable threshold:

Effect size analyses revealed that affected townships contained 142.8% more pasture/hayfield, 0.4% more crops and 26.1% less coniferous forest on average than their matched unaffected townships. The standardized difference between mean R values of affected and unaffected townships was large (Cohen's d = 0.84; Cohen 1988). A LOOCV revealed that 64.3% of the affected and unaffected townships overall were correctly classified (60.0% of affected townships and 68.6% of unaffected townships).

Road density served as a predictor of likely wolf territory establishment in the UP in our second map of relative risk. Thus, we computed effect size for road density between affected and unaffected townships as well. Affected townships contained 24.2% more roads (km/km²) than their matched unaffected townships. However, post-hoc univariate tests of association revealed that road density does not discriminate affected from unaffected townships in the UP (sign test: P = 0.18, t₃₄ = 1.47, P = 0.15).

Risk maps

Inspection of land cover data revealed that 44.9% (N = 276) of the townships in the UP contained 0% pasture/hayfield. Within our sample of 70 affected and unaffected townships, our first map (without incorporating road density) identified three (4.2%) as highest risk, one (1.4%) as moderately-high risk, 44 (62.8%) as moderate risk, 15 (21.4%) as moderately-low risk and seven (10.0%) as lowest risk (Fig. 2A). Of the remaining 544 townships, we identified zero (0%) as highest risk, six (1.1%) as moderately-high risk, 148 (27.2%) as moderate risk, 116 (21.3%) as moderately-low risk and 274 (50.3%) as lowest risk. Lowest and moderately-low risk townships were mainly concentrated in the north and central UP. Highest, moderately-high and moderate risk townships were mostly concentrated in the eastern and southcentral regions of the UP.

Relative risk dropped considerably when we excluded areas unlikely to support wolves based on road density (Potvin et al. 2005; > 0.70 km/km²; Fig. 2B). Within our sample of 70 affected and unaffected townships, zero (0%) were identified as highest or moderately-high risk, 18 (25.7%) as moderate risk, 11 (15.7%) as moderately-low risk and 41 (58.6%) as lowest risk. Of the remaining 544 townships, we identified zero (0%) as highest risk, zero (0%) as moderately-high risk, 68 (12.5%) as moderate risk, 88 (16.2%) as moderately-low risk and 388 (71.3%) as lowest risk. Lowest risk townships were mainly concentrated in the north, central and extreme southern UP. Moderate and moderately-low risk townships occurred in the west and across the central UP.

Efficacy of the 2004 model in Michigan

Our initial results showed that the 2004 model did not predict risk of wolf predation on livestock in the UP (61%) as well as in Wisconsin and Minnesota (76%). There are several possible explanations for this result. The 2004 model was constructed only for large livestock (e.g. excluding white-tailed deer, small stock and poultry), which might explain its lower predictive power when applied to the diverse predations in our data set. Similarly, both the 2004 and adapted models are limited by affected and unaffected sites. Townships which were not in our sample and for which we extrapolated relative risk to, may have very different features from the matched pairs used to construct the model. We suggest that model refinement after additional data are collected. Additionally, in contrast to our results, Treves et al. (2004) found that deer density significantly discriminated affected and unaffected townships. Though natural prey abundance influences carnivore predations on livestock (Gunson 1983, Treves et al. 2004, Bradley & Pletscher 2005), we measured deer density in DMUs which are on average, an order of magnitude larger than a township. Because of variation in deer habitat suitability within DMUs, the assigning average DMU deer densities to townships undoubtedly resulted in errors. Improved estimates of deer abundance at a comparable spatial resolution would likely improve our model performance.

Wolves establish territories in areas with greater amounts of open water and emergent wetlands (Mladenoff et al. 1995, Kohn et al. 1996), as these habitats are important as water sources, prey habitat and pup rearing sites (Kohn et al. 1996). However, Treves et al. (2004) observed that wolves preyed on livestock in townships with less emergent wetland and open water. We found no association of risk and open water or emergent wetland in the UP. Differences in association of risk of livestock predation by wolves with open water and emergent wetland within the GLR may reflect the diversity of habitat types in which livestock farms are located.

Model adaptation

Accuracy of our adapted 2004 model (70%) exceeded that of the original 2004 model (61%), suggesting that our revised model is a better predictor of wolf predation on livestock in the UP. However, a LOOCV did not successfully predict ≥ 70% of townships in our sample. The adapted model was better at identifying townships where conflict is less likely to occur (77%) than in townships where it is more likely (63%) in our known sample, when the latter is perhaps a more useful management application.

Still, wolf predations on livestock in the UP occurred in townships sharing a consistent series of landscape variables. For example, our results concur with suggestions that pasture/hayfield is an important predictor of wolf predations on livestock in the GLR (Mech et al. 2000, Treves et al. 2004). Comparable to observations in Wisconsin and Minnesota (Treves et al. 2004), wolf predations on livestock in the UP were positively associated with townships containing higher proportions of pasture/hayfield. We expected this result, as livestock are associated with pasture/hayfield (Bradley & Pletscher 2005). Townships with larger amounts of pasture/hayfield likely contain more cattle, which on average would increase livestock exposure to wolf predation (Mech et al. 2000).

Adaptation of the 2004 model in the UP identified croplands as only slightly positively associated with higher risk across townships, which contradicts results from Wisconsin and Minnesota (Treves et al. 2004). In the UP, there are approximately 900 livestock farms, which are generally located in geographical clusters because of soil and climatic conditions (MDNRE, unpubl. report). Farms containing croplands are likely located in the same geographical clusters. Percent cropland was significantly correlated with percent pasture/hayfield among the affected townships in the UP (r² = 0.61, P = 0.0001, N = 35). Soil and climatic conditions may impose stronger constraints on suitable agricultural areas in the UP compared to Wisconsin and Minnesota. Thus, the predictive ability of croplands may in part be a consequence of geographic association with pasture/hayfield and a reflection of regional variation in the distribution of those two land covers across the Midwest, rather than a statistical artifact. Finally, more data on fine-scale livestock abundance and distribution are needed for Wisconsin, Minnesota and Michigan's UP to more accurately separate the land covers from their uses.

Our results, as well as those in Treves et al. (2004), suggest that coniferous forest is negatively associated with predation risk at the township scale in the GLR. However, managers must be aware that differences in predation risk might occur at varying spatial scales (e.g. farm vs township). For example, in contrast to what they observed at the township scale, Treves et al. (2004) found that an increase in coniferous forest was associated with higher livestock predation risk at the farm scale. Coniferous forest in close proximity to the farms probably allows wolves to approach livestock closely without detection similar to hunting patterns described for wild prey (Kunkel & Pletscher 2001).

Managers in the Upper Michigan can use our results to identify townships with higher relative risk of livestock predation by wolves. Our adapted model identified important predictor variables, successfully distinguished affected from unaffected townships and improved on the 2004 model for the UP. Managers can target areas of greater risk for outreach, education and non-lethal control techniques to manage wolf-livestock conflicts. Our results can also augment preventative efforts, management plans, indemnification programs and research regarding wolf predations on livestock in Upper Michigan.

Risk maps

Based on adaptation of the 2004 model to the UP, we developed two maps of relative risk of livestock predation by wolves in the UP, which managers can use to identify townships potentially at risk. We created the first map (see Fig. 2A) by assuming a uniform distribution of wolves across the UP. The second map (see Fig. 2B) incorporated road density (Potvin et al. 2005) as an indicator of likely wolf territory establishment and thus focuses consideration of relative risk to those townships most likely to contain wolves (30.1% of all UP townships).

There are certain limitations with both risk maps and managers should consider both maps when planning management actions. For example, neither map differentiated townships with multiple predation events. Our first map (see Fig. 2A) may overestimate the risk, as it occasionally assigned a relatively higher risk to townships that likely contain a relatively lower number of livestock operations due to soil and climatic conditions (e.g. the Keweenaw Peninsula). In contrast, our second map (see Fig. 2B) may underestimate the risk due to its reliance on road density as an indicator of likely wolf territory establishment. Although road density appears to be a useful predictor of wolf occupancy (Thiel 1985, Fuller et al. 1992, Mladenoff et al. 1995, Wydeven et al. 2001, Potvin et al. 2005), this variable may become less important as wolf populations increase (Mech 1995). For example, on 3% of the occasions, our second map predicted the lowest risk assignment for previously affected townships. An extreme example of this was a township that had eight previous predation events on three separate farms from 1996 to 2009, with a road density of 0.97 km/km².

The successful recovery of wolves in the GLR, and their coexistence with humans depends on the ability of people to limit wolf-livestock conflicts (Beyer et al. 2006, Ruid et al. 2009). We anticipate an increase in risk of livestock predation by wolves as populations of humans and wolves increase in the UP. Application and further refinement of this predictive model will aid managers in reducing wolf predation on livestock and help ensure wolf persistence in Michigan.