Estimation of Forest Disturbance from Retrospective Observations in a Broad-Scale Inventory

: Understanding the extent and timing of forest disturbances and their impacts is critical to formulating e ﬀ ective management and policy responses. Broad-scale inventory programs provide key estimates of forest parameters that indicate the extent and severity of disturbance impacts. Here, we review the use of a post-stratiﬁed estimator in a panelized design, in the context of disturbance observations that are collected retrospectively. We further develop a sample weight adjustment that is requisite for proper estimation of the extent and timing of disturbances. Using populations from areas of Arkansas, California, and Maine in the US, the weight adjustment technique was tested in a Monte Carlo simulation. We found that the estimated area of disturbance using the weight adjustment technique had satisfactory agreement with the true population values and performed considerably better than the conventional post-stratiﬁed estimation approach. The proliferation of panelized forest inventory designs globally suggests that accurate estimates of areal extent and timing of disturbances will often require that weighting adjustment techniques be employed in the estimation process.


Introduction
Disturbance is a universal component of forests and quantifying the areal extent and timing of forest disturbances is key to understanding and managing forested lands. Quality data on forest impacts from disturbance events aids in identifying management opportunities, restoration, and more accurately estimating effects on carbon stocks, sequestration rates, and other forest structure parameters. The primary causes of disturbance in the United States are fire, insects and diseases, drought, and harvesting [1], all of which have transformative effects on forest function and total biomass even at low levels of mortality [2]. Abiotic and biotic natural disturbance agents and interactions between them have resulted in amplified levels of forest disturbance [3][4][5], providing a pressing need for information on the extent, timing, and effects of forest disturbance. Likewise, there have been shifts in harvest frequencies over time. While some disturbance events are stand replacing, others result in selective mortality to a species or cohort. This research aims to develop improved techniques for disturbance estimation under rotating panel, repeated measure sampling designs for forest inventory, such as those used in the United States, France, and China [6].
Broad-scale forest disturbance is monitored through two main sources: the field-based component of national forest inventories (NFI) and remotely sensed imagery. Remotely sensed imagery has been used in forestry since the 1950s. More recent remote sensing applications include producing annual disturbance maps based on the Landsat time series stacks [7][8][9]. Since 1982, Landsat has collected spectral reflectance data at a 30 m spatial scale and 16-day temporal density, which is particularly applicable to forestry given the ability to detect changes across relatively short time intervals. With the application of robust models, these data allow for evaluating the timing and geographic extent of major forest disturbances. However, with the exceptions of Moisen and colleagues [10], Schleeweis and colleagues [11], and the Monitoring Trends in Burn Severity (MTBS) project [12,13], there has been little focus on developing models that identify the type of disturbance (e.g., drought, fire, harvest, etc.). Additionally, the impacts on forest structure and biomass are often not quantified in existing remote sensing approaches to disturbance assessment.
Forest inventory data from NFIs provide an alternative to remote sensing approaches or may be combined with remotely sensed data using model-assisted or model-based estimators. Through NFIs, in situ data on disturbance occurrence, disturbance type, and forest structure pre-and post-disturbance are collected, providing a valuable base for disturbance related analytics. However, it is common with longitudinal designs to record mortality-causing disturbances that have occurred since the previous plot measurement. This retrospective observation of disturbance, when coupled with rotating panel repeated measure designs, creates issues for estimating annual disturbance totals. To our knowledge, these issues are not addressed in the classic sampling and estimation literature offered by Cochran [14], Kish [15], and Särndal and colleagues [16], or the NFI specific literature (e.g., [17,18]). However, these estimation issues (described in detail in the following section) must be rectified, not only to construct correct annual estimates based on a direct estimator or model-assisted estimator, but also to construct correct annual estimates when employing small area estimation techniques where remotely sensed information serves as ancillary data (e.g., the Fay-Herriot area model [19]).
The objective of this research is to develop an appropriate approach for estimating the annual disturbance parameters such as areal extent and mortality volume when multiple panels of a field-based inventory in a rotating panel sample design are combined for estimation. To accomplish this objective, we first develop a technique to modify the sample weights of the rotating panel design to account for the retrospective observation and then test the technique in a Monte Carlo framework using three hypothetical populations with different disturbance levels and temporal trajectories. We further demonstrate the technique with a case study.

Materials and Methods
For this research we rely on the statistical design implemented by the USDA Forest Service, Forest Inventory and Analysis (FIA) program. Specifics about the sampling design and estimators are provided by Bechtold and Patterson [17]. Details regarding data collection protocols are provided by the USDA [20].

Plot-Level Disturbance
The FIA program employs a 5 or 7 panel design by state in the eastern United States and a 10-panel design in the western United States. One panel is measured each year, so the remeasurement schedule is approximately equal to the number of panels (P). Up to three observable disturbances (e.g., fire) and the year of each disturbance are recorded. Further, up to three forest management treatments (e.g., timber harvesting) are also recorded, along with the year of each treatment [20]. For simplicity, we refer to both disturbances and treatments simply as disturbances.
Plot-level disturbances are identified in a retrospective fashion, in that a field crew visiting a permanent sample plot determines whether the plot is disturbed relative to the previous visit. Disturbance type and the year of disturbance are the basic data collected, but we must think of the data in a different way for estimation purposes. In fact, under a 5-panel design, there are five different observations collected, one observation for each year since the last measurement. This scenario is expressed in Table 1, which represents a hypothetical example based on a 5-year remeasurement period, Forests 2020, 11, 1298 3 of 14 with "1" indicating a disturbance, "0" indicating no disturbance, and "-" indicating unobservable (i.e., a year in the future from the field visit or a year preceding the previous visit). In this example, plot 1 was visited in 2008, and the data recorded were binary indicators for whether the plot was disturbed for each year since the last measurement. As such, there are five retrospective disturbance observations recorded; one for each year between field visits. Plots from other measurement years, or panels in the five-panel rotation, have the same pattern. In many cases, no disturbances are observed, as illustrated by plot 4 in Table 1, having zeros recorded for each of the possible disturbance years. Table 1. Hypothetical example showing retrospective disturbance data collected for a five-year rotating panel forest inventory. "1" indicates a disturbance was recorded, "0" indicates a disturbance wasn't recorded, and "-" denotes an unobservable.

Post-Stratified Estimator
Bechtold and Patterson [17] (2005) describe the post-stratified estimator used by the FIA program. The technique is post-stratified in that permanent plot locations are used and strata are determined based on ancillary data after the sample was drawn. The ancillary data are typically derived from classified satellite imagery [21] such as the National Land Cover Database (NLCD) land cover map or the NLCD tree canopy cover map [22,23]. The estimate of the population mean for a domain (d) of interest for parameter y is: where, W h is the post-stratification weight, n h is the sample size in stratum h = 1 . . . H, and y i is the observation from the ith plot of variable y. The estimated variance of y d is: The estimate of the total is thenŶ where A T is the total area of the population. The standard error of the estimate is SE( In practice, y hid is either a proportion (e.g., proportion of the ith plot in stratum h and domain d that was forested) or on an areal unit basis (e.g., tree volume m 3 ha −1 ). The sampling frame is area-based and the sample is equal probability, so the sampling weights are A T /n and the estimation weights (assuming W h is known) are A T W h /n h . For estimation, all panels are combined, and the same stratification is applied throughout and, thus, weights are constant across panels and disturbance years for a given stratum.

Sampling Weight Adjustment
The rotating panel design is assumed to produce an equal probability sample of the population for each panel and across panels. However, we contend that the retrospective observation of disturbance  (Table 2). Table 2 shows an example with field measurements (indicated by M) collected from 2006 through 2010 and the year that a disturbance could be reported for each p panel (indicated by x). Note that disturbances that occurred in 2002 can only be observed during the measurement of panel 1 in 2006 and, likewise, disturbances that occurred in 2010 can only be observed in panel 5. However, disturbances that occurred in 2006 can be observed in all panels in the example. Table 2. Schematic depicting the Forest Inventory and Analysis (FIA) panel design where "M" identifies the year each panel (p) is measured and "x" identifies the corresponding years where disturbance can be retrospectively recorded by a field measurement. The total number of plots (n) from panels p = 1 to 5 are decomposed by panel (n p ) and disturbance year (n d ). The overall sampling weight Ω = A T /n is also decomposed into panel sampling weight (Ω p ) and disturbance year sampling weights (Ω d ).
The overall sampling weight Ω = (A T /n) arises from the sample design and is the inverse of the inclusion probability. Given that A T is fixed for any particular population, it is inappropriate to alter the weight without justifying a change in n. However, from Table 2 The effective sample size can be calculated as n d = nD, where n is a 1 × P matrix of the number of field plots in each panel and D is a P × D design matrix, where each row corresponds to the disturbance years that can be observed for each panel. D can take on numerous forms, depending on how many years of retrospective observations are considered. If one assumes that a disturbance and the disturbance year can be collected since the last plot measurement then D = This same approach can be used to determine the sample size for post-stratified estimation of disturbance. In this case n hd = n h D where n h is an H × P matrix of the number of plots by panel for each h stratum and the resulting n hd matrix is H × D: The total sample size for domain d is n d = 1n hd where 1 is a 1 × D matrix of ones. The elements from n hd and n d corresponding to h and d in Equations (2)-(4) are used: The post-stratified estimator as presented in Equations (1)-(4) has not changed; rather, the correct within stratum and across strata sample size is used for the time domain of interest. It is important to note that the original design-based sampling weights (A T /n) and the post-stratified weights (A h /n h ) are incorrect and the correct post-stratified weights are A h /n hd .

Monte Carlo Simulation
We considered two estimation approaches. The PS (post-stratified) approach employed the post-stratified estimator as presented in equations 1 through 4. The WAPS (weight adjusted post-stratified) method employed the weight adjustment as presented in equations 1, and 5 through 7. To compare these approaches, we identified three populations from the forest disturbance time-series maps [11,24], which were considered the true population totals for area of disturbance. These maps identified forest disturbance, causal agent, and disturbance year from 1994-2010 at a 30 m spatial resolution. We used these populations for a Monte Carlo simulation. We selected populations in southwest Arkansas, the northern interior portion of California, and the Aroostook region in Northern Maine (  (Homer et al. 2016) was used to create three strata for each population. For the Arkansas population, the strata were 0-9% canopy cover (h = 1), 10-79% canopy cover (h = 2), and 80-100% canopy cover (h = 3). The strata for the California population were 0-9% canopy cover (h = 1), 10-49% canopy cover (h = 2), and 50-100% canopy cover (h = 3). For the Maine population the strata were 0-9% canopy cover (h = 1), 10-69% canopy cover (h = 2), and 70-100% canopy cover (h = 3). The strata boundaries for each population were identified following the Westfall et al. (2011) recommendations on minimum within strata and total sample sizes. Weights (W h ) were recorded for each stratum in each population. Each population was sampled based on the FIA hexagonal sampling frame (Bechtold and Patterson 2005), where each hexagon was approximately 2400 ha and contained one sample plot. Each hexagon was assigned to 1 panel in a 5 panel design, and, for our purposes, we assumed panel 1 was measured in 2006, panel 2 was measured in 2007, etc., and panel 5 was measured in 2010.
Step 1 was to generate a random plot center location within each hexagon.
Step 2 was to assign the plot to a stratum.
Step 3 was to generate the plot design based on each plot center. The plot design was a four-point cluster with point 1 being subplot 1 and plot center. Three additional subplots were located 36.6 m away from plot center at azimuths of 0 degrees, 120 degrees, and 240 degrees. Each subplot Each population was sampled based on the FIA hexagonal sampling frame (Bechtold and Patterson 2005), where each hexagon was approximately 2400 ha and contained one sample plot. Each hexagon was assigned to 1 panel in a 5 panel design, and, for our purposes, we assumed panel 1 was measured in 2006, panel 2 was measured in 2007, etc., and panel 5 was measured in 2010.
Step 1 was to generate a random plot center location within each hexagon.
Step 2 was to assign the plot to a stratum.
Step 3 was to generate the plot design based on each plot center. The plot design was a four-point cluster with point 1 being subplot 1 and plot center. Three additional subplots were located 36.6 m away from plot center at azimuths of 0 degrees, 120 degrees, and 240 degrees. Each subplot was 0.017 ha in size. In Step 4, the four subplots were then intersected with the populations, and the proportion of area disturbed across subplots by disturbance year was recorded as the observation for the ith plot. In calculating where Y was the true population value for each disturbance year and population. To understand the characteristics of the PS and WAPS approaches we examined the results in the following ways. First, for each approach, we compared the distribution ofŶ from the Monte Carlo replicates to the true population total (Y). Second, we examined the bias of the PS and WAPS approaches. Third, for each approach, we compared the distribution of SE(Ŷ) from the Monte Carlo replicates to the RMSE across Monte Carlo replicates.

Results
There were distinct differences between the PS and WAPS estimate of the total. Overall, the PS method tended to underestimate the true population total where the inter-quartile range (IQR) of PS estimates generally failed to include the true population value (Figure 2). There were two exceptions. Exception 1 was for the 2006 disturbance year. This exception was because the sample was properly weighted for that year. Exception 2 was for the CA population, where estimates in the beginning of the period (2002-2007) exhibited little underestimation. This was because there was little fire disturbance at the beginning of the time-period, and proper weighting of the sample was, hence, less consequential. The underestimation of the PS approach generally increased towards the beginning of the time-period (2002) and the end of the time period (2010) considered. The width of the IQR was smaller for the PS estimates of the total as compared to the WAPS estimates. The WAPS estimates of the total generally followed the true population total where in all cases the true total fell within the IQR of the WAPS estimates. To understand the characteristics of the PS and WAPS approaches we examined the results in the following ways. First, for each approach, we compared the distribution of from the Monte Carlo replicates to the true population total (Y). Second, we examined the bias of the PS and WAPS approaches. Third, for each approach, we compared the distribution of SE( ) from the Monte Carlo replicates to the RMSE across Monte Carlo replicates.

Results
There were distinct differences between the PS and WAPS estimate of the total. Overall, the PS method tended to underestimate the true population total where the inter-quartile range (IQR) of PS estimates generally failed to include the true population value (Figure 2). There were two exceptions. Exception 1 was for the 2006 disturbance year. This exception was because the sample was properly weighted for that year. Exception 2 was for the CA population, where estimates in the beginning of the period (2002-2007) exhibited little underestimation. This was because there was little fire disturbance at the beginning of the time-period, and proper weighting of the sample was, hence, less consequential. The underestimation of the PS approach generally increased towards the beginning of the time-period (2002) and the end of the time period (2010) considered. The width of the IQR was smaller for the PS estimates of the total as compared to the WAPS estimates. The WAPS estimates of the total generally followed the true population total where in all cases the true total fell within the IQR of the WAPS estimates.  Neither the WAPS nor PS method was unbiased (Figure 3). The PS method yielded biased estimates across years and populations, however, estimates for 2002 and 2010 generally showed the largest bias. The exception for the CA population was due to little disturbance in the beginning and end of the time period. Estimates arising from the PS method typically were negatively biased (i.e., underestimates). While the WAPS method produced estimates with less bias than the PS method, the WAPS method did not result in a mean bias of zero for any particular estimate. However, the WAPS method did not tend to either underestimate or overestimate disturbance area. Neither the WAPS nor PS method was unbiased ( Figure 3). The PS method yielded biased estimates across years and populations, however, estimates for 2002 and 2010 generally showed the largest bias. The exception for the CA population was due to little disturbance in the beginning and end of the time period. Estimates arising from the PS method typically were negatively biased (i.e., underestimates). While the WAPS method produced estimates with less bias than the PS method, the WAPS method did not result in a mean bias of zero for any particular estimate. However, the WAPS method did not tend to either underestimate or overestimate disturbance area. The standard errors (SE) of the PS estimates were generally smaller than the SE of the WAPS estimates ( Figure 4). This was because under the WAPS approach the sample size was adjusted downward for each disturbance year, except 2006, when the sample was the same for the two methods. One would generally expect the distribution of the SE of the estimate to be similar to the RMSE from a simulation approach. The IQR of the SE of the WAPS estimate generally contained the RMSE. However, there were large differences between the IQR of the SE of the PS estimate, except for the 2006 disturbance years, and when there was little disturbance in the CA population ( Figure  4). The standard errors (SE) of the PS estimates were generally smaller than the SE of the WAPS estimates ( Figure 4). This was because under the WAPS approach the sample size was adjusted downward for each disturbance year, except 2006, when the sample was the same for the two methods. One would generally expect the distribution of the SE of the estimate to be similar to the RMSE from a simulation approach. The IQR of the SE of the WAPS estimate generally contained the RMSE. However, there were large differences between the IQR of the SE of the PS estimate, except for the 2006 disturbance years, and when there was little disturbance in the CA population ( Figure 4).

Case Study
To demonstrate the WAPS method, we constructed estimates of the extent and timing of disturbances in east Texas as well as tree mortality volume. We selected east Texas as the case study area, owing to the occurrence of several prominent disturbance events that occurred during the period of annual data collection 2001-2018 [25]. East Texas is an area encompassing approximately 9 million hectares, of which about 4.9 million hectares is forest land [26]. The temperate forest biome that covers much of the southeastern US reaches its southwestern extent in east Texas [27,28]. In 2005, Hurricane Rita passed through east Texas with sufficient force to impact 6 percent of the total timber volume in east Texas [29]. In 2008 a second hurricane, Ike, impacted 4 percent of the total timber volume [30]. In 2011, an exceptional drought occurred that resulted in the mortality of approximately 4 percent of all trees in the region [31][32][33].

Case Study
To demonstrate the WAPS method, we constructed estimates of the extent and timing of disturbances in east Texas as well as tree mortality volume. We selected east Texas as the case study area, owing to the occurrence of several prominent disturbance events that occurred during the period of annual data collection 2001-2018 [25]. East Texas is an area encompassing approximately 9 million hectares, of which about 4.9 million hectares is forest land [26]. The temperate forest biome that covers much of the southeastern US reaches its southwestern extent in east Texas [27,28]. In 2005, Hurricane Rita passed through east Texas with sufficient force to impact 6 percent of the total timber volume in east Texas [29]. In 2008 a second hurricane, Ike, impacted 4 percent of the total timber volume [30]. In 2011, an exceptional drought occurred that resulted in the mortality of approximately 4 percent of all trees in the region [31][32][33]. Fire was also prominent in 2011, during the drought event and in subsequent years. The expectation is properly constructed estimates will show spikes of wind disturbance in 2005 and 2008 and drought and fire in 2011.
Based on the design and results of the simulation study one might consider just using the panels needed to construct estimates for the year(s) of interest, however, the execution of the rotating panel design is rarely as simple as described in the simulation. The actual measurement of a panel may span more than one year and, hence, the remeasurement period rarely equals the number of panels. Further, some plots may be measured out of panel because of temporary hazardous conditions and other practical problems that field crews and inventory implementers face.
Trends are similar for both disturbance area and mortality volume ( Figure 5). For consistency with area estimates, volumes reported are mortality (any cause) on the conditions of noted disturbance. Estimates show large spikes of wind disturbance occurring in 2005 and 2008, which coincide with Hurricanes Rita and Ike, respectively. Hurricane Rita is generally believed to have had greater impact on timber resources than Hurricane Ike [29,30], but estimates here suggest similar impacts, at least in terms of area disturbed and mortality volume. The first sign of drought impacts is in 2011, the year in which the worst 1-year drought in Texas history occurred. Drought impacts Based on the design and results of the simulation study one might consider just using the panels needed to construct estimates for the year(s) of interest, however, the execution of the rotating panel design is rarely as simple as described in the simulation. The actual measurement of a panel may span more than one year and, hence, the remeasurement period rarely equals the number of panels. Further, some plots may be measured out of panel because of temporary hazardous conditions and other practical problems that field crews and inventory implementers face.
Trends are similar for both disturbance area and mortality volume ( Figure 5). For consistency with area estimates, volumes reported are mortality (any cause) on the conditions of noted disturbance. Estimates show large spikes of wind disturbance occurring in 2005 and 2008, which coincide with Hurricanes Rita and Ike, respectively. Hurricane Rita is generally believed to have had greater impact on timber resources than Hurricane Ike [29,30], but estimates here suggest similar impacts, at least in terms of area disturbed and mortality volume. The first sign of drought impacts is in 2011, the year in which the worst 1-year drought in Texas history occurred. Drought impacts were more severe than either hurricanes, a finding consistent with another study that looked at the impacts of three events on standing dead trees [25]. The drought affected at least 48% more area and 62% more mortality volume than either hurricane. Substantial drought impacts are estimated for 2012, but much reduced in 2013-2015. Fire is also most prominent in 2011, as expected. Interestingly, fire has continued to affect sizeable area of east Texas 2012-2018, but impacts on mortality volume have not been very large. Fire estimates reported here are combined estimates of three different fire categories: crown and ground fire, ground fire, and crown fire. Crown fires are not that common, with 2011 being a notable exception where approximately 50% of the fires were crown fires. More typical is ground fire, which likely explains the smaller mortality volume even when sizeable areas have been impacted by fire as is the case in recent years (e.g., approximately 93% of fires were ground fires in 2018).
have not been very large. Fire estimates reported here are combined estimates of three different fire categories: crown and ground fire, ground fire, and crown fire. Crown fires are not that common, with 2011 being a notable exception where approximately 50% of the fires were crown fires. More typical is ground fire, which likely explains the smaller mortality volume even when sizeable areas have been impacted by fire as is the case in recent years (e.g., approximately 93% of fires were ground fires in 2018).

Discussion
The goal of this research was to develop an appropriate approach for estimating the annual areal extent of disturbance under a rotating panel sample design when multiple panels are combined for estimation. Our literature review revealed no guidance for this problem, a finding that we contend is because the issue is not with the estimator, but rather the sampling weights. Most of the available literature focuses on estimation assuming that sampling weights are defined during the initial survey design. The design-based weights of the FIA sample are equal within a population; however, considering the nature of retrospective observation, the correct sampling weights are unequal when multiple panels are considered. The change in the sampling weights for estimating annual disturbance area are not immediately clear to the practitioner. Further, standard estimation tools do not take the correct sampling weights into account when constructing estimates of annual disturbance extent. The exception is when the disturbance year of interest is the mid-point (i.e., 2006

Discussion
The goal of this research was to develop an appropriate approach for estimating the annual areal extent of disturbance under a rotating panel sample design when multiple panels are combined for estimation. Our literature review revealed no guidance for this problem, a finding that we contend is because the issue is not with the estimator, but rather the sampling weights. Most of the available literature focuses on estimation assuming that sampling weights are defined during the initial survey design. The design-based weights of the FIA sample are equal within a population; however, considering the nature of retrospective observation, the correct sampling weights are unequal when multiple panels are considered. The change in the sampling weights for estimating annual disturbance area are not immediately clear to the practitioner. Further, standard estimation tools do not take the correct sampling weights into account when constructing estimates of annual disturbance extent. The exception is when the disturbance year of interest is the mid-point (i.e., 2006 in the example). We have demonstrated that the post-stratified estimator is suitable to estimate annual disturbance totals, but the effective sample size and, hence, sampling weights, must be calculated. It turns out to be a rather simple task to correctly calculate the effective sample size and sampling weights for annual disturbance estimates. These results are relevant to other NFIs that use rotating panel design and collect disturbance information in a retrospective fashion.
While we considered the WAPS methods suitable to estimate annual disturbance totals, it is important to note that the estimates were biased based on the simulation results. The percent bias of the WAPS method was 5.8%, across years and populations. The bias arose because the disturbances were rare events or small domains and, under traditional NFI sampling intensities, only a small number of disturbed plots may be observed in a given year. We note that bias was typically minimized for 2006 estimates (Figure 3). This was because the effective sample size included all plots from all panels ( Table 2), which further suggests that increased sample sizes may be required to achieve unbiased estimates of disturbance extent by year. As suggested by Kish [15] and Cochran [14], there are sampling techniques that are appropriate for rare events. However, most NFIs were not designed as rare event surveys and, further, their sample designs and intensities are relatively fixed at the national level. Alternative estimators such as, model-assisted, model based, and small area techniques may prove valuable for reducing the bias and increasing precision of annual disturbance estimates.
The analysis presented examined the performance of the WAPS approach for estimating disturbance extent and timing. The reason for this emphasis was because information was available to serve as a true population in the Monte Carlo simulation. However, our results are applicable to a suite of domain estimates that can be constructed given the variables collected as part of broad-scale forest inventories. For example, one may use the WAPS technique to estimate the volume of trees that died due to fires in a given year, because tree mortality year and cause are also collected, retrospectively. Likewise, annual estimates of removal volumes can also be constructed. There are numerous other estimates relating to disturbance impacts that can be constructed using the proposed technique.
There are other methods to account for the retrospective observations when all panels are combined for estimation. One straightforward approach is to annualize the observations. Observations can be annualized by simply dividing the observed value from the plot by the length of the measurement interval for the plot. This technique is currently used by FIA when constructing estimates of growth, removals, and mortality [34]. If the retrospective observations of disturbance are annualized then no adjustment to the sampling weights is needed. However, the average annual estimate of the disturbed area for any population would cover an approximately 2P year span. Actual annual estimates are more informative for addressing current resource issues than an estimate of average annual disturbance area over a 10+ year time span. Further, as Roesch [35] points out, these average annual estimates, calculated over different measurement intervals, are different population parameters.
Our focus here was on developing an estimation technique for retrospective observations under a rotating panel design, and we did not address measurement error issues related to plot-level retrospective observation. Most NFIs have published quality objectives, and those of the FIA program suggest that the year of disturbance should be within one-year of the actual disturbance year 99% of the time [20]. Pollard and colleagues [36] found that this quality objective is generally met. However, the study by Pollard and colleagues [36] should be updated given the expanding FIA data volume available. In some cases, a practitioner may wish to only consider disturbances (and year of disturbance) recorded within the three-years of plot measurement (for example). In this case, it is relatively straightforward to modify the D matrix to take into account that only three years of retrospective observations are being used to construct the estimates.

Conclusions
Disturbance at varying frequencies, intensities, and spatial scales is an important component of forest systems. Understanding the extent of disturbances and their impacts is critical to formulating effective management and policy responses. A key information source is large-area forest inventory data, such as those produced across the US by the FIA program. Due to the retrospective nature of the observations in a panelized design setting, the use of standard estimators that generate typical forest resource statistics can lead to inaccurate representations of the areal extent and timing of disturbance. In this paper, the post-stratified estimator was extended to account for the specific manner in which disturbance data are recorded. In addition to the theoretical development, the estimates exhibited satisfactory agreement with the true population values when assessed in a simulation environment. Broadly, the increasing use of panelized forest inventory designs globally suggests it will often be critical for researchers to account for the necessary weighting adjustments by observation year of the panel when estimates of areal extent and timing of disturbances are desired. Thus, we recommend the methods presented in this paper be adopted or the concepts generally incorporated into other types of forest inventory estimation frameworks.