EMERGING TECHNOLOGIES Vegetation mapping to support greater sage-grouse habitat monitoring and management: multi- or univariate approach?

Conservation planning for wildlife species requires mapping and assessment of habitat suitability across broad areas, often relying on a diverse suite, or stack, of geospatial data presenting multidimensional controls on a species. Stacks of univariate, independently developed vegetation layers may not represent relationships between each variable that can be characterized by multivariate modeling techniques, leading to inaccurate inferences on the distribution of suitable habitat. In this paper, we examine the role of variable combining in mapping multiple dimensions of greater sage-grouse (Centrocercus urophasianus, GRSG) habitat as a basis for GRSG conservation in the great basin ecoregion within southeastern Oregon. We compare two modeling approaches: a univariate random forest regression model (RF regression) and a multivariate random forest nearest neighbor (RFNN) imputation model , across an array of variables. These include five GRSG habitat descriptor variables: percent cover of trees, juniper, sagebrush, and GRSG food forbs, and the proportion of grasses that are exotic annuals. We also model species distributions of 51 common species in the sage steppe and combine these predictions to estimate alpha diversity. Our results show that RF regression and RFNN can yield univariate predictions with similar performance, but RF regression predictions tend to contain slightly more bias at broader spatial scales. Stacking univariate predictions from RF regression yields covariance errors that manifest as logical errors (juniper cover > tree cover), biases in estimates of GRSG habitat area, and biases in estimates of alpha diversity. Combining variables from the RFNN model does not introduce covariance errors. We conclude that multivariate modeling approaches are better suited to map multidimensional habitat niches at broader spatial scales, and also better suited to provide information for defining multivariable adaptive management triggers at the population level or above.


INTRODUCTION
Sage-steppe vegetation in the western United States dominates many arid and semi-arid landscapes, but land use and land cover change have reduced the capacity of many of these ecosystems to support species of conservation importance, like the greater sage-grouse (GRSG; Davies et al. 2011). Threats to GRSG habitat include expansion of western juniper woodlands (Cook 2015) and exotic annual grass invasions (Arkle et al. 2014). In order to assess key ecosystem attributes across scales, natural resource management and conservation planning for GRSG habitat require spatially explicit vegetation assessment (i.e., mapping). However, major challenges remain in producing maps of vegetation attributes, such as cover by life-form, let alone species. Classified maps of vegetation types, such as LANDFIRE (Rollins 2009) or the USGS GAP vegetation layer (Gergely and McKerrow 2013), can provide useful information on GRSG habitat across its range (Knick et al. 2013, Doherty et al. 2016, but these data are too thematically coarse to fully represent the details needed to thoroughly characterize habitat (there are only six classes of sagebrush vegetation represented in Oregon). Therefore, improved vegetation mapping in arid and semi-arid landscapes is needed to address the management challenges of the sage steppe.
Vegetation mapping efforts designed to inform GRSG habitat designation in arid and semi-arid landscapes must confront significant challenges. Greater sage-grouse habitat is defined by a hyperdimensional niche (sensu Hutchinson 1957). Realistic covariation between many vegetation attributes is a crucial component for illustrating that hyperdimensional niche. Fine-scale studies of GRSG habitat indicate that they respond to a variety of vegetation attributes, including the cover of tall grasses that are used for hiding cover, Wyoming big sagebrush, low sagebrush, the absence of invading juniper trees, and the presence of those forbs that are favored by juvenile birds during the brood rearing season (Peterson 1970, Sveum et al. 1998). Juniper trees serve as perching sites for predators, a landscape feature that has the strongest detrimental effect on grouse in the spring and early summer months while chicks are small. Exotic annual grasses that compete with native species (Rodhouse et al. 2014) directly affect habitat quality. They also lead to larger, more frequent fires (Balch et al. 2013), and wildfires can negatively impact GRSG habitat by reducing sagebrush cover (Davis and Crawford 2015). In addition to the multidimensionality of the GRSG niche, finescale structural variation in woodlands, shrublands, and grasslands caused by small individual plant canopies and phenological variation complicates vegetation mapping (Asner and Heidebrecht 2002). For example, mapping the abundance of different grass species, guilds such as the native bunchgrasses, or even simply grass abundance, has been challenging at fine spatial resolutions (Rodhouse et al. 2014, Ali et al. 2016. Recent research has made substantial progress in mapping vegetation relevant to GRSG habitat. In particular, maps of grouse-relevant variables (e.g., sagebrush cover, shrub height, cover of herbaceous vegetation, and bare ground) have been constructed across Wyoming and more recently the sage-steppe habitat of the western United States (Homer et al. 2012, Xian et al. 2015, Jones et al. 2018. Still, some key vegetation attributes describing GRSG habitat quality have yet to be mapped (e.g., cover of native bunchgrasses, measures of sagebrush cover that are species-specific, and the cover of those native forbs that provide food to young grouse). Therefore, continued research in vegetation mapping in sage-steppe landscapes remains an active field of study.
Despite these recent advances, the use of many related but independently modeled vegetation attributes for wildlife habitat modeling has potential pitfalls. These hazards are exemplified by inference on communities utilizing species list predictions based on multiple univariate, independently fitted species distribution models, otherwise known as stacking (Dubuis et al. 2011, Di Febbraro et al. 2018. Errors in species distribution modeling results are often amplified in maps of species richness as more layers are included (Ferrier and Guisan 2006). When species distribution maps are biased against absences, species presence is over-mapped and diversity estimates are inflated when those maps are stacked (Pottier et al. 2013, Calabrese et al. 2014. Univariate modeling ignores innate covariation in attributes, potentially leading to unrealistic joint predictions, such as co-occurrence of species never observed together (Ovaskainen and Soininen 2011, Pollock et al. 2014. Calabrese et al. (2014) also postulate that missing environmental covariates can lead to characteristic biases in both individual species distribution and macroecological models. Regardless of the source, errors in individual models propagate through stacked species distribution models (SDMs). The accuracy of multivariate predictions can also suffer from missing environmental covariates, but these errors do not compound when variables are recombined. Furthermore, multivariate models can account for some prediction biases for community composition, such as joint species distribution models that model residual covariation in presenceabsence data which can represent both biotic interactions and missing environmental predictors (Pollock et al. 2014, Copenhaver-Parry andBell 2018). Efforts to constrain estimates of composition derived from stacked SDMs have had some success (Guisan and Rahbek 2011, D'Amen et al. 2015, Del Toro et al. 2019). We have not seen an analogous framework that was designed for the process of combining multiple continuous vegetation descriptor variables describing habitat, although D' Amen et al.'s (2015) work combining SDMs with macroecological models describing plant functional traits comes close.
Within the vegetation mapping community, both univariate and multivariate methodologies have attracted substantial interest. Random forest (RF) regression models are well-suited for mapping vegetation attributes because of their ease of use and flexibility with mapping univariate responses (Prasad et al. 2006, Mutanga et al. 2012. They suffer from some constraints with respect to bias. Low and high values are consistently under-represented in model predictions, because standard RF regression model predictions are generated as an average of the predictions of all of the internal regression trees within the ensemble. Fortunately, these biases can be corrected by altering and tuning how the array of values predicted by the members of the RF ensemble are combined into a single summary (Roy and Larocque 2012). For multivariate modeling, k nearest neighbor (kNN) imputation with k = 1 approaches, such as the gradient nearest neighbor (Ohmann and Gregory 2002), have garnered interest because they restrict multivariate predictions in vegetation maps to only those value combinations observed in the data, such that unreasonable pixel-scale combinations of vegetation attributes cannot occur. For example, stacked univariate RF models for species occurrence and abundance in the Western Cascade ecoregion of Oregon, USA, overestimated average plot-level species richness by 29% and 621%, respectively, while a nearest neighbor approach produced estimates no different than plot data on average . For habitat mapping that depends on distinct combinations of multiple vegetation attributes, error accumulation with variable stacking could be problematic, implying a need for joint modeling of species or vegetation attributes. However, little is known about how these types of multivariate models behave as they are transferred from closedcanopy, coniferous forest ecosystems where they perform well (Ohmann et al. 2011) to arid and semi-arid landscapes. Arid and semi-arid landscapes pose different mapping challenges, such as short-term temporal variability (Ali et al. 2016) and high spatial heterogeneity in their structure (Mishra and Crews 2014). Sparsity of observational data is sometimes cited as a motivation for remote sensing of vegetation in arid lands (Lawley et al. 2016), but the scarcity of data also poses a challenge to successful mapping in these regions. Fortunately, this situation is improving, with the inception of the Assessment Inventory and Monitoring (AIM) program within the Bureau of Land Management (BLM; Taylor et al. 2014).
In this study, we leverage both satellite and aerial remote sensing to map vegetation variables often used in defining GRSG habitat, as well as the distributions of 51 common species and an associated measure of alpha diversity across arid and semi-arid landscapes of southeastern Oregon. Specifically, our objectives were to (1) assess the performance of univariate (i.e., RF regression) and multivariate (i.e., nearest neighbor imputation) vegetation mapping approaches and (2) examine the potential consequences of these different modeling approaches for habitat classification. We frame our discussion in the context of generating maps of diversity, species co-occurrence, and habitat classification to reflect the maps' utility to the study of GRSG habitat. The concepts discussed apply generally to applications requiring information on vegetation community condition, which require information on many, non-independent variables simultaneously.

Location
Our study area encompasses 11 million hectares in southeastern Oregon, which includes the majority of the sage-steppe vegetation in the state (Fig. 1). Vegetation types range from salt desert in the rain shadow of the Steens Mountain to the edge of the semi-arid forests and woodlands of the Blue Mountains to the North, and eastern Cascades Mountains to the West. The climate in this region is arid and continental, with hot, dry summers, and cold winters with some snow.

Data
Plot.-Our vegetation dataset for this analysis contains 2709 plots from four different sources. Most of these plots (2531) were drawn from three Bureau of Land Management rangeland monitoring programs. These included 1318 plots from BLM's AIM program, and 981 from their Landscape Monitoring Framework (LMF) program (Toevs et al. 2011), and 232 from the Ecological Site Inventory program (Habich 2001). In addition, we included 313 50-m line-intercept transects collected for the BLM by colleagues within the Institute for Natural Resources, placed to illustrate the effects of recent large fires in the area. AIM and LMF plots are line-intercept transects and species lists, while ESI survey data are ocular estimates within plots that are configured to represent homogeneous patches of vegetation. Although the survey methods differ among the data sources, all are of a ground configuration that is reasonable for use in mapping to 30-m pixels. Plot sizes and transect lengths generally fit within a 3 9 3-or 4 9 4-pixel window, or they are configured to fit within small patches of consistent vegetation on the ground.
Our dataset contained plots collected across a wide range of dates, but we only used recent plots (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) in order to reduce imagery-plot date mismatches, and also because this project aims to represent current vegetation conditions. Most plots (93%) were surveyed between 2012 and 2017. Older plots were included to ensure adequate representation of western juniper (Juniperus occidentalis). This dataset was neither fully random nor regular in its design (although the AIM sub-dataset does have a stratified random design). However, the collection of plots was large enough and well-enough distributed in space to support this analysis.
We performed a manual quality control screen on all of these plots, inspecting the data along with current airphoto imagery to identify plots that appear to mismatch the information in the imagery, either due to an intervening disturbance, noticeable growth in woody vegetation, or due to possible coordinate errors.
Raster.-Our explanatory variables included raster data layers containing information on topography, climate, soil, and remote sensing data (Table 1). Optical remote sensing provided direct measurement of vegetation across all portions of the study area. We used two forms of optical remote sensing: bands and transformations of a normalized mosaic of Landsat 7 (ETM+) imagery showing conditions in 2013 and texture summaries of National Agriculture Imagery Program (NAIP) imagery from 2016. Landsat imagery's 30-m pixels are much larger than the dominant plant functional types in arid and semi-arid landscapes, namely shrubs, grasses, and forbs. For this reason, we also utilized 1-m NAIP imagery that was summarized using texture analysis to describe fine-scale variability (Nielsen and Noone 2014) at the 30-m scale to match Landsat imagery. Thus, optical remote sensing was used to measure both average 30-m spectral properties (Landsat) and sub-pixel-scale variation (NAIP).
Topographic, climate, and soil data layers were also included in order to represent environmentally mediated vegetation variation along with the remote sensing data (sensu Ohmann and Gregory 2002). Topography descriptors included elevation and aspect (Gesch et al. 2002). We derived climate descriptors from 30-yr normals in the PRISM climate data (Daly et al. 2008). The soil raster layers were downsampled from 90 m with bilinear interpolation to 30 m from an array of modeled soil properties' surfaces (Chaney et al. 2019). Due to the large number of raw variables for both the NAIP imagery and also POLARIS soil data layers, variable reduction was performed with principal components analysis to collapse the information content of each set of data layers into fewer raster layers. All raster explanatory variables were extracted at all known plot locations and were available for the entire study region (i.e., full coverage). The data matrix containing both vegetation and explanatory variables formed the basis for subsequent modeling and model assessment. The full list of available variables was reduced with a stepwise reverse selection procedure, removing variables based on a permutation variable importance measure from conditional inference forests (R package party; Strobl et al. 2008) built from each y-variable described above for random forest nearest neighbor (RFNN). Conditional inference forests were chosen for variable selection because they tend to produce more robust estimates of variable importance than ordinary RFs when response variables are highly correlated (Strobl et al. 2007). See Table 1 for a full list and description of raster data variables that were used in all models.

Analysis
Overview. -We modeled vegetation variables using two approaches: (1) univariate RF regression models (one model per response variable) and (2) one multivariate imputation model (one model for all response variables). We assessed the strength of our results for each variable independently and then assessed the covariance of the response variables by recombining them and assessing joint predictions for variable pairs, as well as aggregated predictions of alpha diversity. Finally, we created mapped aggregate predictions of GRSG habitat. Our analysis was conducted within the R Statistical Programming environment (R Core Team 2018). For reference, we provide a diagrammatic representation of our modeling and mapping workflows (Fig. 2).
Response variables.-We built univariate and multivariate model predictions for an array of response variables. The response variables  (3) Imagery-National Airphoto Inventory Program Principal components axes no. 1, 4, 6, 8, 9, 11, 12, 14, 16, 17, 18, 23, 25, 29, 33, 34, and (7) Notes: NAIP, National Agriculture Imagery Program. Sources are (1) Landsat-7 (ETM+) imagery courtesy of the U.S. Geological Survey; (2) Crist and Cicone (1984) include five variables that describe aspects of GRSG habitat, and percent cover for all positively identified species with greater than 2% frequency in our input dataset (51 species). Predicted values from the single-species model predictions were aggregated after modeling to provide predictions of species richness. Four of the habitat-relevant variables describe the combined cover of multiple species: all tree species, all sagebrush species (genus Artemisia), the proportion of all grass cover that is exotic annual species (variable hereafter referred to as "exotic annual grass"), and the percent cover of all likely food forbs for GRSG. The list of likely food forbs for GRSG was estimated from the literature (Klebenow and Gray 1968, Peterson 1970, Drut et al. 1994, Sveum et al. 1998; see Appendix S1). The fifth GRSG habitat variable was percent cover of western juniper.
Random forest regression. -Random forest models are a machine learning ensemble approach, forming a single prediction by aggregating many predictions from an array of internal classification or regression trees (Breiman 2001). We generated univariate predictions of continuous variables with the random forest algorithm in regression mode (hereafter, RF regression), creating one model for each response variable. Random forest models are most often used as a classifier in vegetation mapping (Lawler et al. 2006, Evans and Cushman 2009, Rodriguez-Galiano et al. 2012 Dr agut ß 2016), but RF regression can also be a flexible method for prediction of continuous variables (Prasad et al. 2006, Beck et al. 2011, Mutanga et al. 2012, Gessner et al. 2013. Unfortunately, RF regression predictions often fail to encompass the full range of observed values (Roy and Larocque 2012) and are also biased against predicting absences (Savage et al. 2015). Some studies have used single regression trees to predict sage-steppe vegetation summaries (Homer et al. 2012, Xian et al. 2015. These predictions under-represent low and high values (see Fig. 5 in Homer et al. 2012). Recently, RF regression models have built rangewide maps of sage-steppe indicator variables. The biases against low and high values are also evident in these new maps (see Fig. 5 in Jones et al. 2018). In this section, we describe our workflow, applying bias corrections to RF regression models.
We fit RF regression models to plot data and geospatial predictors (Table 1) using the quan-tregForest function (Meinshausen 2017;Fig. 2a.ii). We set the number of explanatory variables tested for each regression tree split equal to 13 (the default) and used 1000 trees in each model. These choices were based on exploratory data analysis of four common species (Poa secunda, Bromus tectorum, Artemisia tridentata var. wyomingensis, and J. occidentalis).
After initial model fitting, we used two correction procedures, adjusting both range and abundance ( Fig. 2a.iii). We extracted two predictions from each RF ensemble, which were corrected and combined: (1) average value predicted by all regression trees in the ensemble and (2) an absence probability (the proportion of regression trees predicting zero). We applied a linear correction to the average-value prediction (Legendre 2018) describing the relationship between observations and predictions for those average-value predictions that were above zero. Corrected values were truncated at zero to avoid unrealistic negative values. After the linear correction, we estimated likely absences from the absence-probability prediction. We identified one absence-probability threshold for each model that balanced false-positive and falsenegative errors using the precision-recall f-measure (Sing et al. 2005), and used those thresholds to create binary layers. The corrected averagevalue predictions were constrained by the binary layers to yield absence-corrected abundance predictions ( Fig. 2a.iv).
Random forest nearest neighbor imputation.-We used a nearest neighbor imputation approach that leverages information extracted from RF models to build multivariate predictions (Crookston and Finley 2008). The forestry community has used multivariate maps to illustrate forest structure over large areas for many years (Tomppo 1991, Ohmann and Gregory 2002, Hudak et al. 2006. Imputation has been a mainstay for forest inventory and monitoring precisely because it generates unbiased multivariate predictions (when k = 1) that maintain the covariance structure of the input data (Tomppo 1991, Eskelson et al. 2009, Ohmann et al. 2011. Imputation models generate predictions in a different way from RF regression models. Random forest regression models summarize all predicted values from many regression trees within the RF ensemble, while imputation models summarize predicted values from chosen source observations. When only one source observation is used (k = 1), response variable values are simply the values present in the chosen source observation (our approach, illustrated in Fig. 2b.iv, v).
In its simplest form, a Euclidean distance matrix defined by the values of all explanatory variables is used to select source observations (kNN; Tomppo 1991). The distance matrix can also be more complex, sometimes constructed from ordination models (e.g., MSN, GNN; Moeur andStage 1995, Ohmann andGregory 2002), or from the nodes matrices of one or more random forest models (our approach, hereafter referred to as RFNN, using R package yaImpute; Crookston and Finley 2008).
Our distance matrix was constructed from two RF models predicting categorical vegetation summaries. One categorical y-variable indicated plant community composition, and the other represented vegetation structure (Fig. 2b.ii). The first was constructed from a species by cover matrix (Wisconsin, square-root transformed, R package: vegan; Oksanen et al. 2018). We constructed the classes with agglomerative hierarchical clustering, and Ward's linkage method, applied to a distance matrix derived from the species 9 cover matrix (R package: cluster; Maechler et al. 2018). The second y-variable was created with the same general approach, creating a distance matrix from a cover by life-form matrix, to inform the clustering approach described above. Each internal RF model contained 1000 classification trees. We tested six variables at each split in each tree, the default (Fig. 2b.iii). As with the RF regression, we used the full suite of geospatial predictors (Table 1).

Assessment
To compare RF regression and RFNN, we performed a series of univariate and multivariate assessments of model performance. For both methods, we extracted predictions of all response variables. For RFNN, we extracted the second nearest neighbor plot for all of our input plots, roughly equivalent to, but dramatically more computationally efficient than, a leave-one-out cross-validation strategies (Ohmann and Gregory 2002). For RF regression, we extracted model predictions for our input plots using the out of box predictions for each plot. All assessments were based on these extracted predictions.
Individual response variables.-We assessed the accuracy of the five GRSG habitat summary variables through regression modeling describing geometric mean functional relationships (GMFR; Ricker 1984, Riemann et al. 2010. We reported three statistics derived from the model predictions for the plot data: unsystematic agreement coefficient (AC_uns), systematic agreement coefficient (AC_sys), and agreement coefficient (AC). To characterize precision, AC_uns describes scatter, or noise, around the regression line between observations and predictions (analogous to R 2 ). To characterize bias, AC_sys describes how closely the modeled regression line matches a 1:1 relationship between observations and predictions. To assess overall performance, AC integrates information from the two statistics to give an indication of prediction quality. For AC_uns, AC_sys, and AC, values closer to one indicate less noise, less bias, and greater overall agreement, respectively. These metrics convey a wellrounded understanding of the nature of different model strengths and weaknesses in vegetation mapping predictions (Wilson et al. 2012. We also constructed summaries of our observations and model predictions across larger hexagons in order to estimate the map's accuracy at two spatial scales (Fig. 1), following the methods described by Riemann et al. 2010. Selected hexagons for assessment summaries contained at least three plots. For each assessment hexagon, we calculated the average value for the observations and predictions for all plots falling within that hexagon and then generate AC_uns, AC_sys, and AC metrics.
To assess the capacity of the differing models to predict species occurrence patterns, we calculated sensitivity, specificity, and true skill statistics (Fielding andBell 1997, Allouche et al. 2006) for presence and absence predictions of 51 common plant species. For all of the accuracy statistics reported, values approaching one indicate perfect model performance, while values near zero indicate model failures. We examined the differences in model performance for binary predictions with paired t tests comparing the suite of statistics grouped by the two modeling approaches across all 51 species.
Species richness. -Additionally, we assessed model performance at estimating the community-level attribute of species richness from the 51 species mentioned above. We calculated predicted and observed species richness by summing the number of species observed or predicted as present at each plot. Relationships between the observations and predictions were characterized by AC_uns, AC_sys, and AC.
Multivariate habitat classifications.-We explored how modeling technique affected the classification of habitat based on multiple predicted vegetation attributes (i.e., joint distributions) with density difference plots. These plots were constructed by calculating two 2-variable histograms (modeled and original data) and subtracting one from the other, allowing us to examine differences in the joint distributions of predictions. We chose variable pairings to highlight (1) the possibility of logical errors in the model predictions (i.e., tree cover < juniper cover) and (2) the relationship of three pairings of GRSG habitat variables as they relate to conditions that are favorable to GRSG habitat (variables included juniper, sagebrush, exotic annual grasses, and food forbs). For the GRSG habitat variables, we selected threshold values delineating likely habitat based on the literature (Table 2). To summarize the overall model performance in multivariate classification, we calculated the proportional difference (hereafter referred to as the relative percent difference, or RPD) between the observed (o) and predicted (p) number of plots meeting each criteria as (p À o) 9 o À 1 9 100.

Mapping
We also examined the effect of stacking all four habitat variable maps. We used the stacked maps to show differences in how each modeling approach might represent GRSG habitat, as defined by the thresholds for each variable. Each stacked habitat map (one for each modeling approach) indicated areas meeting all criteria for GRSG habitat. We summarized the stacked allcriteria maps across GRSG Priority Areas for Conservation (PAC, U.S. Fish and Wildlife Service 2013), calculating the proportion of each PAC that met all GRSG habitat criteria.
The two modeling approaches also differed in the percentage of plots with predicted values that met our criteria for GRSG-suitable habitat (Table 3). For single-variable, RFNN-based estimates of GRSG-relevant habitat summaries, discrepancies in estimated habitat percentage were quite small. The absolute difference between observed and predicted percentage of the observations that meet individual habitat criteria for RFNN was always less than 3%. Single-variable, RF regression-based estimates of GRSG habitat showed larger divergences from GRSG habitat estimates indicated by the observations. Most notable, for food forbs, the RF regression model showed an 8.35% increase over the observations in the percentage of plots with predictions meeting the habitat criterion. For the combined estimates of habitat, RFNN's estimate was 4.46% less than the plot-based estimate, while RF regression's estimate was 5.02% greater than the plotbased estimate.

Species presence/absence and richness
The two modeling approaches differed in their capacity for discerning species presence and absence at the plot and 50,000-ha scales (Fig. 4), with the greatest differences observed at the 50,000-ha scale. At the plot scale, paired t tests (2-tailed) indicated that plot scale true skill statistics (mean difference = 0.14, P < 0.001) and sensitivities (mean difference = 0.22, P < 0.001) were greater for RFNN than for RF regression, but specificities were slightly greater for RF  Sveum et al. (1998), Connelly et al. (2000) † We used a definition for this variable threshold from early drafts of the document cited. Currently, this variable is defined to identify habitat when the ratio of invasive annual grass to perennial grass is less than one.
‡ Literature supports the importance of food forbs for summer grouse habitat, but no threshold is currently indicated. We have adopted a low threshold to show areas where their presence is more than incidental. Fig. 3. Geometric mean functional relationships between observed and predicted values for five greater sagegrouse habitat descriptor variables, comparing the two modeling approaches. The solid gray lines show a 1:1 line, while the dotted blue lines show the GMFR regression line describing the relationship between the observations and predictions. Panel labels along the left margin apply to all panels in the row. Axis names for panel m apply to panels a through m, while axis names in panel q apply to q through t. regression than for RFNN (mean difference = À0.08, P < 0.001).
At the 50,000-ha scale, the same pattern held true. Random forest nearest neighbor was stronger overall (mean difference true skill statistics = 0.11, P < 0.001), and with respect to sensitivity (mean difference = 0.35, P < 0.001), but the RF regression models yielded higher specificities (mean difference = À0.25, P < 0.001).
Predictions of species richness showed a related pattern. At the plot scale, the aggregated prediction of species richness from the 51 RF regression models performed better in terms of AC and AC_uns than RFNN, but RF regression predictions were more biased than RFNN (RFNN: AC_sys = 1.00, RF regression: AC_sys = 0.49; Fig. 5a, b). RF regression had a tendency to underestimate species richness at both scales (Fig. 5b, d). At the 50,000-ha scale, the RF regression's prediction had a greater correlation (RFNN: AC_uns = 0.50, RF regression: AC_uns = 0.90), but RFNN prediction of species richness performed better overall than the RF regression prediction due to its lack of bias (Fig. 5c, d). Random forest regression consistently underestimated species richness.

Joint distributions and covariance in GRSG habitat variables
Our examination of the joint distributions began with a comparison of total tree canopy cover and tree cover attributed to western juniper. This analysis showed that RFNN more closely mirrored our observations than RF regression (Fig. 6). The two-dimensional histogram for the joint RFNN prediction showed frequency divergences from the two-dimensional histogram of the observations of À16 to 11, whereas RF regression prediction frequency divergences ranged from À65 to 17. Additionally, RFNN never predicted juniper > tree cover, whereas RF regression predicted juniper > tree cover on 5.8% of the plots (a logical error since we calculated tree cover as the sum of all species that potentially reach tree stature). Most of these logical errors were of small magnitude (95% of the errors were of <7% difference). Similarly, the range in frequency divergences shown in the two-dimensional difference histograms for GRSG  habitat variables for RFNN was less than that for RF regression, especially in cases where one of the two response variables equaled zero (i.e., along the x-axis or y-axis; Fig. 7).
When we aggregated these divergences based on GRSG habitat thresholds (Table 2), we found that RF regression overestimated and RFNN underestimated the number of plots meeting all GRSG habitat criteria (Table 3). The two-variable aggregated percent divergences ranged from À3.5% to 1.0% for RFNN and À1.4% to 10.1% for RF regression (Fig. 7). Despite the over-predictions by the variable pairings shown here, aggregated RF regression model predictions still showed a relatively small absolute difference in the plot-based estimates of GRSG habitat (RF regression classified an extra 1.3% plots as habitat than were present in the source observations, which is a 5.02% relative difference; Table 3). Random forest nearest neighbor, on the other hand, underestimated the percentage of plots meeting all habitat criteria (1.2% fewer plots predicted as habitat than were observed, which is a À4.46% relative difference; Table 3).

Mapped distributions of GRSG habitat
We found substantial differences between the maps generated by RFNN and RF regression for individual variables, and the stacked GRSG habitat maps. Maps of tree cover indicated broadscale agreement between the two methods. The RFNN and RF regression maps of tree cover <5% were the most similar of all habitat variables considered. Random forest nearest neighbor and RF regression maps indicated that 76% and 66% (respectively) of the landscape met this criterion (Fig. 8a, b). The two maps agreed on the tree criterion over 87% of the landscape (Fig. 8c). For the sagebrush maps, RFNN and RF regression indicated that 62% and 73% (respectively) of the Fig. 5. GMFR regression assessment of species richness predictions at the plot and 50,000-ha scale for random forest nearest neighbor (RFNN) (a, c), and stacked random forest (RF) regression models (b, d). Species richness calculations only include the 51 species that were modeled for this paper. landscape met the sagebrush criterion for GRSG habitat (Fig. 8d, e). The two maps were in agreement over 72% of the landscape (Fig. 8f). For the remaining 28%, lands were mapped as habitat by one approach, but not the other. These areas of disagreement were generally interspersed within the areas mapped as habitat for both models. The RFNN and RF regression maps of the areas with food forb cover >2.5% differed more noticeably than those describing the tree and sagebrush criteria. The RFNN map of food forbs showed small patches of habitat scattered throughout the landscape (49% of the map; Fig. 8g). The RF regression map showed larger regions where predicted food forb cover exceeded the habitat threshold (74% of the map; Fig. 8h). These larger regions were concentrated at higher elevations. The two maps agreed about the food forb habitat criterion over 67% of the landscape (Fig. 8i). The two mapping approaches also yielded differing pictures of exotic annual grasses in the landscape. Both maps show the majority of the region meeting the GRSG habitat criterion for exotic annual grasses (76% for RFNN; Fig. 8j, and 71% for RF regression; Fig. 8k). The two maps were in agreement about the exotic annual grass criterion over 77% of the landscape (Fig. 8l) as a whole. Differences were especially noticeable through the eastern half of the map where RF regression showed much less area meeting this criterion.
When the maps of all four GRSG habitat criteria were stacked, the resultant habitat maps showed different amounts of habitat present in the landscape (20%, and 25% for RFNN and RF regression, respectively), and it was mapped to noticeably different locations (Fig. 9). The stacked Fig. 6. Joint-difference histograms comparing original data, and modeled predictions for total tree cover and western juniper cover for RFNN (a), and RF regression (b). Legend indicates the difference in the number of plots that were predicted to each cell in the joint ranges from the predicted data to the original data. Negative values (red) indicate under-prediction within a cell, while positive values (blue) indicate over-prediction in a cell. Numbers indicate the percentage of model predictions that fall within each half of the graph, as delineated by the 1:1 dotted line.  (Table 2). Relative percent difference (RPD) numbers within boxes indicate the percent difference between the number of plots predicted and observed that fall within the delineated variable-space rectangle. GRSG habitat maps for RFNN and RF were in agreement over 76% of the map. In the northeastern corner of the region, RFNN mapped far more habitat than did RF regression, while in the northwestern quarter of the map, RF regression indicated more habitat than did RFNN.
In general, both RF regression and RFNN map the most GRSG habitat in the PACs located in Fig. 8. Single-variable maps describing potential greater sage-grouse habitat for random forest nearest neighbor (RFNN; yellow; a, d, g, j), random forest regression (RF; blue; b, e, h, k), and both (green; c, f, i, l). the northeast, northwest, and southwest corners of the region. The RF regression-based habitat map consistently showed more GRSG habitat throughout all of the PACs (average absolute difference: 14%). Differences in the spatial distribution of mapped habitat between the two methods were highlighted by zonal summaries of the proportion of area within each PAC that met all four habitat criteria in the maps (Fig. 10).

DISCUSSION
This study highlighted a common challenge for modeling vegetation to support applications where multiple descriptor variables are needed: Independently modeled vegetation attributes inherently ignore the correlated nature of vegetation attributes. The correlated nature of vegetation is relevant to many practical uses of the vegetation map, such as habitat mapping as we have illustrated here. We further discuss the implications of small attribute biases, and errors in attribute covariance in the context of mapping to support conservation planning and land-management policy.
When single variables were considered independently, multivariate imputation (RFNN) often out-performed univariate models (RF regression) by a small to modest margin at fine spatial scales (Fig. 3). The small biases identified for RF Fig. 9. Unified map showing areas that meet all four criteria from Fig. 8 according to each model type. The combination of yellow and green shows areas highlighted by random forest nearest neighbor (RFNN) as potential habitat in panels a, d, g, and j from Fig. 8. The combination of blue and green shows areas highlighted as potential habitat in panels b, e, h, and k in Fig. 8.   Fig. 10. Proportion of each Priority Areas for Conservation that is potential greater sage-grouse habitat, based on the four maps in Fig. 8, as mapped by random forest nearest neighbor (RFNN) (a), random forest (RF) regression (b), and the difference between the two estimates (RFNN À RF regression, c). regression became important when they coincided with GRSG habitat threshold values. This was most evident in the food forb variable. The RF regression models often inflated small values enough to cross the habitat threshold. This pattern may have stemmed from a combination of three factors: (1) The statistical distribution of raw values for this variable (zero-inflated, with a few large values), (2) missing relevant explanatory variables (e.g., competition, facilitation, or very recent fluctuations in rainfall), and (3) the resultant weak correlation in the raw model prediction may have performed poorly with our correction procedure. Regardless of the cause of the bias, the net effect is that the RF regression binary prediction of habitat based on the food forb variable yielded an inflated area estimate for suitable habitat in this dimension. When RF regression model predictions were combined, covariance errors often emerged in the joint predictions. This did not happen with the RFNN predictions. The difference between the two techniques with respect to error aggregation was most effectively shown with the example of species richness. Species richness was underestimated by the stacked RF regression models, but not by RFNN (Fig. 5). Aggregation errors were also evident in occasional illogical joint predictions of tree and juniper cover (juniper cover > tree cover) for the RF regression model, while this error never occurred in the RFNN model (Fig. 6).
Despite similar point-based accuracy metrics for broadscale model performance (Fig. 3), we still observed noticeable differences between the two modeling approaches in the broadscale patterns that emerged in three of the four GRSG habitat variables (Table 3, Fig. 8f, i, l). Unsurprisingly, those discrepancies were also evident in the stacked GRSG habitat maps (Table 3, Fig. 9). These trade-offs (precision, accuracy, threshold performance, and the capacity for post-modeling variable recombination) hold implications for the utility of vegetation maps for different applications. Maps of sage-steppe vegetation are needed to support a variety of activities, including habitat suitability modeling, conservation reserve design, and rangewide conservation programs for GRSG. Features that define a good map vary from application to application. Recent advances in change detection mapping for rangelands at very large spatial scales are very encouraging (Jones et al. 2018), but careful consideration of the error and bias structure is vital for determining how maps such as these can inform a planning and regulatory framework.
For vegetation data underlying habitat modeling, strong correlations between observations and predictions can be quite important (Hines et al. 2005). This need is balanced by other priorities, such as the availability of relevant variables (Johnson and Gillingham 2005) and accurate depiction of spatial configuration (Johnson andGillingham 2005, Rondinini et al. 2006) at appropriate spatial scales (Graf et al. 2005, Timm et al. 2016. For example, much work has been done in PNW forests to describe habitat for the Northern Spotted Owl, where the array of variables available in multivariate maps of forest condition have supported several studies mapping habitat across the region (McComb et al. 2002, Davis et al. 2011, Ackers et al. 2015. Published studies of regional-scale GRSG habitat mapping have relied heavily on classified vegetation maps such as LANDFIRE (Knick et al. 2013), on unclassified remotely sensed imagery (Donnelly et al. 2016), or on stacked univariate layers (Fedy et al. 2014). The multivariate approach to vegetation mapping that we used here expands the vegetation summaries that can be made available as explanatory variables to inform habitat models. For example, currently available broadscale data include information on tree and sagebrush cover (Xian et al. 2015), 250-m maps of herbaceous annual cover (grass; Boyte and Wylie 2018), and most recently, 30-m maps of percent cover of annual forbs and grasses, perennial forbs and grasses, shrubs, and bare ground over the western United States (Jones et al. 2018). In contrast, our joint modeling with RFNN can provide increasingly refined structural and compositional data. Here, we presented just five vegetation summary variables and diversity derived from presence/absence predictions for 51 plant taxa. However, other indicators of rangeland health (e.g., those described by Pyke et al. 2002) can also be summarized and mapped to the RFNN raster layer without the need for generating a new model and map. New variables can be formulated, shown spatially, and assessed independently of the RFNN modeling and mapping process. Covariance between new and old variables will remain consistent with the input plot data. The current study reiterates the importance of joint vegetation modeling for applications requiring species-or taxon-specific information to describe habitat, as is the case for GRSG. Our multivariate maps do contain a fair amount of fine-scale noise, a property shared with their forested cousins (Bell et al. 2015), but that noise dampens at broader spatial scales Fig. 5). Greater sage-grouse respond to vegetation at local and landscape spatial scales (Doherty et al. 2010, Arkle et al. 2014, and thus map performance across multiple scales merits consideration. For conservation reserve design, geographically marking variations in composition and diversity may be important (e.g., conservation of diversity hotspots, Myers et al. 2000). Species richness is a community-level metric that can be useful in this context. It is also a useful avenue for illustrating the hazards of aggregating predictions from many models. Our finding that stacked univariate models underestimated diversity is unusual in the literature. Many studies that map diversity by stacking multiple species distribution models find that species richness is overestimated (Dubuis et al. 2011, Guisan and Rahbek 2011, Pottier et al. 2013), a pattern we observed in predictions produced by stacked, uncorrected RF regression models (data not shown). In this paper, the under-prediction of species richness with the RF models was likely due to a consistent bias in favor of species' absences in our RF regression model predictions (indicated by lower model sensitivity; Fig. 4), which stems from our bias correction procedure. The binary correction that we used to generate our responses was unbiased, but the model II regression correction that we used to stretch the range of predicted values introduced a few additional absences to each prediction. Compared to stacked univariate models that may substantially over-or underpredict species richness, the multivariate RFNN approach can provide a less biased estimation of diversity.

Habitat assessment and adaptive management for GRSG in Oregon
The work we present here was conducted while building an imputed vegetation map to support GRSG conservation efforts in the state of Oregon. In Oregon, GRSG conservation efforts and plans have been facilitated by the Sage-Grouse Conservation Partnership, a collaborative group that includes the BLM and other federal agencies, state agencies, county-level entities, and conservation organizations and industry groups as well as other stakeholders. Through this partnership, the state adopted the Oregon Sage-Grouse Action Plan (Sage-Grouse Conservation Partnership 2015), which outlined a vision for cross-ownership cooperation in the state related to sage-grouse conservation, including assessing and mapping GRSG habitat. The Habitat Assessment Framework, developed by the BLM, and the western association of fisheries and wildlife agencies, was designed to facilitate understanding GRSG habitat across four spatial scales of understanding: species range, populations, home-range, and seasonal habitats/daily movements (Stiver et al. 2015). The BLM Approved Resource Management Plan outlines additional strategies for monitoring and managing GRSG habitat on Oregon's federal land (U.S. Bureau of Land Management 2015). All of these efforts rely on high-quality GIS data describing GRSG habitat, and effective management is hindered when these datasets are unavailable. Several relevant variables are not used in federal frameworks for addressing GRSG habitat because consistent and accurate maps describing them are unavailable across the GRSG's full range. Instead, derivatives from classified vegetation maps are often used (Chambers et al. 2017).
Our work is unique in that it provides detailed habitat information that is relevant to seasonal habitats in a spatially consistent way across the full range of GRSG in Oregon. BLM's adaptive management framework requires information from just two vegetation attributes within management areas for making strategic decisions (sagebrush canopy cover and tree canopy cover). The Oregon Sage-grouse Action Plan calls for the use of additional layers of information (e.g., exotic annual grasses and native grasses; Sage-Grouse Conservation Partnership 2015), because they are available in an imputed dataset much like the one created for the research presented here (created by the primary author to support region-wide modeling of landscape resources; Halofsky et al. 2014, Creutzburg et al. 2015. In both the rangewide and the within-Oregon context of GRSG monitoring and conservation, quantitative data layers are generally thresholdtransformed to binary and then combined into a single layer to provide information for making decisions. For example, the BLM defines their adaptive management triggers as multivariate, threshold-based decisions that are based on two data layers (e.g., a soft trigger is described: "When the area with at least 5% sagebrush canopy cover and less than 5% tree canopy cover drops below 65% of the sagebrush capable area within an individual Oregon PAC, but remains above 30%," U.S. Bureau of Land Management 2015: Appendix D). When conditions described by a trigger are met, management and conservation measures that are prescribed to the land in question become more protective of GRSG habitat.
Our analysis shows how aggregations of threshold-defined data layers can be quite sensitive to small response variable biases (e.g., food forbs and exotic annual grasses see Table 3 and Fig. 7e, f), yielding highly biased binary maps. When biased binary maps are combined, those biases accumulate. These discrepancies were clearly manifest at the scale of GRSG PAC, especially in terms of habitat amount (Figs. 9, 10). The differences in the aggregate map predictions were likely driven by the difference in the two mapped predictions for food forbs and exotic annual grasses ( Fig. 8g-l), our weakest variables. Reliability of combined estimates of habitat location and quantity are crucial components of adaptive management or regulatory frameworks.

CONCLUSIONS
As a foundation for wildlife habitat mapping and species conservation planning, the multivariate RFNN approach had clear advantages over univariate RF regression. At fine spatial scales, for single-variable predictions, the two approaches had comparable performance. At coarser spatial scales, both approaches had comparable precision, but RF regression often contained more bias. Recombination errors were consistently observed in stacked RF continuous regression maps, but not RFNN. For food forbs and exotic annual grasses, the most biased portions of predicted values from RF regression corresponded with habitat definition thresholds, yielding biased predictions of habitat quantity that were amplified when the two variables were combined. Because of its performance with respect to bias and the absence of covariance errors, we conclude that RFNN is wellsuited to inform mid-to broadscale strategic conservation/management decisions.