Predicting Occupancy for Pygmy Rabbits in Wyoming: An Independent Evaluation of Two Species Distribution Models

Species distribution models are an important component of effective natural resource management and conservation planning efforts. Species distribution models are most often developed by relating observation data (presence or presence–absence) for a focal species with environmental predictor variables thought to best represent the environmental conditions required by that species. Selection of predictor variables is informed by professional knowledge and field-based experience with the system in which the focal species interacts, as well as through published investigations into species–habitat relationships; and selection often includes climatic, topographic, edaphic, floristic, and habitat structure variables (Franklin 2009). Actual models may then be structured upon quantitative statistical models (e.g., Phillips et al. 2006; Fitzpatrick et al. 2007; Jensen et al. 2008), empirical models (e.g., Estes et al. 2013), or models based on deductive reasoning and expert opinion (e.g., Yamada et al. 2003; Kiesecker et al. 2009). Species distribution models can then be displayed using Geographic Information System (GIS) software, to digitally map the suite of environmental conditions associated with occupancy throughout a predefined geographic area and resulting in a spatially explicit prediction of species distribution, or more appropriately, of its fundamental niche (Guisan and Zimmermann 2000; Franklin 2009).

Knowing the predictive accuracy of a species distribution model is important to resource managers, who may designate critical habitat, conduct reintroductions, or restrict certain activities in an area based on mapped distributions for species of high conservation concern. Problems associated with validating the predictive accuracy of species distribution models have received considerable attention over the past 20-plus y (e.g., Verbyla and Litvaitis 1989; Rykiel 1996; Fielding and Bell 1997; Barry and Elith 2006), and it is generally advised that models developed for the purpose of predicting species distributions be validated using independent data (Chatfield 1995; Araújo and Guisan 2006; Barry and Elith 2006). Although predictive species-distribution models are often assessed by cross-validating with an independent sample of the original data (reviewed by Franklin 2009), this method can provide overly optimistic estimates of a model's predictive power because of spatial and temporal autocorrelation inherent in the training and evaluation data (Araújo et al. 2005; Araújo and Guisan 2006; Randin et al. 2006).

Two pygmy rabbit Brachylagus idahoensis species-distribution models have been developed in the state of Wyoming, and both are intended to guide landscape-scale conservation management actions. The first model was developed by the Wyoming Natural Diversity Database (WYNDD) as part of efforts to map Wyoming's Species of Greatest Conservation Need (WGFD 2010), and the second model was developed by The Nature Conservancy (TNC) for use in landscape-scale mitigation planning (Kiesecker et al. 2009). Our primary objective was to use current, spatially comprehensive and independently collected pygmy rabbit survey data to evaluate the predictive accuracy of both models across the potential range of habitat on public lands in Wyoming. We viewed the models and their associated maps as representative of the fundamental niche of pygmy rabbits in Wyoming rather than of the actual distribution of pygmy rabbits (i.e., the realized niche). Also, because both models have binary interpretations, we hypothesized that areas predicted to be habitat would be occupied more often than randomness would predict, nonhabitat would be occupied less often than random, and that models that characterized the fundamental niche well would have higher accuracy assessments than those that did not. We then compared the models' footprints to identify regions of congruence and disagreement.

Anthropogenic activities have caused distinct disturbance patterns to be present in most, if not all landscapes on earth (e.g., Antrop 2005; Drummond and Loveland 2010), yet these patterns are not considered in many distribution models. Southwestern Wyoming is undergoing rapid and widespread oil- and gas-energy development, and although anthropogenic disturbances of this type have high potential to influence distributions of small vertebrates (e.g., LoBue and Darnell 1959; Merriam 1995), oil- and gas-field development was not considered in either model. We therefore assessed the specific influence of gas-field-related disturbance on each model's predictive success to ascertain whether energy development influenced the predictive ability of either model.

The pygmy rabbit is a Species of Greatest Conservation Need in Wyoming (WGFD 2010) and was recently considered for federal listing protection with a finding of “not warranted” pursuant to the U.S. Endangered Species Act (Figure 1; USFWS 2010). The species is present in six counties in Wyoming (Carbon, Fremont, Lincoln, Sublette, Sweetwater, and Uinta). However, there remains an imperfect understanding of the species' distribution and habitat associations (Purcell 2006). Pygmy rabbits rely year-round on sagebrush for both food and cover and are patchily distributed throughout their occupied range. Pygmy rabbits are associated with relatively tall, dense stands of big sagebrush Artemisia tridentata occurring on deep, friable soils (Green and Flinders 1980; Katzner and Parker 1997; Purcell 2006); habitats meeting these criteria are also patchily distributed. Maximum sagebrush heights in occupied patches range from 0.5 to 2.1 m, while canopy cover immediately surrounding burrow systems generally exceeds 30% (reviewed in U.S. Federal Register 2010 [USFWS 2010]).

Our study area spanned a six-county area in southwestern Wyoming bounded on the east by a longitude approximately equivalent to the town of Casper, westward to the Utah state line, and from the Colorado state line northward to a latitude approximately equivalent to Dubois, Wyoming (Figure 2). Elevation of survey plots ranged from 1,793 to 2,787 m. Big sagebrush, which was the dominant vegetation species, has three codominant subspecies (A. t. wyomingensis, tridentata, and pauciflora). Saltbush Atriplex argentea and greasewood Sarcobatus nees were present in alkaline areas, and other associated vegetation such as western wheatgrass Pascopyrum smithii, needle-and-thread grass Hesperostipa comata, Sandberg bluegrass Poa secunda, prickly pear cactus Opuntia spp., scarlet globemallow Sphaeralcea coccinea, and rabbitbrush Chrysothamnus spp. occur throughout the area (Knight 1994).

Big sagebrush typically occurs on deep, well-drained soils (Thatcher 1959; West 1988), but a combination of soil moisture, depth, texture, salinity, and local climate variation all influence the fine-scale distribution and structure of big sagebrush, resulting in a distribution ranging from broad coverage in more mesic basins to patchy distributions in drier areas (Knight 1994). The geology of our study area is primarily Tertiary shales and sandstones. Soils on the Tertiary bedrock are an association of Haplocambids and Torriorthents, with Fluvents along ephemeral channels, Mollisols on favorable sites, Psamments on stabilized sand dunes and salinized soils in playas. Sodium-affected soils (Natrargids) occur on alluvial fans on high-sodium parent materials. Uplifted areas of cretaceous and older rock add to the complexity of the area (Munn and Arneson 1998).

Based on climate data for the 30-y period 1971–2000, averaged across each of our survey plots, average minimum and maximum temperatures occurred in January (−9.6°C) and July (17.7°C), respectively (PRISM Climate Group 2006). Peak average monthly precipitation (4.2 cm) occurred during May for 96% of the plots. The average annual precipitation total across all of the plots was 31.8 cm (PRISM Climate Group 2006).

Here, we briefly summarize how each model was originally developed. Model development occurred prior to, and independently of, the evaluation work addressed in this paper. The WYNDD model (Figure 3) was developed using data drawn from WYNDD's central species-occurrence database, which is the principal repository for sensitive-species observation records in Wyoming. Occurrence data included all observations of pygmy rabbits reported as a result of surveys conducted between 1981 and 2007, for which location precision was within 700 m. For observations meeting these criteria, model developers used the most precise location information that was provided. To avoid issues of local sampling bias, (i.e., clusters of points resulting in intensive local sampling efforts) a minimum separation distance of 3,200 m between occurrences was enforced by randomly choosing a single point from each cluster. This resulted in a final occurrence data set of 164 locations, of which 80% (n  =  132) were used to build the model and 20% (n  =  32) were withheld to test the model. Eighteen environmental predictor variables (Table 1) were selected for modeling based on published pygmy rabbit habitat preferences (e.g., Ulmschneider et al. 2004) and exploratory analyses of pygmy rabbit occurrences. Pygmy rabbit occurrence data were contrasted with 10,000 random background points using the default parameterization of Program MaxEnt (Phillips et al. 2006). Model performance was assessed based on test data using several metrics, including Area Under the Curve of Receiver Operating Characteristic Plots and confusion matrices of binary model output. Model resolution was 60 m.

The model developed by TNC (Kiesecker et al. 2009; Figure 4) was deductive, incorporating expert-based pygmy rabbit habitat preferences identified from published information; this information was then used to establish a range of values associated with suitable conditions for each of six environmental predictors. These raster data were combined in a GIS to produce a binary map that explicitly identified locations of pygmy rabbit habitat based on the landscape characteristics of the modeled region (Kiesecker et al. 2009). The model developers used a focal mean function to smooth the result and produce the final 30-m-resolution model. The final model was subject to expert review and validated using independent survey data (provided by J. Dahlke, Wyoming Wildlife Consultants, personal communication), which were collected on and adjacent to the Jonah gas field (Table 2). In the study reported here, we validated this model over a much broader geographic area.

We defined our sampling frame by first considering the geographic extent of predicted habitat in both models. We then buffered the area of predicted habitat in the WYNDD model by 5 km because doing so resulted in a more even distribution of area of habitat and nonhabitat in both models than did any other method we considered. We further restricted our sampling domain to federal or state managed land parcels >0.4 km². Within these constraints, we then selected 100 random survey locations in areas predicted to contain pygmy rabbits and 100 random survey locations in areas predicted to not contain pygmy rabbits. To ensure spatial independence, we imposed a minimum distance of 2 km between all survey plots, and these plots were not allowed to intersect gas wells or well pads. We conducted all GIS analysis using Environmental Systems Research Institute ArcGIS Desktop version 10.0 (ESRI 2010).

During 2008 and 2009, we conducted pygmy rabbit occupancy surveys at 187 (we could not access 13 of the original 200 plots) survey plots (Figure 2). At each plot, we surveyed a 400 m × 400 m area by walking parallel transects spaced at 50-m intervals while guided by handheld GPS units. Surveyors deviated up to 25 m from each transect to thoroughly search all patches of potential pygmy rabbit habitat, and searched all potential habitat patches within 50 m of the external plot boundary. Pygmy rabbit sign consisted of observations, pellets, or burrows, and all surveyors were trained by an individual with prior pygmy rabbit survey experience. Pellets and burrows were screened to determine whether they met a set of criteria specific to pygmy rabbits (UDWR 2003; Ulmschneider et al. 2004; MNHP 2008) and all confirmed sign and corresponding Universal Transverse Mercator location coordinates were noted on field forms. In most cases where pellets were observed, a sample was collected and filed at the U.S. Geological Survey Fort Collins Science Center. Once sign was found, surveyors searched the remainder of the surrounding vegetation patch for additional sign and the survey concluded. A plot was considered occupied when a pygmy rabbit, pygmy rabbit pellets, or an active pygmy rabbit burrow was found. Plots were searched until completed or until occupancy was confirmed.

Another component of interest is the effect that anthropogenic disturbance may have on the predictive ability of each model. To evaluate this, we surveyed 120 additional plots on four major gas fields. We randomly selected 40 plots within each of three strata representing low, moderate, and high road densities (<5.0, 5.0–10.0, >10.0 km per km² area, respectively) on the Continental Divide/Creston, Jonah and Pinedale Anticline Project Area (pooled together), and Moxa Arch gas fields (Figure 5). Site-occupancy surveys were conducted during July–November of 2012, using the survey procedures described above.

At each occupied plot, the coordinates associated with pygmy rabbit sign defined the location of the occupied point. At unoccupied plots, the centroid of the plot was used. We compared both models with the field survey data by overlaying evaluation points onto each predictive map and buffering by three distances (75-, 150-, and 300-m radii). Within each resulting polygon, we calculated the percent of cells classified as habitat. We also evaluated the predicted values of the models directly underlying each evaluation point (Data Files S1 and S2, Supplemental Material). As described above, we used a buffered version of the WYNDD model to guide distribution of survey plots. We reviewed the distribution of plots among predicted classes for each original model to verify that neither would be biased by uneven distribution of plots (Table 3).

We used Yule's Q measure of association to identify the combination of buffer size (i.e., search radius) and percentage of predicted habitat contained therein that optimized the predictive performance of each model. To do so, we evaluated all combinations of the four buffer radii mentioned above and minimum predicted habitat levels of 1, 10, and 25%. Yule's Q measures the proportional reduction in error contributed by each combination of variables (i.e., buffer radius and percent predicted habitat), and evaluates how each model performed relative to a random assignment of cases to each outcome class (Wilkinson et al. 1996). After identifying the optimal combination for each model, we assessed the strength of the association between each optimally parameterized model and the evaluation data using Yule's Q test for significance (Wilkinson et al. 1996; Systat 2007).

To assess agreement between the two models, we generated a confusion matrix. Prior to summarizing levels of agreement, the WYNDD model was downsampled (i.e., each 60-m raster cell was divided into four 30-m raster cells) and reprojected into the same coordinate system as the TNC model. “Snapping” was applied in the GIS processing such that cell values within the raster grids were properly aligned to minimize error introduction. Both models were then clipped to the region with predicted scores from both models and a direct “at cell” comparison was made at each predicted pixel in the raster layers.

To evaluate how the effects of anthropogenic disturbance might affect model predictive performance, we compared pygmy rabbit occupancy rates in areas predicted to be habitat but having one of three increasing levels of road density present (Data Files S3 and S4, Supplemental Material). We did so using Goodman–Kruskal's Gamma tests (Bowman and Shetty 2007) and an H_o of no differences across levels of disturbance.

Of the 187 plots surveyed for our model evaluation analysis, 91(48.7%) were occupied by pygmy rabbits. There were no occupied plots in the extreme eastern or northern ends of our study area, suggesting that both models were generally successful at bounding the occupied range of pygmy rabbits in Wyoming. Throughout the remainder of the surveyed area, we saw no obvious pattern in the distribution of occupied and unoccupied plots—both occurred throughout the sampled area.

Although the two models performed similarly overall, they differed markedly in their interpretation of the spatial distribution of habitat at both the local scale and across the predicted range of the models. For the TNC model at all three buffer sizes, a large proportion of plots contained no habitat, and predicted habitat was evenly distributed across the range of possible values in the remaining plots, excepting the smallest buffer radius where approximately 15% of the plots contained 100% predicted habitat. The WYNDD model was bimodally distributed at all three buffer sizes, with a large number of plots containing no habitat, a large number of plots containing 100% predicted habitat, and relatively few plots containing intermediate amounts.

For the TNC model, the optimal combination of buffer radius and minimum percent area classified as habitat was 150 m and 1%, respectively, and this combination of values reduced error in the model by >10% more than any other combination we evaluated. The optimal combination of values for the WYNDD model included the value that directly overlaid the evaluation point, and this predictor reduced 32.7% more error than the next-best combination of values (Figure 6). The TNC model correctly classified 39 of 77 (50.6%) unoccupied evaluation plots and 65 of 88 (73.9%) occupied plots for an overall classification success of 63.3% (Yule's Q  =  0.487, Z  =  3.839, P < 0.001). The WYNDD model correctly classified 53 of 95 (55.8%) unoccupied plots and 59 of 88 (67.0%) occupied plots for an overall classification success of 61.2% (Yule's Q  =  0.439, Z  =  3.550, P < 0.001). Based on 95% asymptotic confidence intervals bounding the Yule's Q estimates, classification success of the two models did not differ. Sensitivity (no. true positives / [no. true positives + false negatives]) of the TNC model was 0.739 and of the WYNDD model was 0.670 (Table 4).

At the landscape scale, the TNC model classified 24.6% of our study area as pygmy rabbit habitat and the WYNDD model classified 38.8% of the area as habitat. The models jointly classified 10.8% of the area as habitat and 47.4% of the area as nonhabitat, but were discordant in classifying the remaining 41.9% of the area (Table 5; Figure 7). The primary reason for the TNC nonhabitat–WYNDD habitat discordance was the expanse of area where the existing vegetation cover type did not fall within the set of classes defined in the TNC model criteria. Alternatively, the TNC habitat–WYNDD nonhabitat discordance occurred primarily where the TNC model predicted habitat in areas beyond the periphery of the region identified as habitat by the WYNDD model. Last, classification success declined sharply for both models as road-density level increased (TNC model Goodman–Kruskal Gamma  =  −0.627, P < 0.001; WYNDD model Goodman–Kruskal Gamma  =  −0.626, P < 0.001; Table 6). For both models, the shape of the function that best fit the data was linear, suggesting a steady decline in predictive success as road density increased.

Species distribution models have been used effectively in terrestrial, freshwater, and marine environments (Elith and Leathwick 2009). Because of their potential to influence policy and management decisions, it is increasingly important that the accuracy of species distribution models be thoroughly validated (Ottaviani et al. 2004; Araújo et al. 2005; Sinclair et al. 2010). Our evaluation of two models designed to predict the distribution of habitat of pygmy rabbits on public lands in Wyoming found that both models were effective at predicting habitat, but ineffective at predicting nonhabitat (based on classification success scores). Although this resulted in Yule's values of overall association that were not particularly strong, the reduction in error achieved by both models was significantly better than random. Neither model was designed to account for the effects of anthropogenic disturbance, and disturbance influenced the performance of both models by increasing commission error in disturbed areas. Key differences existed between models when mapping the geographic distribution of habitat at both local and landscape spatial scales.

Factors that commonly keep species distribution models from achieving high levels of classification success include missing covariates, small sample sizes, lack of absence data, model misspecification, and the presence of unsaturated suitable habitat (Fielding and Bell 1997; Araújo and Guisan 2006; Barry and Elith 2006); and several challenges existed in developing the TNC and WYNDD models. Pygmy rabbits in Wyoming (and elsewhere) are known to associate closely with stands of relatively tall and dense big sagebrush on deep friable soils (Katzner and Parker 1997; Keinath and McGee 2004; Purcell 2006), but the requisite spatial distribution of habitable vegetation patches throughout the landscape is poorly understood. Even were it known, representative, high-resolution digital vegetation data were not available at the time of model development, and high-resolution digital soil-texture maps still do not exist, but represent an important component of pygmy rabbit habitat. Also, population-level effects resulting from predation, disease, habitat isolation, and competition (for burrows and other resources) likely resulted in the presence of unsaturated habitat, but data do not exist with which to model these factors. Collectively, these factors probably impacted the performance of both models.

Other intractable factors may subtly impact post hoc classification assessments of species distribution models if not carefully dealt with. The data set used by WYNDD to develop their model contained only observations of live pygmy rabbits, while our evaluation data consisted primarily of pygmy rabbit sign. Also, the WYNDD data set originated from a large number of observers over a 25-y period. We dealt with this through stringent filtering of the historical data, and thorough training and close supervision of the evaluation data collection crew. However, we acknowledge that error in the WYNDD source data could contribute to modeling error, and error in the evaluation data could affect the reported performance of either model.

The models differed in their interpretation of the distribution of habitat, and this had implications for model interpretation. Based on our local-scale sensitivity analysis, choosing an optimal buffer size and threshold level of predicted habitat with which to evaluate each model was nontrivial; for example, had we arbitrarily used the optima for either model, the other model would have performed poorly in that assessment. At the landscape scale, the TNC model interpreted pygmy rabbit habitat as associated with ephemeral drainages throughout a large portion of the occupied landscape, while the WYNDD model predicted habitat over larger, more contiguous areas containing both ephemeral drainages and intervening upland areas. It is fortuitous to have two species distribution models for a species in the same geographic area; however, the degree to which upland areas within the predicted range should be considered habitat remains unclear. Resource managers who wish to know the amount of pygmy rabbit habitat present in an area or wanting to direct gas drilling efforts away from pygmy rabbit habitat may want to consider both models in an ensemble manner (e.g., Stohlgren et al. 2010; Jones-Farrand et al. 2011), where more confidence is placed in mapped areas (i.e., pixels) for which both models agree than for areas where there is model disagreement.

Neither model was designed to account for anthropogenic disturbance, but disturbance is manifesting throughout southwestern Wyoming in the form of exurban expansion, off-road recreational vehicle use, and energy development (Finn and Knick 2011; Rowland and Leu 2011). The effect of these stressors on sensitive wildlife should be quantified because most species do not have sufficient legal protection to ensure that development will be directed away from areas where they are predicted to occur. Although our study was not specifically designed to relate gas-field road network density with pygmy rabbit occupancy rates, the observed effect of road density on model prediction success demonstrates a negative relationship. The fact that both models had significantly higher commission error rates in areas with higher road densities suggests that habitat quality, and thus the likelihood of pygmy rabbit presence, declined as road density increased. This offers support for another study that showed pygmy rabbits to be sensitive to disturbance at a more localized spatial scale (e.g., Lawes et al. 2012). The observed relationship between road density and model prediction success also demonstrates the importance of incorporating relevant anthropogenic stressors when developing habitat models intended to accurately reflect current landscape conditions.

In conclusion, two species distribution models exist for pygmy rabbits in southwestern Wyoming. The models did not differ in overall classification success rates, both models were better at predicting habitat than nonhabitat, and the predictive success of both models were negatively influenced by road networks associated with energy development. Because of the level of anthropogenic disturbance occurring in southwestern Wyoming, models that incorporate anthropogenic disturbance are needed in order to accurately account for disturbance effects. Until one is produced, when pygmy rabbit distributions are a management concern, users may want to consult both current models. In areas with model agreement (58.2% of the jointly mapped area), confidence in their joint prediction will be higher than when using either model alone (e.g., Brieman 1996). Where models disagree, it may be useful to consider that the TNC model had a higher measure of sensitivity and so may be the model of choice when mapping habitat is the goal. However, it is equally useful in many management scenarios to know where nonhabitat is predicted, and the WYNDD model had a higher specificity rate and therefore may be the model of choice in these situations.

Please note: The Journal of Fish and Wildlife Management is not responsible for the content or functionality of any supplemental material. Queries should be directed to the corresponding author for the article.

Data S1. Data set containing the model evaluation survey results, predicted habitat class, and percent area predicted to be habitat within each of three buffer radii for The Nature Conservancy of Wyoming 2009 pygmy rabbit Brachylagus idahoensis species-distribution model.

Found at DOI: http://dx.doi.org/10.3996/022014-JFWM-016.S1 (14 KB XLSX).

Data S2. Data set containing the model evaluation survey results, predicted habitat class, and percent area predicted to be habitat within each of three buffer radii for the Wyoming Natural Diversity Database 2008 Wyoming pygmy rabbit Brachylagus idahoensis species-distribution model.

Found at DOI: http://dx.doi.org/10.3996/022014-JFWM-016.S2 (14 KB XLSX).

Data S3. Data set containing site and plot identifiers, evaluation-survey-based pygmy rabbit occupancy results, and road density (km/km²) levels at plots predicted to be occupied by pygmy rabbits by The Nature Conservancy of Wyoming 2009 pygmy rabbit Brachylagus idahoensis species-distribution model for southwestern Wyoming.

Found at DOI: http://dx.doi.org/10.3996/022014-JFWM-016.S3 (12 KB XLSX).

Data S4. Data set containing site and plot identifiers, evaluation-survey-based pygmy rabbit occupancy results, and road density (km/km²) levels at plots predicted to be occupied by pygmy rabbits by the Wyoming Natural Diversity Database 2008 Wyoming pygmy rabbit Brachylagus idahoensis species-distribution model for southwestern Wyoming.

Found at DOI: http://dx.doi.org/10.3996/022014-JFWM-016.S4 (13 KB XLSX).

Reference S1. Keinath DA, McGee M. 2004. Species assessment for pygmy rabbit Brachylagus idahoensis in Wyoming. Wyoming Natural Diversity Database Report to U.S. Department of the Interior Bureau of Land Management, Cheyenne, Wyoming.

Found at DOI: http://dx.doi.org/10.3996/022014-JFWM-016.S5; also available at http://www.uwyo.edu/wyndd/_files/docs/reports/speciesassessments-/pygmyrabbitmar2004.pdf (3538 KB PDF).

Reference S2. Ulmschneider H, Hays D, Roberts H, Rachlow J, Forbes T, Himes J, Sequin E, Haworth M, Katzner T, Kozlowski A, Rausher R, Lauridson P. 2004. Surveying for pygmy rabbits Brachylagus idahoensis. 4th draft. Boise, Idaho: U.S. Department of the Interior Bureau of Land Management.

Found at DOI: http://dx.doi.org/10.3996/022014-JFWM-016.S6; also available at http://sagemap.wr.usgs.gov/docs/DraftPygmy%20RabbitProtocol6_10_04.doc (149 KB PDF).

We thank K. Hughes, J. Jewell, S. Meloy, D. Probasco, and A. Weiwel for collecting evaluation survey data; D. Woolwine for assisting with survey training; and M. O'Donnell for providing GIS support. B. Collier, C. Jarnevich, and two anonymous reviewers provided valuable editorial comments that greatly improved the manuscript. We also thank Subject Editor S. Jones for facilitating the manuscript review. Funding for this study was provided through the U.S. Geological Survey Wyoming Landscape Conservation Initiative (WLCI).

Any use of trade, product or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government.