Species distribution models (SDMs) use information about where a species occurs to determine the relationship between species occurrence and environmental variables and thus create maps that predict the species' distribution. SDMs are useful for studying poorly understood species, species of conservation concern, and invasive species. These models provide information about what environmental variables influence the species' distribution, identify both core and marginal habitat important to the conservation (or eradication) of the species, and pinpoint areas where surveys or monitoring should occur (Elith et al. 2006; Urbina-Cardona and Flores-Villela 2010; Gormely et al. 2011; Peterman et al. 2013; Murray et al. 2015).

The Oscura Mountains Colorado chipmunk (Neotamias quadrivittatus oscuraensis) is a rare subspecies of the Colorado chipmunk endemic to the Oscura Mountains in central New Mexico (Sullivan 1996). The Colorado chipmunk is a montane species that occurs at high elevations (1,380 to 3,360 m) in coniferous forests of Colorado, New Mexico, and Arizona (Best et al. 1994). Studies have suggested that rocky areas, conifer trees, shrubs, and herbaceous plants, are important to this species (Bergstrom 1986; Best et al. 1994), including in the Oscura Mountains (Sullivan 1996; Rivieccio et al. 2003). The Oscura Mountains Colorado chipmunk is completely isolated from other populations of the Colorado chipmunk because its montane habitat on the Oscura Mountains is surrounded by lower-elevation desert and grassland (Ditto and Frey 2007). Colonization of this habitat island likely occurred during the Wisconsin glaciation (33,000–12,000 years ago), when cooler conditions resulted in downward elevational displacement of vegetation communities and connectivity of coniferous habitat between the Oscura Mountains and the Southern Rocky Mountains (Patterson 1980).

The Oscura Mountains Colorado chipmunk is listed as threatened by the state of New Mexico, primarily because it is an endemic taxon with a small, isolated habitat and the potential for continuing habitat loss (New Mexico Department of Game and Fish 2016). The most basic information needed to manage and conserve a species is its distribution and habitat selection. However, knowledge about the ecology of the Oscura Mountains Colorado chipmunk is limited, and most studies and monitoring efforts have not been published or adequately documented. Studies of the Colorado chipmunk suggested that coniferous woodlands and rocky cliffs may be two of the most important factors influencing the chipmunk's distribution on a landscape scale (Best et al. 1994; Sullivan 1996; Sullivan and Wilson 2000), but these hypotheses have not been rigorously tested.

We created an occupancy model for the Oscura Mountains Colorado chipmunk using camera-trapping data to address the foregoing hypotheses and assess what ecological factors influence its distribution (Perkins-Taylor and Frey 2018). Occupancy modeling uses repeated surveys at sites to estimate the likelihood that a species is present at a site while accounting for imperfect detection during surveys (MacKenzie et al. 2006). In general, these models can examine differences in occupancy probability associated with many kinds of environmental variables, including both microhabitat data collected at each site and landscape-scale data that can be spatially represented as layers in a geographic information system (GIS). However, when only landscape-scale covariates are used, the resulting occupancy probabilities can be mapped and the model can be considered a SDM. These models require carefully designed sampling methods to assess the occupancy status of sites, which can require significant time and resources. Like many other hierarchical models, occupancy models are more “data hungry” than more simple regression models, i.e., the complexity of estimating occupancy and detection probabilities simultaneously requires a larger sample size to obtain precise estimates of coefficients and standard deviations, particularly when many variables are included. Our occupancy analysis identified three important landscape-scale features that influence the chipmunk's distribution: piñon (Pinus edulis) woodland, proximity to an escarpment, and elevation (Perkins-Taylor and Frey 2018). However, we were concerned that our sample size and the data hungry nature of occupancy modeling restricted the number of variables we could test and the model's statistical power to identify significant variables. In particular, we were concerned that: (1) our “proximity to escarpment” variable was very crude because of our inability to test multiple topographic variables related to rocky cliffsides, and (2) the occupancy probabilities estimated in piñon woodland areas far from an escarpment were unexpectedly low for a species that relies on coniferous woodlands. These concerns are well addressed using a different SDM method, MaxEnt.

MaxEnt is a program that uses the principles of maximum entropy to create SDMs from presence-background data (Phillips et al. 2006). MaxEnt compares occurrence points to background points using a machine learning algorithm to select environmental variables that best discriminate between them. Because true absences are not needed, MaxEnt can use nonsystematically collected occurrence records from a variety of sources (Elith et al. 2006). The output is a relative occurrence rate (ROR) across the landscape, i.e., the relative probability that a cell is contained in a collection of presence samples (Fithian and Hastie 2013), which can be interpreted as an index of habitat suitability (Merow et al. 2013). MaxEnt's machine learning algorithm is capable of testing a large number of covariates and has been shown to perform well with small sample sizes (D'Elia et al. 2015). Machine learning methods also perform better than traditional methods when analyzing complex, nonlinear relationships (Shoemaker et al. 2018). The strengths of MaxEnt therefore are well suited to investigating the validity and potential weaknesses of our occupancy model.

Our goal was to create a SDM for the Oscura Mountains Colorado chipmunk using MaxEnt and compare the resulting relative likelihood of occurrence and the ecological covariates identified with those of our previous occupancy model (Perkins-Taylor and Frey 2018). Our aim was to either strengthen our original assessment of the relevant ecological factors responsible for the chipmunk's distribution or to identify additional factors that were not captured by our occupancy model, thereby providing an additional tool for management of this taxon of conservation concern. To generate the MaxEnt models, we designed a novel model selection process for tuning model parameters while also selecting environmental variables based on model performance. Previous studies have cautioned that MaxEnt can obtain biased results when spatial bias and imperfect detection are ignored (Rota et al. 2011; Lahoz-Monfort et al. 2014). We therefore created two MaxEnt models with and without data from nonsystematically collected occurrence records that are subject to greater spatial bias and bias from imperfect detection to examine how detection probability may have influenced our results.

Materials and Methods

Study area.—The Oscura Mountains are encompassed within White Sands Missile Range, in Lincoln and Socorro counties, New Mexico ( Supplementary Data SD1). They are north-south trending and approximately 40 km long and 6–10 km wide. Elevations range from approximately 1,500 to 2,627 m. The Oscura Mountains are part of the fault block associated with the Rio Grande Rift, which results in a steep west face with limestone outcrops, and a more gently sloping east face dominated by conifer woodlands ( Supplementary Data SD1). Our study area consisted of the Oscura Mountains and adjacent areas bounded by northing 3747027.456–3702297.456 and easting 360415.165–404275.165 in UTM zone 13S (datum: WGS 1984).

The vegetation communities present on White Sands Missile Range were classified and mapped using satellite imagery (Muldavin et al. 2000a, 2000b). There are five vegetation community types present within our study area on the Oscura Mountains and immediately surrounding areas (Muldavin et al. 2000b). (1) Piñon woodland is characterized by moderately closed canopies (25–60% cover) dominated by piñon pine (P. edulis), with oneseed juniper (Juniperus monosperma) subdominant or codominant, but never more than 50% of the canopy, and is located at higher elevations (1,540 to 2,627 m), mainly along the gradual eastern slopes. (2) Juniper woodland is dominated by oneseed juniper and occurs at lower elevations (1,494 to 2,498 m) than the piñon woodland, along the eastern flank of the mountains. Juniper woodlands are savannah-like, with open canopies (10–40% cover) and grassy inter-tree spaces; piñon pine usually is entirely absent, although a few individuals occasionally may be present. (3) Montane scrub is dominated by mountain mahogany (Cercocarpus montanus) with wavyleaf oak (Quercus undulata) as a common associated species, and usually occurs on slopes and ridges (1,506 to 2,606 m), particularly the western slope and bottom of the eastern slope. (4) Interior chaparral is characterized by shrub live oak (Quercus turbinella) communities, usually found on low- to mid-elevation slopes (1,524 to 2,321 m) near juniper woodlands or grasslands. (5) Foothill-montane temperate grassland is characterized by blue grama (Bouteloua gracilis), New Mexico needlegrass (Achnatherum perplexum), and western wheatgrass (Pascopyrum smithii) dominated communities and occurs on mountain slopes or in mountain valleys (1,384 to 2,595 m).

Occurrence points.—We surveyed sites for the Oscura Mountains Colorado chipmunk during May–August 2015 and 2016 using baited camera traps in an occupancy model design (Perkins-Taylor and Frey 2018). We placed sites within 500 m of roads due to accessibility, but at least 120 m away from roads to minimize impact of traffic noise (Chen and Koprowski 2015), and at least 120 m away from other camera sites. Using these spacing rules, we selected locations to survey using a random design stratified by vegetation community type (Muldavin et al. 2000b). Cameras were deployed at each site for a minimum of 5 days in 2015 and 7 days in 2016. Identification was noncontroversial, because the only co-occurring sciurids in the Oscura Mountains are the rock squirrel (Otospermophilus variegatus) and Texas antelope squirrel (Ammospermophilus interpres). The Oscura Mountains Colorado chipmunk is distinguishable easily in photographs by its black and white striping along the face and back, as well as its body size and brownish coloration. All field methods conformed with the ASM guidelines (Sikes et al. 2016) and were approved by the New Mexico State University institutional animal care and use committee.

In addition to our camera surveys, we recorded opportunistic sightings of the Oscura Mountains Colorado chipmunk while traveling to sites, both while driving and walking. We also obtained occurrence records collected by White Sands Missile Range from 2004 to 2009 using various methods (D. W. Burkett, Eco Inc., pers. comm.). However, because most of these records lacked physical evidence, we only used those for which we were fairly confident of the species identification and the accuracy of the geographic coordinates based on the survey methodology (livetrapping or visual identification during walking surveys from Burkett, or our own opportunistic sightings). Visual surveys are likely subject to much greater influences of imperfect detection than our baited camera surveys (Perkins-Taylor 2017). We controlled for false presences by removing occurrence records that had the potential for misidentification based on the methodology, but false absences during visual surveys could skew the models if detection probabilities vary across the landscape. Chipmunks are easier to observe when perched on rocks or cliff edges, or when running across roads, meaning that these opportunistic records may exhibit more spatial bias than the camera data.

To determine whether the occurrence locations from our camera surveys were sufficient to assess the relationship between the environmental variables and the distribution of the Oscura Mountains Colorado chipmunk or whether including the opportunistic occurrence records produced different results, we created models with (n = 49) and without (n = 36) the opportunistic records (Fig. 1). For both models, presence points were rarefied at 200 m to reduce clustering and spatial bias. This method of spatial filtering is one of the best methods for correcting sampling bias (Kramer-Schadt et al. 2013; “systematic sampling” in Fourcade et al. 2014). After rarefying the presence points, the model with both camera trapping and opportunistic records (n = 49) contained 29 records from camera trapping, six from our opportunistic sightings while deploying cameras (Perkins-Taylor and Frey 2018), four from livetrapping conducted by White Sands Missile Range, and 10 from visual observations during walking surveys by White Sands Missile Range.

Background extent.—The background extent limits where background points are drawn from, and thus the range of values of the environmental variables that MaxEnt will compare to the occurrence points. The choice of background extent can influence the shape of response curves and even which variables are selected by the machine learning algorithm (Phillips et al. 2009; Anderson and Raza 2010; Barve et al. 2011); thus, the background extent should have an ecological justification based on the questions being asked (Merow et al. 2013). To use the results of species distribution modeling to understand how environmental factors are related to habitat selection, background points should be drawn from areas accessible to the species (Barve et al. 2011), thereby comparing used sites to accessible sites. We therefore used a background extent based on a 2-km buffer around all presence points, which we considered a reasonable maximum dispersal distance based on information for this and similar species of chipmunk (Meredith 1974; Loew 1999; Sullivan and Wilson 2000). Like our prior occupancy model, our SDM is an order 1 resource selection function, which is selection of the location of populations throughout the range of the species (sensu Meyer and Thuiller 2006). The camera points were designed to encompass the entire range of the Oscura Mountains Colorado chipmunk (Perkins-Taylor 2017), and our background extent also covered the whole known range; the model therefore will compare occupied areas within the geographic range to background points (available but presumed to often be unoccupied) within the geographic range.

Fig. 1.

(A) Map of occurrence records used for MaxEnt Model 1 for Neotamias quadrivittatus oscuraensis in the Oscura Mountains, New Mexico. All points were taken from camera surveys conducted in an occupancy modeling framework from 2015 to 2016. (B) Map of occurrence records used for MaxEnt Model 2. Points in gray were taken from camera surveys conducted in an occupancy modeling framework from 2015 to 2016. Points in black were taken from opportunistic sightings while traveling to survey sites, or from previous, nonsystematic surveys by White Sands Missile Range from 2004 to 2009. All points were rarefied at a 200-m scale. Black lines represent roads.

Environmental variables.—We tested 35 environmental variables related to vegetation community type, topography, climate, and landscape productivity; we hypothesized these variables to be associated with habitat for the Oscura Mountains Colorado chipmunk (Table 1) based on previous studies (Sullivan 1996; Rivieccio et al. 2003; Perkins-Taylor and Frey 2018). All rasters were scaled to 30 m using ArcGIS 10.3 (ESRI, Redlands, California). The vegetation community type and elevation rasters were provided by White Sands Missile Range, and we calculated other topographic rasters from these data using tools in ArcGIS (see  Supplementary Data SD2). We calculated the number of steep neighbor cells, topographic position index (TPI—Weiss 2001), and vector ruggedness measure (VRM—Sappington et al. 2007) using a moving window approach with 1, 3, and 10-cell windows. We calculated hillshade values for 9 a.m., 12 p.m., and 3 p.m., on the 2016 solstices and equinoxes, and took the average of these values to account for daily and seasonal differences in local illumination. In order to calculate hillshade values, we obtained data describing the sun's altitude and azimuth for each of these dates and times from the United States Navy Astronomical Applications Department ( http://aa.usno.navy.mil/data). We obtained 19 standard bioclimatic variables from the WorldClim website ( http://www.worldclim.org/tiles.php). We obtained annual normalized difference vegetation index (NDVI) data, a measure of greenness often interpreted as representing productivity, from the United States Geological Survey ( https://weld.cr.usgs.gov/) for 2003–2012 and averaged these values.

Model complexity and variable selection.—MaxEnt (Version 3.3.3) has several modifiable parameters, and while many studies use the default parameter settings (Morales et al. 2017), there is a large amount of evidence that using the default settings may not generate the best model (Shcheglovitova and Anderson 2013; Syfert et al. 2013; Radosavljevic and Anderson 2014; Morales et al. 2017). There are two parameters that control model complexity: feature classes and the regularization multiplier. Features are transformations of environmental variables used to build the model, including linear, quadratic, threshold, hinge, product, and categorical features (Merow et al. 2013). The default settings allow certain feature classes to be used based on the sample size of occurrence points (Phillips and Dudík 2008): linear, quadratic, and categorical features are always tested, hinge features are tested when n ≥ 15, and threshold and product features are tested when n ≥ 80. The regularization multiplier, or beta parameter, influences the degree of l₁-regularization, which penalizes a model's maximum likelihood relative to the number of features in the model. The default value for beta is 1 (Phillips and Dudík 2008). The beta and feature class parameters influence the complexity of the model by controlling the number and complexity of features in the final output. However, use of default parameter settings may result in overly complex or simple models (Warren and Seifert 2011), and recent findings suggest that different combinations of feature classes and regularization multipliers should be tested to choose the best parameter settings for the species and data in question (Merow et al. 2013; Syfert et al. 2013; Radosavljevic and Anderson 2014). We therefore tested all possible combinations of feature classes and regularization multipliers from 0 to 15 in increments of 0.5 (Jueterbock et al. 2016).

Table 1.

Environmental variables used in MaxEnt modeling that we identified as potentially related to the distribution of Neotamias quadrivittatus oscuraensis in the Oscura Mountains, New Mexico.

In addition to tuning model parameters, another often overlooked recommendation when building MaxEnt models is to remove correlated variables using objective criteria, such as identifying the variable with the highest percent contribution to the model and then removing any other variables that are correlated with it (Warren et al. 2014). MaxEnt is more stable in the face of correlated variables than other methods, such as stepwise regression (Elith et al. 2011), but MaxEnt models with strongly correlated variables still can behave erratically (Warren et al. 2014). Many studies simply calculate correlations between environmental layers and rely on expert opinion to decide which of the correlated variables to remove. However, this circumvents one of the major advantages of using MaxEnt to create SDMs, which is that the machine learning process can evaluate many variables and pick the ones that best fit the data. Relying on expert opinion also can be difficult when studying a poorly understood species, such as the Oscura Mountains Colorado chipmunk. Therefore, our model selection process was designed to remove correlated variables in a systematic fashion based on their contribution to the model.

Fig. 2.

Diagram depicting our four-step model selection process for the MaxEnt models. This process adjusted model complexity by tuning the regularization parameter (β) and feature classes. Variable selection was done iteratively based on variable contribution to the model instead of a priori. AIC_c = Akaike's Information Criterion with small sample size correction.

Based on the above recommendations from the literature, we designed a model selection process based on Warren et al. (2014) and Muscarella et al. (2014) that allowed us to tune both the regularization multiplier and the allowable features, while also selecting environmental variables based on model performance. To our knowledge, this is the first method to tune both the regularization multiplier and feature classes as well as to undertake variable selection in a systematic fashion based on the MaxEnt output. Muscarella et al. (2014) described a model selection process for tuning the regularization multiplier and feature classes, and Warren et al. (2014) described a process for tuning the regularization multiplier and performing variable selection, but all three have never been done simultaneously, presumably because it is computationally intensive. Our model selection process had four steps (Fig. 2). We used Akaike's Information Criterion with small sample size correction (AIC_c—Burnham and Anderson 2002) as our model selection criteria because this method balances model complexity with goodness-of-fit (Warren and Seifert 2011; Galante et al. 2018). We conducted all analyses in R statistical software (R Version 3.1.3,  www.r-project.org, accessed 4 April 2015) using the “ENMEval” package (Muscarella et al. 2014).

In the first step, we removed extremely correlated (r > 0.99) variables from the data set. Variables that are correlated with each other to such an extreme degree have the potential to greatly influence the models and could sway the choice of model parameters, so we elected to remove these extremely correlated variables first. To do this, we created models using all possible feature class combinations (n = 31) and regularization multipliers from 0 to 15 in increments of 0.5 (n = 961 models). We chose the most informative model, while controlling for overfitting, using AIC_c (Muscarella et al. 2014; Warren et al. 2014; Jueterbock et al. 2016; Morales et al. 2017). From this model, we retained the variable with the highest contribution and removed any variables that were extremely correlated with it (r > 0.99). Of the remaining variables, we then retained the variable with the second highest contribution and removed variables that were extremely correlated with it, and so on until all variables either were retained or removed. This resulted in a reduced set of environmental variables that were not extremely correlated with each other.

In the second step, we tuned the feature classes and regularization multiplier. Again, we created models using all possible feature class combinations and regularization multipliers with the reduced set of variables from step 1. We chose the most informative model using AIC_c (i.e., ΔAIC_c = 0). We then used the feature class combination and regularization multiplier from this model as the model parameters for variable selection in step 3.

In the third step, we selected the best set of uncorrelated variables from the reduced set of environmental variables obtained in step 1. We did this by creating a model using the reduced set of variables from step 1 and the model parameters from step 2. We then retained or removed variables using the following guidelines. Any variables with less than 5% contribution were removed. The variable with the highest percent contribution to the model was retained, and any variables that were highly correlated with it (r > 0.7) were removed, resulting in a reduced set of variables. We then created a new model using this reduced set of variables. We again removed any variables with less than 5% contribution, retained the variable with the highest contribution (apart from the previously retained variable), and removed any remaining variables that were highly correlated with it (r > 0.7). We continued this process until all variables either were retained or removed, resulting in a final reduced set of variables that all contributed to the model but were not highly correlated. We did not measure correlations between categorical and continuous variables, as the Pearson correlation (r) is not appropriate when there are more than two categories. Our one categorical variable, vegetation type, therefore was only subject to removal if it contributed less than 5% to the model.

In the fourth step, we used the final reduced set of variables from step 3 to retune the model parameters. Using only the variables retained during step 3, we created models using all possible feature class combinations and regularization multipliers, and selected the model with the lowest AIC_c value as the final model.

We created separate models for each set of occurrence points. We calculated the relative likelihood of occurrence on the logistic scale and calculated the area under the receiver operating curve (AUC) to assess model fit, although this statistic can be falsely inflated by spatial autocorrelation (Veloz 2009). We separated occurrence points into training and test data using 4-fold geographically masked cross-validation following the checkerboard2 method of Muscarella et al. (2014). This method of data partitioning accounts for potential spatial bias when partitioning the data into training and test data (Veloz 2009; Radosavljevic and Anderson 2014), which is necessary to calculate statistics for model evaluation, such as AUC and the percent contribution of each variable.

Model comparison.—We projected the models created using both sets of occurrence records across the entire study area. We compared these two maps using three similarity statistics. Schoener's D statistic (Schoener 1968) measures niche overlap. It has a history of use in ecological literature, allowing for direct comparison to traditional measures of niche similarity; however, it may suggest a biological interpretation of cell values that is unwarranted (Warren et al. 2008). Warren et al. (2008) proposed a niche similarity statistic that carries no biological assumptions, but simply treats the cell values of each map as a probability distribution: the I statistic. The relative rank statistic (Warren and Seifert 2011) is an estimate of the probability that the relative ranking of any two patches of habitat is the same for the two models irrespective of the quantitative difference in relative likelihood (i.e., habitat suitability) estimates, and thus assesses the similarity of the overall patterns of each map.

We also used the relative rank statistic to compare the MaxEnt outputs to a map of occupancy probability for the Oscura Mountains Colorado chipmunk from our previous occupancy model based on the camera-trapping data (Perkins-Taylor and Frey 2018). We compared the variables selected by each model to examine whether the three models identified similar ecological factors as influencing occurrence. It is important to note that unlike occupancy models, presence-background methods cannot estimate a true probability of occurrence. At best, they estimate a relative likelihood of occurrence that is proportional to the actual probability of occurrence, given that detection probability is constant and there is no sampling bias (Guillera-Arroita et al. 2015). The detection probability for our surveys was not constant, but because it was positively correlated with occupancy probability (Perkins-Taylor and Frey 2018), our MaxEnt models will correctly rank the suitability of sites, although the shape of response curves may not represent the true relationships between environmental variables and probability of occurrence (Guillera-Arroita et al. 2015). Because both the occupancy model and the MaxEnt models will correctly rank likelihood of occurrence, it is appropriate to compare the two methods using the relative rank statistic. It is inappropriate to compare them using Schoener's D or the I statistic because the two modeling methods are computing different values (i.e., true probability of occurrence for the occupancy model and relative likelihood of occurrence for the MaxEnt models).

Results

The best MaxEnt model built with only the camera occurrence points (Model 1) had a regularization multiplier of 3.5, allowed only for linear and categorical features, and had AUC = 0.738. The best MaxEnt model built with both the camera points and opportunistic records (Model 2) had a regularization multiplier of 2.5, allowed linear, quadratic, and categorical features, and had AUC = 0.737. Both models included elevation (linear feature in Model 1, quadratic feature in Model 2) and vegetation community type (categorical feature) as the variables with the highest contribution. The third variable in Model 1 was the linear feature topographic position index with a 10-cell window (TPI-10), which had a percent contribution and permutation importance similar to elevation (Table 2). The third variable in Model 2 was the linear feature vector ruggedness measure with a 3-cell window (VRM-3), which had the lowest percent contribution and permutation importance of the three variables in Model 2. Model coefficients and response curves for Model 1 showed that relative likelihood of occurrence was positively related to elevation, the presence of piñon woodland vegetation type, and TPI-10 (Table 2; Fig. 3A). Model coefficients and response curves for Model 2 showed that relative likelihood of occurrence was positively related to elevation, presence of piñon woodland vegetation type, and VRM-3, and negatively related to montane scrub vegetation type (Table 2; Fig. 3B).

For both MaxEnt models, the highest relative likelihood values were in high-elevation areas near the west-facing escarpment (Fig. 4). The eastern slopes, which are generally covered in piñon woodland, had moderate relative likelihood values, and the low-elevation basin areas had low relative likelihood values. The relative likelihood values predicted across the study area by both MaxEnt models were similar based on Schoener's D and the I statistic (Table 3). The overall geographic patterns of likelihood of occurrence predicted by both MaxEnt models were also similar based on the high relative rank statistic (Table 3).

Table 2.

Percent contribution, permutation importance, and coefficients of the variables in the best MaxEnt models for Neotamias quadrivittatus oscuraensis in the Oscura Mountains, New Mexico, built using only occurrence points from cameras (Model 1) or occurrence points from cameras and opportunistic observations (Model 2). For vegetation type, coefficients are only reported for the categories used by the model. Vegetation types are: 0—other, unsuitable vegetation types such as roads and nonmontane categories, 1—montane scrub, 2—interior chaparral, 3—piñon woodland, 4—foothills-montane temperate grasslands, 5—juniper woodland.

Fig. 3.

Response curves for variables included in (A) MaxEnt Model 1 and (B) MaxEnt Model 2 for Neotamias quadrivittatus oscuraensis in the Oscura Mountains, New Mexico. Relative likelihood of occurrence is predicted for the full range of values of each variable at the background points while holding all other variables at their median values. Vegetation types are: 0—other, unsuitable vegetation types such as roads and nonmontane categories, 1—montane scrub, 2—interior chaparral, 3—piñon woodland, 4—foothills-montane temperate grasslands, 5—juniper woodland. TPI-10 = topographic position index calculated using a 10-cell moving window. VRM-3 = vector ruggedness measure calculated using a 3-cell moving window.

Our occupancy model (Perkins-Taylor and Frey 2018) included occupancy covariates for elevation, piñon woodland vegetation community type, and proximity to an escarpment (Table 4). The model had sharp distinctions between different portions of the Oscura Mountains, with areas near an escarpment having almost 100% occupancy probability, areas of piñon woodland far from an escarpment having ∼30% occupancy probability, and low-elevation basin areas having < 1% occupancy probability (Fig. 4C). In contrast, both MaxEnt models predicted a more gradual decrease in relative likelihood of occurrence associated with decreasing elevation, which resulted in much more nuanced maps of the species' distribution (Fig. 4). The overall geographic pattern of occupancy probability predicted by the occupancy model was similar to the MaxEnt models based on the high relative rank statistic (Table 3).

Fig. 4.

Relative likelihood of occurrence on the logistic scale for Neotamias quadrivittatus oscuraensis in the Oscura Mountains, New Mexico predicted by (A) MaxEnt Model 1 and (B) MaxEnt Model 2, and (C) occupancy probability predicted by the final occupancy model (Perkins-Taylor and Frey 2018). Areas of high likelihood are represented by warmer colors, and areas of low likelihood are represented by cooler colors. The background extent for each MaxEnt model is outlined in black.

Table 3.

Statistics comparing the MaxEnt and occupancy models for Neotamias quadrivittatus oscuraensis in the Oscura Mountains, New Mexico, projected across the study area. Model 1 is the MaxEnt model created using only occurrence points from cameras; Model 2 is the MaxEnt model created using the occurrence points from cameras and opportunistic observations.

Discussion

Similarities among models.—The general pattern of distribution predicted by the MaxEnt models and the occupancy model conformed to our expectations for a montane chipmunk species. All three models included covariates describing the influence of elevation, piñon woodland vegetation community type, and a topographic variable associated with steep escarpments. All three models also predicted the same spatial patterns of increasing occupancy probability or relative likelihood of occurrence with increasing elevation, with the highest values along the top ridge of the mountain range. These patterns are consistent with the montane distribution of Colorado chipmunks and their close association with vegetation communities dominated by conifers in the Pinaceae (Lechleitner 1969; Best et al. 1994; Sullivan 1996). In addition, other studies have found that the Colorado chipmunk is most common in areas with rocky cliffs and ledges (Bergstrom 1986; Best et al. 1994; Sullivan 1996).

The similarities between our MaxEnt and occupancy models lead us to the conclusion that we can be confident that our understanding of the general distribution of the Oscura Mountains Colorado chipmunk and the ecological factors that influence the species are not being skewed too greatly by the weaknesses of one particular modeling method. One concern regarding our previous occupancy model (Perkins-Taylor and Frey 2018) was that our limited sample size and the data hungry nature of occupancy modeling may have resulted in an overly simple model that disregarded the influence of other ecological factors on probability of occupancy. However, the fact that MaxEnt, which uses a machine learning algorithm capable of testing a large number of covariates and has been shown to perform well with small sample sizes (D'Elia et al. 2015), also selected similar environmental variables increased our confidence that we have a sufficient understanding of the major ecological factors influencing the distribution of the Oscura Mountains Colorado chipmunk.

Our novel method of both tuning model complexity and selecting variables based on model importance greatly reduced bias in what variables were selected by ensuring that our MaxEnt models were not overfit and by providing an objective method of removing correlated variables. These two steps control for bias in the model-building process in different ways. We therefore believe that a robust model selection process should include both model-building steps. If model complexity is not tuned according to current standards (Morales et al. 2017), then the environmental variables chosen by MaxEnt's machine learning algorithm may overfit the occurrence points and not accurately describe the ecological factors that more generally influence the species' distribution (Warren and Seifert 2011; Shcheglovitova and Anderson 2013). And had we not created an objective method for removing correlated variables, then we would not know if the similarities between our MaxEnt and occupancy models were due to our own subjective choices during the model selection process. Thus, the model selection process played an important role in verifying the similarities between our models.

Table 4.

Parameter estimates and standard errors for the final occupancy model for Neotamias quadrivittatus oscuraensis in the Oscura Mountains, New Mexico (Perkins-Taylor and Frey 2018). Mature piñon is a categorical variable describing presence or absence of any mature piñon trees (root crown diameter ≥ 20 cm) within 30 m of the camera. Piñon vegetation type is a categorical variable describing if the camera site was classified as piñon woodland based on the vegetation community type layer. Proximity to escarpment is a categorical variable describing if the camera was within 150 m of an escarpment. Values are on the logistic scale.

One critique of this model selection method is that we evaluated more models than the number of occurrence points, which is contrary to the traditional approach when using AIC (Burnham and Anderson 2002). We do not believe that this biased our results. The greatest concern from evaluating a large number of models is that by chance will randomly happen to fit the data well. However, we found that in each model selection step, models that were similar to the best model also had lower AIC_c values and that the combinations of model parameters chosen as the best model remained fairly consistent throughout each step of the model selection process. Our model selection process was based on current methods used in the machine learning context (Muscarella et al. 2014; Warren et al. 2014), and we believe that the alternative of arbitrarily relying on expert opinion presents far greater biases. Still, the issue of evaluating a large number of models in the MaxEnt model selection process may warrant further investigation.

For our MaxEnt models, the similarities with the occupancy model suggest that spatial bias of the occurrence points did not greatly influence the results. Of the three models, we would have expected MaxEnt Model 2 to be most influenced by spatial bias because it used opportunistic occurrence points that were not systematically collected (20 of 49 points). However, the relative rank statistics showed that the pattern of distribution in Model 2 was similar to both Model 1 and to the occupancy model. This indicates that the opportunistic occurrence points in MaxEnt Model 2 did not greatly influence the model's spatial bias (after rarefying the occurrence points and controlling for spatial bias during data partitioning), at least no more so than our systematically collected camera-trapping data. The MaxEnt models also did not appear to be greatly influenced by imperfect detection, based on the overall similarity of the patterns of distribution and the environmental variables selected by the occupancy model and the MaxEnt models.

Differences among models.—Although our MaxEnt models and occupancy model exhibited some broad similarities, there were differences. The maps generated from our MaxEnt models were more detailed and nuanced than those based on the occupancy model, and thus may serve as a more useful tool for scientists making management decisions for this taxon, such as identifying core areas for protection, protecting potential movement corridors, or planning prescribed fires. The MaxEnt models showed gradual, continuous changes in relative likelihood of occurrence across the landscape, whereas the occupancy model showed sharp distinctions between areas of high, moderate, and low occupancy probability. Based on the occupancy model alone, we would have concluded that areas near an escarpment had the highest occupancy probability, and that all other areas classified as piñon woodland vegetation type had a moderately low occupancy probability of ca. 30%. This lower occupancy probability in piñon woodlands was one of the results that made us skeptical of the occupancy model, as the Colorado chipmunk is predominantly associated with coniferous woodlands (Best et al. 1994) and piñon pinecones are hypothesized to be an important food source for the Oscura Mountains Colorado chipmunk (Sullivan and Wilson 2000). The MaxEnt models suggest that there may not be such a stark difference between areas near escarpments and areas farther away from escarpments, particularly areas in piñon woodland. Places where the models differ in predicting occurrence of populations could be evaluated in future surveys.

The differences in nuance of the distribution maps are related to the differences in the variables selected by each model and their coefficients. Both MaxEnt models had a positive relationship between relative likelihood of occurrence and elevation, which was the main factor in both of these models contributing to the pattern of gradual change in likelihood values across the landscape. The occupancy model also had a positive coefficient on the elevation covariate, but this effect was almost entirely cancelled by the negative coefficient on the interaction between elevation and piñon woodland vegetation type. The occupancy model predicted very similar occupancy probabilities in almost all areas of piñon woodland regardless of elevation, which contributed to the sharp distinctions in occupancy between piñon woodland and other vegetation types.

Another important factor contributing to the differences among models was the different topographic variables used in each modeling method. Because of the strengths of MaxEnt's machine learning algorithm, we could test a variety of complex topographic variables related to slope, ruggedness, and aspect, which are characteristics that have been used to describe the Colorado chipmunk's habitat (Bergstrom 1986; Best et al. 1994; Sullivan 1996). For each MaxEnt model, a continuous topographic variable was selected by the variable selection process (TPI-10 for Model 1 and VRM-3 for Model 2), which contributed to the gradual changes in relative likelihood values. When creating the occupancy model, we did not have—because of our limited sample size—the same statistical power to test a large variety of variables to account for the influence of escarpments and rocky outcrops. We therefore chose to use a simple categorical variable that we believed would adequately capture the influence of escarpments. This categorical variable ended up exerting the strongest influence on occupancy probability out of any of the covariates, resulting in the difference between the nearly 100% occupancy probability of areas near an escarpment and the moderate to low occupancy probability of areas farther from an escarpment. Because the differences among the distribution maps were partly due to MaxEnt's ability (and occupancy modeling's inability) to test a large number of variables, we believe that the more gradual patterns of distribution from the MaxEnt models should not be dismissed because they are based on relative likelihoods rather than probabilities of occurrence.

Conservation implications.—For rare and poorly studied species, understanding distribution, and the factors that shape those distribution patterns, can offer analytical challenges and be resource-intensive to evaluate. And yet, it is this very information that is fundamental to enact sound conservation and management plans. For the Oscura Mountains Colorado chipmunk, our SDMs based on maximum entropy modeling unequivocally supported the influence of high elevation, piñon woodland, and rugged terrain on its distribution. However, the MaxEnt models also suggested that the distribution of the chipmunk has less area of high-quality habitat, and is more fragmented, than indicated by the occupancy model. Both of these factors can increase extinction risk.

Our study was able to evaluate the contemporary distribution of the Oscura Mountains chipmunk, but this may not be identical to its historical distribution. The “niche reduction hypothesis” (Scheele et al. 2017) recognizes that anthropogenic disturbance does not occur randomly in geographical space, so the patterns of distribution observed for a species in decline may be influenced by how anthropogenic effects are related to the environmental variables (McDonald et al. 2018). There are not enough historical data available for the Oscura Mountains Colorado chipmunk for us to compare with the results of this study, so we can only speculate how the contemporary patterns of distribution may be due, in part, to anthropogenic effects. It is plausible that declining piñon pine health (Johnson and Smith 2006; Johnson et al. 2014), increased temperatures due to climate change, and possible exposure to plague from other rodent species, could contribute to the contemporary patterns of distribution we observed. However, piñon woodlands and rugged, rocky terrain are the characteristics that we expected to be important based on past observations and the ecology of the species (Bergstrom 1986; Best et al. 1994), so the contemporary distribution is likely to be similar to the historical distribution.

While the occupancy model suggested that escarpments were the most important drivers of the chipmunk's distribution, the MaxEnt models indicated that piñon woodland was a more important driver of its distribution. Piñon woodlands in the southwestern United States are declining due to multiple stressors, including wildfires, drought, climate change, and bark beetle (Ips confuses) outbreaks (Williams et al. 2010, 2013; Meddens et al. 2015). Significant piñon disease and mortality is already occurring in the Oscura Mountains (Johnson and Smith 2006; Johnson et al. 2014). The threats to piñon woodlands, along with their importance to the chipmunk, suggest the need to more carefully manage the piñon woodlands to support the needs of the chipmunk. Additional research is needed to better understand how these chipmunks use the escarpment and piñon woodlands, particularly at the finer macrohabitat and microhabitat scales. Lastly, while occupancy models are useful for developing monitoring strategies, the resulting SDMs may be overly simple due to the data hungry nature of the method. MaxEnt provides an alternative method for creating SDMs that has particular analytical strengths—such as ability to test large number of variables and use small sample sizes—that make it ideally suited for rare and poorly known species.va