US National Maps Attributing Forest Change: 1986–2010

: National monitoring of forestlands and the processes causing canopy cover loss, be they abrupt or gradual, partial or stand clearing, temporary (disturbance) or persisting (deforestation), are necessary at fine scales to inform management, science and policy. This study utilizes the Landsat archive and an ensemble of disturbance algorithms to produce maps attributing event type and timing to > 258 million ha of contiguous Unites States forested ecosystems (1986–2010). Nationally, 75.95 million forest ha (759,531 km 2 ) experienced change, with 80.6% attributed to removals, 12.4% to wildfire, 4.7% to stress and 2.2% to conversion. Between regions, the relative amounts and rates of removals, wildfire, stress and conversion varied substantially. The removal class had 82.3% (0.01 S.E.) user’s and 72.2% (0.02 S.E.) producer’s accuracy. A survey of available national attribution datasets, from the data user’s perspective, of scale, relevant processes and ecological depth suggests knowledge gaps remain. attribution event types to forest canopy cover loss new dataset, the American Dynamics Attribution (NAFD-ATT) data, attributes the event types event, this major types of fast

Datasets produced to inform scientists, land managers and policy makers about forest status and trends have progressed over time with technology and user needs. Historically, statistical inventories, such as the US Department of Agriculture's Forest Service's (USFS) Forest Inventory Analysis (FIA) and the National Resource Conservation Service's (NRCS) National Resource Inventory (NRI), provided the nation's forest monitoring and land use change data, respectively [41,42]. These inventory programs, once decadal, have responded to needs for statistically sound estimations derived from more frequent sample periods [43,44]. The USFS and NRCS have always used high-resolution aerial imagery to view plot locations at the time of ground-based sampling, to help reduce costs and improve efficiencies.
A proliferation of national-and continental-scale, satellite-based forest maps have become available, due to the legacy of moderate resolution imagery from spaceborne sensors such as Landsat, free public access to the Landsat archive, and the availability of high-end computing capabilities [45][46][47]. Some mapping approaches are more manually intensive, balancing expert knowledge with automatic algorithms and expending large efforts on final cartographic cleanliness [48,49], while others put a high premiums on full automation [50,51]. Some national maps of forest change are produced from single collections (a collection of imagery from a sensor processed with a single method through time) using multi-year time steps (4+ years between sequential image observations) [52,53], single collections for annual (or higher step) time series stack analysis [54,55], and multiple collections for multi-year time steps [49,56,57]. In addition, each change detection algorithm has its own merits and limitations [58][59][60]. In dense time series image stacks, abrupt changes in forest cover are aptly found by change algorithms focused on statistical anomalies [61,62], while slower and more subtle changes are more readily captured by segment-or trajectories-based algorithms [63,64]. Algorithms based on fitting curves to the time series signal and analyzing the residuals can capture subtler, faster changes between and within seasons [65][66][67][68]. The capturing of a change and the attribution of a causal process are often achieved at different steps in the workflow or processing stream.
This paper focuses on research into the attribution of event types to forest canopy cover loss events in the contiguous US. We introduce a new spatial dataset, the North American Forest Dynamics -Attribution (NAFD-ATT) data, which attributes the most common event types (no event, removal, fire, stress, conversion and wind) to forests in the conterminous US (CONUS) between 1986 and 2010. A goal of this study was to employ the outputs from multiple forest change algorithms to differentiate all major types of fast (removals, fire, wind and conversion) and slow (stress) disturbances, and to differentiate between temporary (disturbance) and persistent canopy loss (deforestation) processes. The North American Forest Dynamics (NAFD) project was a long-term collaboration between the University of Maryland, NASA and the US Forest Service, Forest Inventory and Analysis (USFS FIA) program [54,69,70]. The primary NAFD goal was to increase the spatial and temporal accuracy of national-scale assessments of forest disturbance dynamics, for carbon cycle science under NASA's North American Carbon Program. Table 1 lists available CONUS attribution datasets, highlighting salient differences between the products from a user's perspective of scale, processes and metrics of ecological significance. It is known that areal measures derived from estimation and mapping do not always coincide [71]. Table  1 highlights estimated products (light grey), mapped products (dark grey) and hybrid products (medium grey). To evaluate which dataset might be more appropriate for their needs, based on possible reporting and analysis scales, users can look under the temporal and spatial row tabs in Table 1. Referencing this table, data users can easily look to see if a data product captures withinclass changes, such as degradation or partial forest removals, or a particular process of interest. Monitoring Trends in Burn Severity (MTBS) attributes and maps only one process: annual fire events on all lands [72]. A testament to its efficacy in mapping this one process, MTBS data is ingested by most of the other projects that attribute fire, including NAFD-ATT [50,53,73,74]. Other users may want to focus on ecologically relevant metrics to monitor and describe disturbance regimes, such as frequency, succession, magnitude and duration [75,76]. Table 1. Publicly available the conterminous United States (CONUS) datasets attributing forest change processes. * denotes operational projects that are expected to maintain data legacy. ** denotes more than one process that was aggregated into a single type class for attribution. Y is short for "Yes"; N is short for "No".
The NAFD-ATT map and data products were created from statistically sampled training data, across all ownerships of forested land, and a nationally consistent machine learning approach across the CONUS, yielding 25 years of forest event type history. These products will be available from the Oak Ridge National Laboratory (ORNL) Distributed Active Archive Center (DAAC) (Available online (http://daac.ornl.gov) from ORNL DAAC, Oak Ridge, Tennessee, USA. https://doi.org/10.3334/ORNLDAAC/1799) and the USFS geospatial data platform. According to the results presented here, between 1986 and 2010, 70.6% of all CONUS forests remained stable without loss events. Temporary canopy loss from removals, fire and stress were attributed to the majority of canopy loss events (70.7%, 23.6%, 3.6%, and 1.4% of all forest land respectively) and only 0.6% of all forest land had canopy loss events attributed to land use/cover change (persisting forest loss). Further spatio-temporal results, accuracy measures, limitations and next steps are reported in the manuscript.

Training Data
We used a sample of 7200 plots across the CONUS, collected using TimeSync software for the North American Forest Dynamics (NAFD) project as the response data for the classification model [77]. A stratified two-stage sample design was used to collect the plots. The five strata used are coarse generalizations combining different ecoregions that share similar forest types. Selecting more samples in strata with more abundant and diverse forest increased the likelihood that response plots would fall on forestland and, by extension, that more forest change would be captured. Theissen polygons, calculated using the center points of the Landsat World Reference System path-row grid, comprised the non-overlapping sample units. Of the 442 Theissen polygon sample units that overlap CONUS boundaries, 180 were randomly selected from within the five strata. Within each sample unit, 40 plots were randomly distributed. A total of 2700 of the training plots fell on land that had forest cover (planted or naturally vegetated land likely to contain 10% or greater tree cover at some time during the time series sequence). Full details of the sample and response design of the training dataset and related visuals are published [73].
The plot observations were collected on each plot (30 × 30 pixel) using annual cloud-filled Landsat time series stacks (nominally from 1984 to 2012) within TimeSync software, ancillary datasets and high resolution aerial imagery in Google Earth where available [54]. Analysts temporally segment the data (image series) by labeling and recording visual observations of land cover, land use and change processes at the plot level [77]. Plots originally labeled "Mechanical" change processes and which occurred in forest cover were relabeled to "removals" or "conversion," depending on the land cover/use and timing observations in the plot records. At the end of the imagery time series, it is harder to determine with confidence whether a forest clearing activity was related to conversion or removals. We cautiously labeled all mechanical clearings as removals, if they occurred within six years of the end of the time-series stack, unless evidence of specific land use change was recorded.
In addition to these plots, further training data, gathered for the pilot study, included plots from the LANDFIRE project, FIA ground data and analyst interpretation of Landsat/aerial time-series images using TimeSync software [78]. These additional, opportunistically gathered data, not appropriate for inference, were not used in calculating (OOB) accuracy measures [79].

Predictor Variables
Forest disturbance algorithms have unique sets of costs and benefits in terms of errors omitted (false negatives) and committed (false positives) for each forest disturbance event type [58]. We used outputs from three change algorithms as predictor variables to mitigate omissions. The vegetation change tracker (VCT) algorithm excels with abrupt disturbance processes of varying magnitude [61,80]. The Shapes algorithm is a spline approach constrained by a priori defined spectral trajectory "behaviors", based on forest ecology principles, including long slow declines [81][82][83]. MTBS burn severity raster data reflect known fire events that meet a certain minimum size criterion [84]. The Shapes and VCT algorithms were run on NAFD Landsat time-series stacks including nominal imagery from 1984 to 2011. We subset all algorithm outputs to the years 1986-2010.
As part of the North American Forest Dynamics (NAFD) project, access to NASA's high end computing capability (HECC) and the NASA Earth Exchange (NEX) collaborative platform allowed for efficient national scale data prep, storage and processing on 25 years of Landsat data [85]. Over 60,000 Landsat images from the 1984-2011 time period (Landsats 4-7 Thematic Mapper) were evaluated for processing and analysis as part of the NAFD processing stream. Landsat image processing, completed using various Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) and VCT modules, included conversion to surface reflectance, atmospheric adjustment, individual scene masking (for clouds, shadow, snow, bad data and water), intra-annual growing season clear view compositing, vegetation index calculation, assembling and clipping individual year images into a time-series stack, and between-year data interpolation (in areas of no observations). After clear view compositing, 62% of the annual images in the time-series were clear-view, 27% were clear composites, 9% had remaining cloud contamination and 2% had no acceptable images. Full descriptions of image selection, quality and the image processing streams are available in previous publications [54,86].
We ran Shapes on the Landsat time-series stacks (see above) of Normalized Difference Vegetation Index (NDVI), band 5 surface reflectance, Forestness Index (FI) and Normalized Burn Ratio (NBR), yielding band specific temporal, magnitude and classified "shape" metrics for each band.
Annual VCT rasters, generated as part of the NAFD project and referred to here as NAFD-NEX maps,included annual maps of the year of disturbance onset and multiple magnitude indices. We converted VCT disturbance year rasters to annual polygon records and their calculated area, perimeter, shape and fractal indices. Using a "most recent" event per pixel rule, the annual polygon data and related metrics were composited into single date, geometry and magnitude layers. The same rule was applied to the annual magnitude rasters to create a composite layer for each metric.
Annual MTBS burn severity raster data were composited into a presence/absence layer using the low, moderate and high intensity values. In the case of more than one fire event in a pixel through time, the highest magnitude event date (year) was selected. The most recent event was selected in cases where multiple highest magnitude events occurred in the same location. Selected pixels were composited into a binary fire layer and a separate corresponding year of fire layer.
Local geographic predictors included aspect, slope, latitude and longitude, elevation and transformed aspect. Vegetation predictors include USFS forest group type and probability, down sampled from 250 to 30 m [87], and LANDFIRE Existing Vegetation Type (EVT v1.4.0) data collapsed to three classes: agriculture, forest and developed [88].

Modeling Process
A pilot study of 10 Landsat scene locations demonstrated that a two-step "flat-to-annual" process reliably separates multiple classes of canopy change activities [78]. In the first step, we construct a categorical Random Forest [79] for event type response classes using the R ModelMap package [89]. In the second step, disturbance year (and duration in the case of the stress category) were calculated using a rule-based approach and temporal metrics from the three predictor disturbance algorithms (see Section 2.2).
Modeling continuously across the nation brings different methodological and validation challenges from modeling in localized regions [90]. We tested building and predicting models over different geographic combinations of training data strata, aggregating response classes and down sampling to mitigate imbalanced training data. We visually compared predicted maps from model runs to published geospatial data sets on different process types to select the most representative model predictions. To provide insight into variable importance, we compared per class commission and omission from models constructed with different suites of training data.

Accuracy Assessment
We followed the best practices for robustly estimating map accuracy from Olofsson et al. [91] where possible. Reference data for multiple disturbance types is difficult to attain, so we used all training data to build the model, rather than hold out for validation. We used the R survey package [92,93] to calculate OOB errors for the stratified two-stage design for the five reference data strata. This yielded estimates of overall accuracy, user's accuracy (for commission errors) and producer's accuracy (for omission error) statistics, each with associated sampling variance, standard errors and confidence intervals for each response class. The analysis conducted was consistent with the sampling and response design of the training data. For additional information on the map quality, we created two confidence layers reflecting qualitative measures of model confidence and dominance. Random Forest creates a vector of values with proportions of votes from the forest of trees for each of the response classes and selects the class with the majority proportion of votes as the winning prediction class for each prediction unit (pixel). To document the strength or confidence of the model call, we subtracted the second highest vote proportion from the highest proportion (from the winning response class) [94]. The highest possible score is 100, representing the largest possible separation between the winning and the next common class, and therefore more model "confidence." We calculated Simpson's dominance index [95] by summing the squares of the vote proportions for each response class on each prediction unit (pixel). This dominance index provides a measure of the dominance of the winning class. Values closer to 100 reflect a higher dominance (of one class over the others). A decreased or lower dominance would be related to a related rise in the equitable distribution of the class proportions, especially if the number of classes with votes increases; while this type of diversity is desirable in natural systems, in our categorical map predictions it is not.

CONUS Map of Event Type
The spatial patterns of the most dominant types of canopy cover loss events are visible when mapped (Figure 1). Stable forest-forests with no loss event between 1986 and 2010-is the largest class. Though stable from a loss perspective, these forests may be in phases of regrowth, recovery or reforestation. Removals occurred in 612,200 km 2 (23.6%) of CONUS forests, with regional centers of large commercial forestry activity apparent even in national scale maps. Nearly 41% of the southern forest area had removal events over the 25-year period. Large clusters of fire activity are distinct even at this small scale (1:26,000,000), mostly in the West, with a few large fires visible along the eastern seaboard. Fire occurred in 94,483 km 2 (3.6%) of CONUS forests. The training data defined stress as any event resulting in slow gradual loss of forest canopy, including insect damage, drought and disease, which often co-occur. Most of the 35563 km 2 (1.4%) of CONUS forests affected by stress events were in the Rocky Mountain (RM) FIA region. Note, this area is calculated with onset-based (extensification) accounting, where an event is counted once at the year of onset/detection, rather than with duration (intensification), where an event is counted again in each year of the decline. Stress usually affects a narrow subset of tree species in the forested landscape, leaving a dispersed spatial patter that is difficult to visualize on a national map. The relatively rare (across regional to national landscapes) footprints of conversion and wind events are also hard to view on a map of this scale. The forestland area affected by each of these event types is 16,655 km year ⁄ and 166 km year ⁄ , respectively.

Annual Rates by Event Type
Annually, the prominence of a change type can vary widely within and between regions, though the inter-region dominance of each change class by region is mostly stable (Figure 2). For example, removals occur most in southern forests through the period of observation, though the annual rate of southern removals fluctuates a good deal. The magnitude, direction and variability of the southern removal rate (area year ⁄ ) drives trends in the national removal rate, and makes up more than half of the total national forest canopy loss rate in most years. The multi-year trends in the southern rate are echoed in the national removal rate. Noticeable are the peaks of activity in 1989, 1998 and 2006, the latter followed by a precipitous decline likely due to the economic recession of the late 2000s. By the end of the record (2010), the removal rates returned to almost 1986 levels. Northern forests contribute the second highest area per year to national numbers, rising in 1998 to the steady rate that remained through 2010.

Figure 2.
Forest area (km 2 ) affected annually by event type between 1986 and 2010. Annual rates are stacked and color filled by FIA region. Note that the scale of the y-axis varies per event type. A second y-axis on the right labels the area as a percentage of total CONUS forestland.
Annual removal rates reported here are lower than FIA decadal harvest estimates. FIA estimates report consistent harvest rates of 45,000 km year ⁄ (3500 km 2 is in interior Alaska), with variability in the proportions of "clear cut" and "partial harvest" events by region through time [96,97]. Tao et al. [98] used NAFD and FIA plot data to map the percentage of basal area removed, and found that the fluctuations in overall rates were driven primarily by variability in the rates of partial (< 80% basal area) removal events, while rates for clearing events (> 80% of basal area) were stable over the 25year record. Harvest intensity and management objectives, i.e., commercial versus sanitation harvest, are critical information for evaluating the relationship between harvest rates and social, economic and environmental drivers [70,99,100]. Likewise, information on the type and method of cut, which can help determine potential product lines, decay rates and residue ratios, is important for carbon cycling, but difficult to know from remote sensing data alone [101,102]. In areas with large commercial operations and/or quick stand rotations, where stands often undergo more than one treatment before the "final cut" or clearing, our approach may underestimate annual removal areas. For example, annual NAFD-NEX maps record that 15.3% of disturbed forest locations had two distinct canopy loss events in the 25-year record, and another 1.3% had more than two [80].
The National Interagency Coordination Center reports steep spikes in national wildfire acreage in 1988, 1990, 1994, 1996, 1999, 2000, 2002 and 2004-2007 [103], similar to the peak forest fire years in these data (Figure 2). Interior West forests have the highest variability in annual forest fire area, with multiple increases in successive years (i.e., 2004-2007). As expected for event types driven by large climatic patterns, forest fire rates in the neighboring Rocky Mountain (RM) and Pacific Coast (PAC) forests have similar fire ecology and peak fire years. Southern forest fire rates have different peak years than the peak years for western forest fires, suggesting different regional climatic drivers than in the West. The similarly timed crests in the overall national rates of canopy loss are evidence of the large areas affected in peak western forest fire years.
Counting locations with repeat fires (using low, medium and high severity pixels) from MTBS, we found 10.1% of locations had two fire events and 2.27% had more than two fire events (1986-2010). These MTBS estimates are for all lands, and so will be higher than for forest fire events only. Still, the composite approach used in this study, which only records one event through time, is likely to underestimate sequentially occurring canopy loss events.
Studies suggests that large tracts of forest have been in a state of declining health for decades [73,[104][105][106]. The annual rates for stress record a large spike at the beginning of the record, which is orders of magnitude larger than the mostly consistent rates in other years, and is almost all attributed to RM forests ( Figure 2). It is plausible that, in this study, the spike in observations of events characterized by long slow declines in canopy cover are correctly attributed. Most time-series analyses discount or discard the first and last year of analysis due to confusion between signal and noise in the first and last year of the time-series, due to spurious clouds or the weaker quality control measures, that often make change detection more difficult and increase false positives [94,[107][108][109]. This study likewise discarded the first and last year of the analyses to limit spurious false positives. The image stacks included nominal imagery between 1984 and 2011, but the analyses only report changes detected between 1986 and 2010. Further, the slow long duration decline "shape" in the Shapes algorithm, which is mostly associated with the stress class, cannot be "triggered" by a year or two of erroneous data, even at the beginning of the record of observation. The decline "shape" requires an accumulation of multiple years in the spectral signature [83]. Decline events that began before the period of observation in 1984-for example, an insect infestation or severe drought that began in 1980, and continued or intensified through the beginning of the time-series of observations-are simply attributed the first year label available, which in this case is 1986. It is important to note that 1986 may not be the peak year of the decline. Still, the challenges of correct temporal attribution (year of onset) are separate from the plausibility that the event is correctly attributed as decline and accurately reflects a steady slow decline in the spectral response of the forest canopy. The temporal requirements of the decline "shape" in the algorithm make true change less detectible (perhaps yielding more false negatives) at the end of the temporal record. In the case of slow decline events, problems of false negatives at the end of the time-series are a reality for all change detection algorithms. Training data for the stress class did not include abrupt insect mortality events, such as the widespread mortality from Southern Pine Beetle, so this class cannot be used for investigations into outbreaks of this endemic pest [110].
Accounting methods (see Section 3.1) for long duration canopy loss processes can produce widely different interpretations. This study mapped peak annual rates for the stress class of 0.3% of CONUS forest area. Cohen et al. [73] estimate two values for peak annual stress rates in 2001, depending on whether intensification or extensification accounting was used. When only one observation on the same plot was counted, regardless of the number of sequential years of decline, they estimated a peak rate of decline in 2011 of less than 0.5% of all CONUS forest area affected. When observations on the same plot were counted in each sequential year that the decline persisted, they estimated that the peak rate of decline (2011) occurred on roughly 3.0 % of all CONUS forest area. Though rates between these two studies are similar if the same temporal metric is used (0.3% of CONUS forest land vs. 0.5% of CONUS forest land), comparisons are confounded by differences in area of detected stress (numerator) and/or total forest area (denominator), which are difficult to tease apart when reported as ratios (percent of total forest area). Furthermore, the studies use different methods for calculating onset year and/or duration, and different regional reporting boundaries. This study used regional boundaries defined by FIA regions. Cohen et al. used regional reporting boundaries defined by their custom strata (see Section 2.1). Other studies offer mixed reporting methods for insect-related canopy loss through time, with the majority using onset accounting [111][112][113][114][115][116].
Rates of forestland conversion first show an increase from the beginning of the time-series to a peak in 1998, at 1,073 km year ⁄ (0.04% of CONUS forest year ⁄ ), then a steady and steep annual decline to 250 km 2 in 2010-half of the 1986 rate ( Figure 2). The timing and direction of these trends agree with the findings of both the attribution pilot study [117] and the NRI estimates, the latter inventory pertaining to mostly rural, non-Federal lands [118]. The NRI shows a peak rate in newly developed lands from forest lands of 3,804 ± 76 km year ⁄ (1992-1997 inventory), then a sharp decline to a rate of 1,133 ± 38 km year ⁄ (2007-2012 inventory), a little less than half of its estimated 1982-1987 annual rate.
Land use conversion and wind activities affect relatively small percentages of forest area at the national level, but locally are important drivers of persistent and temporary forest change, particularly in southern forests [4,31,[119][120][121][122][123][124]. Gathering training data on rare classes across large landscapes is difficult with statistical sampling designs, as is using/validating imbalanced data in statistical multi-class models [125][126][127][128][129]. Opportunistically sampled FIA data (see Section 2.1) yielded training data on mortality caused by Hurricane Hugo (1989) and the Boundary Waters-Canadian Derecho (1999). These data enabled the models to attribute plausible type and timing for these events, but not for extreme wind events, such as Hurricane Katrina (2005), outside of the local space and time neighborhood of the training. There is little national data available for wind-related forest damage in the West.

Accuracy, Confidence, and Variable Importance
User's and producer's accuracy are lowest for the stable, no loss class, with a national average of 81.3% (± 0.01 S.E.) and 94.4% (± 0.01 S.E.), respectively. The removal class performed similarly in the East and West models, with a weighted national average of 82.3% (± 0.01 S.E.) user's and 72.2% (± 0.02 S.E.) producer's accuracy. For the remaining classes, the models' predictive performance varied geographically, with higher class accuracies where the event types are most common and more training plots were available. In the fire-prone West, user and producer's fire accuracy are 91.2% (± 0.03 S.E.) and 70.0% (± 0.05 S.E.), respectively. Eastern forest fires, not including prescribed burns, which often occur under the canopy and are difficult to detect with Landsat, have lower accuracies: 41.2% user's (± 0.13 S.E.) and 31.4% producer's (± 0.1 S.E.). The stress class performs better in the West (66.7% user's (± 0.08 S.E.) and 26.9% producer's (± 0.05 S.E.)), again where the bulk of the training data occurs.
The attribution model conservatively predicts each event type compared to no change demonstrated by the number of combined accuracy metrics that graph above the 1:1 line (Figure 3). In other words, the results for each event type have fewer false positives (commission errors) than missing true positives (omission errors). The comparisons between this study's rates (area/time) for each class (not including stable) and those of other published studies (Section 3.2) support the finding that omission is higher than commission across classes, and suggest under-prediction of canopy loss by class. Due to the limited number of samples, we cannot make inferences about the accuracy of conversion in the West, nor about wind events. Differences in definitions and scales of observation of forest and forest loss events (in spatial, temporal or magnitude terms) can muddle accuracy metrics, by confounding errors of definition and observation with errors of measurement [130][131][132]. The training data had no minimum magnitude threshold, a minimum duration threshold of intra-seasonal loss (defoliation) and a minimum area threshold of a single tree (labeled as change if the analyst had visual evidence of canopy loss in a time-series of Landsat and/or high resolution aerial image chips). For the VCT and Shape algorithms, a minimum duration of two concurrent observations of change was required to detect change, to minimize false positives from clouds and spurious noise. The minimum Landsat-based canopy cover quantity, for reliable detection, may be between 10% and 20% canopy cover depending on environmental context, processing methods used and the definition of tree crown cover used in the training data [133,134]. The sensitivity of Landsat-based change detection algorithms to low magnitude or low intensity canopy cover loss events (percentages in either absolute or relative terms) is even more poorly characterized [54,58,135]. Low magnitude changes may have a lower accuracy than land clearing events. However, within (forest) classes, low magnitude cover changes will have smaller long-term implications or effects compared to high magnitude disturbance events, or high magnitude events that convert between forest cover and other land use/covers, which cause longterm shifts in ecosystem dynamics. . Models were built with "All" and separately with subsets or suites of predictor variables, to assess their individual contributions to omission and commission accuracy metrics for the major event type classes of removal, fire and stress. Subsets of predictor from "Local" geographic variables, MTBS disturbance variables, "Shapes" disturbance algorithm variables, the "VCT" disturbance algorithm variables and "Veg(etation)" variables are described in Section 2.2. User's Accuracy is 1 -Omission. Producer's Accuracy is 1 -Commission.
To evaluate variable importance, we plotted the class-based omissions and commissions for models built using all predictor variables and using suites of related groups of predictors for only local geography, vegetation, VCT, Shapes or MTBS groups of predictors (Figure 4 and see Section 2.2). While variable importance metrics are available for each individual predictor used in a Random Forest model, the value of these metrics are quickly diluted in the presence of correlation with other predictor variables, thus convoluting the variable importance analysis [136]. Model predictions built using all predictors had the lowest omission and commission errors in the removal and stress classes. Models built from VCT predictors alone had the second lowest omission and commission in the removal class. The model built from only MTBS predictors had slightly better (lower) omission and commission for fire than the model built with all predictors, but it had no predictive ability for other classes. For the stress class, the second lowest omission was achieved using a model built with the local geography group of predictors, but the Shapes group of predictors had the second lowest commission.

Conclusions
This study attributed both abrupt and slow forest loss event types, including removal (removals), fires, stress and conversion, annually across the nation, at the spatial scale of forest stands (hectares). The data and maps from this study used a unique combination of methods, including a consistent approach across the 25 years of Landsat data, an ensemble of predictors from multiple disturbance algorithms, and training/reference data sampled with an independent (not mapdependent) sampling scheme. Care was taken to harmonize definitions where possible, but residual errors remain due to differences in the definitions, the spatial support sizes of observations and the sensitivity (to the full range of canopy loss magnitude/severity) of reference, training and predictor data.
A 25-year record, in terms of forest disturbance regimes, is quite short [26,137]. This study found that over 25 years, CONUS forest area was dominated by stable forests (75%), and disturbance and conversion events were relatively rare (< 25%). We suggest that, overall, map accuracies are not a relevant metric where attribution of a minority (< 25%) class or event type is the goal. In such cases, class-level metrics yield more evidence of the method's efficacy. Here, all classes except stable had higher omission than commission, suggesting a conservative estimate of forest areas affected by canopy loss events due to fire, removal, stress and conversion. Stable, fire and removal classes had good class accuracies, with reasonably small sampling errors (for accuracy metrics). Very rare events in the CONUS, such as fire in the East, stress and conversion, had substantially lower training data amounts, accuracy metrics and sampling errors (for accuracy metrics). Wind events and conversion in the West were so infrequently captured in the training data that accuracy metrics could not be generated.
Sampling schemes designed to proportionally target rare events have some advantages when the primary purpose is categorical disturbance and attribution modeling [50,94,107]. The training data used in this categorical modeling effort was collected for a different purpose, using a sampling approach that did not explicitly target forest change events [54]. However, the study illustrates the importance of training data that can allow separation of disturbance and deforestation. A key component of this training data, and of many comprehensive forest monitoring programs, is the longterm remote sensing records, collected consistently at adequate scales, to observe fine-grain processes.
This national attribution map can inform sampling efforts, including pre-and post-stratification, to improve precision in design-based and small area estimation [117]. As with data of similar extent and grain, per pixel analysis is cautioned against [138]. Current applications of these data include moving window analyses, exploring relationships between gross forest loss, fragmentation and proportions of attributed change processes, and combining these observations with FIA plot attributes, such as ownership and forest type.
The survey of currently available CONUS attribution products (Table 1) suggests that a lack of ecological depth, and a difficulty in characterizing rare classes, persists. Using meta-analysis to combine attribution datasets has helped to evaluate the impacts of multiple processes on forest ecosystem services, such as carbon storage [116,[139][140][141]. The best methods for comparing and harmonizing multiple maps describing area and location of forest change events, in lieu of historical reference data, are under research [58,142]. Quantifying map uncertainty and its spatial characteristics would also aid user's ability to choose and/or combine maps that support their specific needs [135,143]. It is important to note that only operational national monitoring programs (* in Table  1) can provide data legacy, a key element for many users in natural resource management and decision support frameworks. The funding for this work ended in 2016, preventing further updates to the maps and methodology. However, the next generation of operational products, such as the USFS Landscape Change Monitoring System (LCMS), the Landscape Change Monitoring, Assessment and Projection (LCMAP), the LANDFIRE Remap (v2.0.0) and the National Land Cover Database (NLCD), promises to push the boundaries of national-scale ecological monitoring [48,49,64,144].