Pine and mixed‐conifer forests at low/mid‐elevations, where historical fires were relatively frequent, are broadly distributed across several ecoregions in the western United States (Fig. 1 ; Appendix S1 : Table S1). Although ponderosa pine often dominates these forests, they can also include Jeffrey pine ( Pinus jeffreyi ), which in places intermix with, and are similar to, ponderosa pine forests, and Madrean pine–oak ( Quercus spp.) forests with a diversity of pines. Mixed‐conifer forests at low/mid‐elevations are also broadly distributed across multiple ecoregions (Fig. 1 ). They can include additional pines (e.g., lodgepole pine, Pinus contorta ; sugar pine, Pinus lambertiana ), true firs ( Abies spp.), Douglas‐fir ( Pseudotsuga menzeisii ), and incense‐cedar ( Calocedrus decurrens ).

However, the fundamental premise for this fire management strategy has not been rigorously tested across broad regions. We broadly assessed the influence of forest protection levels on fire severity in pine and mixed‐conifer forests of the western United States with relatively frequent‐fire regimes to test this assumption. We used vegetation burn severity data from all fires >405 ha over a three‐decade period, 1984–2014, in forests with varying levels of protection.

It is a widely held assumption among federal land management agencies and others that a lack of active forest management of some federal forestlands—especially within relatively frequent‐fire forest types such as ponderosa pine ( Pinus ponderosa ) and mixed conifers—is associated with higher levels of fire severity when wildland fires occur (USDA Forest Service 2004 , 2014 , 2015 , 2016 ). This prevailing forest/fire management hypothesis assumes that forests with higher levels of protection, and therefore less logging, will burn more intensely due to higher fuel loads and forest density. Recommendations have been made to increase logging as fuel reduction and decrease forest protections before wildland fire can be more extensively reintroduced on the landscape after decades of fire suppression (USDA Forest Service 2004 , 2014 , 2015 , 2016 ). The concern follows that, in the absence of such a shift in forest management, fires are burning too severely and may adversely affect forest resilience (North et al. 2009 , 2015 , Stephens et al. 2013 , 2015 , Hessburg 2016 ). Nearly every fire season, the United States Congress introduces forest management legislation based on this view and aimed at increasing mechanical fuel treatments via intensive logging and weakened forest protections.

Methods

We used Gap Analysis Program (GAP) protection classes (USGS 2012), as described below, to determine whether areas with the most protection (i.e., GAP1 and GAP2) had a tendency to burn more severely than areas where intensive management is allowed (i.e., GAP3 and GAP4). We compared satellite‐derived burn severity data for 1500 fires affecting 9.5 million hectares from years for which there were available data (1984–2014) among four different forest protection levels (Fig. 1), accounting for variation in topography and climate. We analyzed fires within relatively frequent‐fire forest types comprised of pine and mixed‐conifer forests mainly because these are the predominant forest types at low to mid‐elevations in the western United States, there is a large data set on fire occurrence, and they have been a major concern of land managers for some time due to decades of fire suppression. We defined geographic extent of forest types from the Biophysical Settings data set (BpS) (Rollins 2009; public communication, http://www.landfire.gov) that derived forest maps from satellite imagery and represents plant communities based on NatureServe's Ecological Systems classification. Baker (2015) noted that some previous work found ~65% classification accuracy of this system with regard to specific forest types and, accordingly, he analyzed groups of related forest types in order to improve accuracy. We followed his approach (see Appendix S1: Table S1). The categories selected from the Biophysical Settings map were ponderosa/Jeffrey pine and mixed‐conifer forest types with relatively frequent‐fire regimes (e.g., Swetnam and Baisan 1996, Taylor and Skinner 1998, Schoennagel et al. 2004, Stephens and Collins 2004, Sherriff et al. 2014), compared to other forest types with different fire regimes such as high‐elevation forests and many coastal forests not studied herein. Forest types in our study totaled 29.2 million hectares (Fig. 1; Appendix S1: Table S1). We used the BpS data to capture areas that were classified as forests before fire, because postfire vegetation maps can potentially show these same areas as temporarily changed to other vegetation types. We sampled our response and predictor variables on an evenly spaced 90 × 90 m grid within these forest types using ArcMap 10.3 (ESRI 2014). This created a data set of 5,580,435 independent observations from which we drew our random samples to create our models. The 90‐m spacing was chosen because it was the smallest spacing of points that was computationally practical with which to operate.

Fires The Monitoring Trends in Burn Severity project (MTBS, public communication, http://www.mtbs.gov) is a U.S. Department of Interior and Department of Agriculture‐sponsored program that has compiled burn severity data from satellite imagery, which became available in 1984, for fires >405 ha, and was current up to 2014 (Eidenshink et al. 2007). The MTBS Web site allows bulk download of spatial products that include two closely related indices of burn severity: differenced normalized burn ratio (dNBR) (Key and Benson 2006) and relative differenced normalized burn ratio (RdNBR) (Miller and Thode 2007). Both indices are calculated from Landsat TM and ETM satellite imagery of reflected light from the earth's surface at infrared wavelengths from before and after fire to measure associated changes in vegetation cover and soil characteristics. We defined burn severity with the RdNBR index because it adjusts for prefire conditions at each pixel and provides a more consistent measure of burn severity than dNBR when studying broad geographic regions with many different vegetation types (Miller et al. 2009a, Norton et al. 2009). RdNBR values typically range from negative 500 to 1500 with values further away from zero representing greater change from prefire conditions. Negative values represent vegetation growth and positive values increasing levels of overstory vegetation mortality. The RdNBR values could be used to classify fires into discrete burn severity classes of low, medium, and high but this was not performed in our study, as we desired to have a continuous response variable in our models. We intersected forest sampling points with fire perimeters downloaded from MTBS to determine fires that occurred in our analysis area, and censored fires with <100 sampling points (81 ha). The remaining points represented sampling locations from 2069 fires (Fig. 1). We extracted RdNBR values at each sampling point as our response variable as well as predictor variables that included topography, geography, climate, and GAP status. These sampling points were used to investigate the relationship between forest protection levels and burn severity (Appendix S1: Tables S2 and S3). We chose topographic and climatic variables based on previous studies that quantified the relationship between burn severity, topography, and climate (Dillon et al. 2011, Kane et al. 2015).

Topographic and climatic data To account for the effects of topographic and climatic variability, we derived several topographic indices (Appendix S1: Table S2) from seamless elevation data (public communication, http://www.landfire.gov/topographic.php) downscaled to 90‐m2 spatial resolution due to computational limits when intersecting sampling points. These indices capture categories of topography, including percentage slope, surface complexity, slope position, and several temperature and moisture metrics derived from aspect and slope position. We used the Geomorphometry and Gradient Metrics Toolbox version 2.0 (public communication, http://evansmurphy.wix.com/evansspatial) to compute these metrics. We also computed several temperature and precipitation variables (Appendix S1: Table S3) by downloading climatic conditions for each month from 1984 to 2014 from the PRISM climate group (public communication, http://prism.oregonstate.edu). Climate grids record precipitation and minimum, mean, and maximum temperature at a 4‐km grid scale created by interpolating data from over 10,000 weather stations. To determine the departure from average conditions, we subtracted each climate grid by its 30‐yr mean monthly value. These “30‐yr Normals” data sets were also downloaded from the PRISM Web site and reflected the mean values from the most recent full decades (1981–2010). We determined mean seasonal values with summer defined as the mean of July, August, and September of the year before a given fire; fall being the mean of October, November, and December of the previous year; winter the mean of January, February, and March of the current year of a given fire; and spring the mean of April, May, and June of the current year.

Protected area status and ecoregion classification We used the Protected Areas Database of the United States (PAD‐US; USGS 2012) to determine forest protection status, which is the U.S. official inventory of protected open space. The PAD‐US includes all federal and most State conservation lands and classifies these areas with a GAP ranking code (see map at: http://gis1.usgs.gov/csas/gap/viewer/padus/Map.aspx). The GAP status code (herein referred to interchangeably as GAP class or protection status) is a metric of management to conserve biodiversity with four relative categories. GAP1 is protected lands managed for biodiversity where disturbance events (e.g., fires) are generally allowed to proceed naturally. These lands include national parks, wilderness areas, and national wildlife refuges. GAP2 is protected lands managed for biodiversity where disturbance events are often suppressed. They include state parks and national monuments, as well as a small number of wilderness areas and national parks with different management from GAP1. GAP3 is lands managed for multiple uses and are subjected to logging. Most of these areas consist of non‐wilderness USDA Forest Service and U.S. Department of Interior Bureau of Land Management lands as well as state trust lands. GAP4 is lands with no mandate for protection such as tribal, military, and private lands. GAP status is relevant to the intensity of both current and past managements. We made one modification to GAP levels by converting Inventoried Roadless Areas (IRAs) from the 2001 Roadless Area Conservation Rule (S_USA.RoadlessArea_2001, public communication, http://data.fs.usda.gov/geodata/edw/datasets.php) to GAP2 unless these areas already were defined as GAP1. We considered most IRAs as GAP2 given they are prone to policy changes and because they allow for certain limited types of logging (e.g., removal of predominately small trees for fuel reduction in some circumstances). However, we note that very little logging has occurred within IRAs since the Roadless Rule, although there occasionally have been proposals to log portions of some IRAs pre‐ and postfire, and fire suppression often occurs. We modified level III ecoregions (U.S. Environmental Protection Agency (EPA) 2013) to create areas of similar climate and geography (Fig. 1). We did this by extracting ecoregions and combining adjacent provinces in our study region.

Random Forests analysis We investigated the relationship between protection status and burn severity using the data‐mining algorithm Random Forests (RF) (Breiman 2001) with the “randomForestSRC” add‐in package (Ishwaran and Kogalur 2016) in R (R Core Team 2013). This algorithm is an extension of classification and regression trees (CART) (Breiman et al. 1984) that recursively partitions observations into groups based on binary rule splits of the predictor variables. The main advantage of using RF in our study is that it can work with spatially autocorrelated data (Cutler et al. 2007). It can also model complex, nonlinear relationships among variables, makes no assumption of variable distributions (Kane et al. 2015), and produces accurate predictions without over‐fitting the available data (Breiman 2001). Our independent observations were a random subset of our 5.5 million points, from which we drew three random samples of 25,000 points each. Each sample consisted of 500 fires randomly selected without replacement from the pool of 2069 fires. Fifty points were then randomly selected within each of the 500 fires. Our dependent variables were all continuous (Appendix S1: Tables S2 and S3) except for the main variable of interest, protected area status, which included the four GAP levels. The three observation samples were used to create three RF model runs, each consisting of 1000 regression trees. We conducted three RF model runs to assess whether our random samples of 25,000 points produced fairly consistent results. The RF algorithm samples approximately 66% of the data to build the regression trees, and the remaining data are used for validation and to assess variable importance. We used this validation sample to determine the amount of variance explained and variable importance. The algorithm also produces individual variable importance measures by calculating differences in prediction mean‐square‐error before and after randomly permuting each dependent variable's values. Variable importance is a measure of how much each variable contributes to the model's overall predicative accuracy. Unlike linear models, RF does not produce regression coefficients to examine how a change in a predictor variable affects the response variable. The analogy to this in RF is the partial dependence plot which is a graphical depiction of how the response will change with a single predictor while averaging out the effects of the other predictors, such as the climatic and topographic variables (Cutler et al. 2007). We used this approach, in addition to using RF to determine overall variable importance as described above, in order to determine the effect of GAP status, in particular, on fire severity, while averaging out effects of climate and topography.

Mixed‐effects analysis We performed a linear mixed‐effects analysis using the “nlme” add‐on package in R (Pinheiro et al. 2015). We used a random intercept model and identified year of fire (n = 31) and ecoregion (n = 10) as random effects. Similar to our RF models, our independent observations were a random subset of our 5.5 million points but for these models we drew three random samples of 50,000 points each. Each sample consisted of 500 fires randomly selected without replacement, and within each of those fires, 100 points were randomly selected. Our dependent variables were the same used in our RF models, and we log‐transformed the non‐normal variables of slope, surface roughness, and topographic radiation aspect index. We removed dependent variables that were correlated with each other (Pearson's r > 0.5), retaining 21 of 45 candidate dependent variables, and centered these on their means. Model reduction was performed in a stepwise process using bidirectional elimination with Bayesian information criterion selection criterion.