Comparison of Model-Assisted Endogenous Poststratification Methods for Estimation of Above-Ground Biomass Change in Oregon, USA

Abstract

Quantifying above-ground biomass changes, ΔAGB, is key for understanding carbon dynamics. National Forest Inventories, NFIs, aims at providing precise estimates of ΔAGB relying on model-assisted estimators that incorporate auxiliary information to reduce uncertainty. Poststratification estimators, PS, are commonly used for this task. Recently proposed endogenous poststratification, EPS, methods have the potential to improve the precision of PS estimates of ΔAGB. Using the state of Oregon, USA, as a testing area, we developed a formal comparison between three EPS methods, traditional PS estimators used in the region, and the Horvitz-Thompson, HT, estimator. Results showed that gains in performance with respect to the HT estimator were 9.71% to 19.22% larger for EPS than for PS. Furthermore, EPS methods easily accommodated a large number of auxiliary variables, and the inclusion of independent predictions of ΔAGB as an additional auxiliary variable resulted in further gains in performance.

1. Introduction

National forest inventories, NFIs, collect measurements and provide the necessary infrastructure to estimate the state and dynamics of the forested areas of their respective countries [1]. To satisfy the demands of information posed by society, national legislation, or international agreements, NFIs have to allow the estimation of means and totals of a large and increasing number of variables of a very different nature [1,2,3]. While NFI sampling designs can vary substantially between countries (see [3]), several features are common to most NFIs. Sampling designs of NFIs typically use systematic or near-systematic approaches [4] based on regular grids, and measurements of field plots are organized over time in a systematic and cyclic fashion. This ensures that samples are spatially balanced [5] and at least a fraction of the field plots are revisited at a frequency that enables NFIs to estimate changes in ecologically important attributes [1,3]. The resulting sampling designs are generalist and standardized within countries; and making possible the estimation of the state and dynamics for multiple forest attributes and facilitating comparisons between spatial domains or time periods is typically considered a priority.

The large degree of standardization and the generality that allow NFI sampling designs to estimate many attributes also implies that they are not optimized in terms of precision for any particular variable [3]. Therefore, the uncertainty of estimates for specific variables can be large. This problem is typically addressed using estimators that incorporate auxiliary information to gain precision in the estimation of key variables. While model-based estimators [6] could be used for this purpose, model-assisted estimators such as generalized regression estimators, GREG, post-stratification estimators, PS, or calibration estimators are typically the preferred choice. This is because model-assisted estimators: (1) allow using abundant sources of auxiliary information to improve the precision in the estimation process [7] and (2) do not involve strong assumptions about the structure of the population of interest in the shape of regression models being approximately unbiased even if the assisting models are not correctly specified [7].

In the NFI of the United States of America, the US Forest Service Forest Inventory and Analysis program, FIA, PS is typically used to estimate both state and dynamic variables. In the Pacific Northwest FIA region, PNW-FIA (i.e., states of Alaska, Hawaii, Oregon, Washington, and California, and the Pacific Island territories), the PS currently used for all variables was developed using expert opinion based on ownership, ownership size, canopy cover, and a classification into forest, non-forest and unknown. PNW-FIA has used the resulting PS during the last 7 years, and it has proven to satisfy the reporting standards of FIA for state variables such as above-ground tree biomass, AGB, or volume, even when disaggregating by broad size or species classes or by large geographical domains such as states or large counties. The auxiliary variables used to develop this PS are primarily static descriptors, and when this PS is used to estimate the change in forest AGB, $\mathrm{\Delta }AGB$, the precision of the PS estimators declines substantially. A potential solution to this lack of precision is developing new model-assisted estimators with dynamic auxiliary variables specially tailored to improve the precision for change variables. However, several features must be analyzed before a model-assisted estimator and a set of auxiliary variables are incorporated into the FIA protocols.

Model-assisted techniques yield improvements in precision that are directly related to the correlation between auxiliary variables and target response. Thus, a cornerstone for developing a new estimator for $\mathrm{\Delta }AGB$ is the identification of auxiliary variables potentially correlated to this variable. Positive $\mathrm{\Delta }AGB$ (i.e., an increase) can be expected to be the result of steady growth dynamics. In contrast, negative $\mathrm{\Delta }AGB$ can result from steady decline dynamics (i.e., senescence or slow diseases) or from events of varying intensity, such as insect outbreaks, fires, or harvest operations such as thinning or clear-cuts. Thus, rather than using a single auxiliary variable, a model-assisted estimator for $\mathrm{\Delta }AGB$ will require multiple and complementary auxiliary variables capturing information about the different drivers of change.

Numerous sources of auxiliary variables potentially correlated to factors affecting $\mathrm{\Delta }AGB$ are currently available. However, the potential correlation with $\mathrm{\Delta }AGB$ cannot be the only factor under consideration. Before a source of auxiliary information is considered appropriate for developing new model-assisted estimators, it must meet some standards in terms of spatial coverage, temporal extent, and update frequency. Otherwise, compatibility with the 10-year rotating panel design of the PNW-FIA region and harmonization with current reporting practices could become infeasible. These considerations are especially relevant for recent remote sensing products that: use Landsat time series and topographic and climate variables to provide yearly predictions of AGB and that are developed using ground observations that are independent of the FIA sampling design [8,9]. By differencing AGB predictions from specific years, it is possible to derive predictions of ΔAGB that are highly correlated with actual ΔAGB. These predictions of ΔAGB can be used as auxiliary variables that encode, in a single auxiliary variable, the most relevant information provided by the remote sensing datasets used in their construction. These datasets are very recent, but their potential correlation with $\mathrm{\Delta }AGB$ and the fact that they were derived with ground data independent of FIA, call for investigations about their use in FIA estimation.

Some model-assisted estimators have been developed and could be used to improve the precision of estimates of $\mathrm{\Delta }AGB$. However, we will focus the remainder of this manuscript on PS estimators. These estimators can be seen as a sub-type of GREG estimators where the auxiliary variables are stratum-membership indicator variables [7]. From an operational perspective, PS estimators are the most appealing alternative because they imply a continuation of the methods currently used in the region, which is an important advantage of PS over more general formulations of GREG estimators. In addition, PS estimators can account for differences in sampling intensities between strata that occur after denial of access to plot locations or hazardous plot conditions. Additionally, PS estimators have a large degree of flexibility and ensure compatibility (i.e., positive weights are ensured) with the design of the FIA database and associated software [10].

Several alternatives exist to develop a PS estimator. Traditional PS methods in PNW-FIA have used expert knowledge and geographic information systems, GIS, and data about forest types, ownership, administrative boundaries, etc., to define homogeneous units or post-strata. The key aspect of this traditional approach to PS is that it is developed “exogenously” or independently of the sampled data [11]. Recent advances in survey sampling have proposed new methods to develop PS in a sample data-driven fashion [11,12,13,14]. Estimators resulting from these approaches are called endogenous PS, EPS estimators to denote that post-strata are defined considering the information contained in the sample (i.e., “endogenously”). To differentiate endogenous and exogenous methods, we will use PS for the former and EPS for the latter in the remainder of the manuscript. EPS estimators have the potential to (1) yield larger improvements in performance than PS because they fully use the sample information and (2) reduce analyst involvement and subjectivity because automatic model selection methods can be used in the development of an EPS. However, an evaluation of EPS methods for the estimation of $\mathrm{\Delta }AGB$ is needed.

Objective

The main objective of this study was to test different EPS estimators for forest $\mathrm{\Delta }AGB$. To achieve this objective, we compiled a spatial database with auxiliary variables that were considered appropriate for their incorporation into FIA protocols. The database included variables currently used by FIA and also new variables that are expected to reflect changes in forests $AGB$. Then, the following EPS estimators were subject to analysis:

• A GREG estimator, GREG-EPS, based on the interactions of multiple categorical variables [13].
• An EPS estimator using the method proposed by Breidt and Opsomer [11] that uses a generalized linear model to form strata, GL-EPS.
• An EPS estimator using McConville and Toth’s [14] method based on recursive partitioning trees, TREE-EPS.

These estimators were selected because published results proving their consistency were available and because they allow using a large number of auxiliary variables, which, a priori, makes them suitable candidates to improve over traditional PS methods. We focused our analyses on the estimation of ΔAGB for 10-year periods, as that is the re-measurement cycle in the PNW-FIA rotating panel. For each EPS estimator, we estimated the gains in precision that they provided with respect to both: the Horvitz-Thompson estimator, HT, and the PS estimators currently used by FIA in the region, FIA-PS, and analyzed the interpretability, ease of use, and potential for incorporation into the FIA estimation protocols of each EPS variant. As secondary objectives, we (1) examined the improvements in performance provided by each estimator when the target was to estimate running means of ΔAGB for consecutive 10-year periods, and (2) assessed the possibility of further improving the performance of these EPS estimators by incorporating as additional auxiliary variables, predictions of ΔAGB derived in a recent Carbon Monitoring System, CMS, project [8,9]. These predictions were derived using remote sensing data and ground observations impendent of the FIA sampling design and will be referred to as ΔCMS.

2. Material and Methods

2.1. Study Area

This study was developed using the US state of Oregon as a testing area (Figure 1). Oregon is included in the PNW-FIA region, and its size, geographic position, and diversity make it well-suited for testing purposes. Oregon is 247,963 Km² in size, and 48.4% is forested land, of which 64% is publicly owned [15]. Oregon has very strong west-to-east environmental gradients. The area west of the Cascade Range is moist and mild, with forest dominated by Douglas-fir, Pseudostuga menziesii (Mirbel) Franco, red alder, Alnus rubra Bong, western hemlock and Tsuga heterophylla (Raf.) Sarg. Areas east of the Cascade Range have a dry, continental climate, and the forests are dominated by ponderosa pine, Pinus ponderosa Lawson & C. Lawson, lodgepole pine, Pinus contorta Douglas ex Loudon, and western juniper, Juniperus occidentalis Hook.

2.2. Auxiliary Information

After the initial literature review, we compiled a set of auxiliary variables available for the PNW-FIA region that were potentially correlated with $\mathrm{\Delta }AGB$. This database contained proxies for potential forest AGB productivity, disturbances and AGB removals, and forest management regimes. All auxiliary variables were rasterized using a common 30 × 30 m grid using the NAD83 2011 reference system with conic Albers equal area projection to ensure that pixel counts can be directly translated to areas. Auxiliary variables contained variables considered static (i.e., did not change over time) and variables considered dynamic. For example, all proxies for potential AGB productivity, such as climate or topographic position, are not expected to change substantially within the period considered in this study and, therefore, were considered static. Proxies for management regimes were based on ownership maps. Ownership can change over time; however, in our analyses, we only considered extensive ownership categories that are not expected to have major changes within the study period. Thus, ownership was also considered a static variable. Proxies for disturbance and AGB removals and the annual CMS AGB maps were all dynamic and their preprocessing focused on obtaining variables that referred to 10-year periods matching the PNW-FIA panel design. For auxiliary variables computed as changes in land cover maps and changes in the annual CMS AGB map, it was impossible to calculate values for some of the 10-year periods considered in the analysis. In those cases, we allowed small temporal mismatches of up to two years with respect to the 10-year cycle under consideration. Finally, some auxiliary variables were categorical, and some were continuous. In Supplementary S1 of the supplementary material, we describe in detail: (1) the processing steps necessary to include each dataset in the database, (2) the derived auxiliary variables, (3) their categorical/continuous nature, and (4) whether they were considered as static or dynamic variables. This information is summarized in

Table 1. Auxiliary information database summary. Categorical variables (Ca), Continuous variables (Co). Categorical variables derived from continuous ones are marked with an asterisk (i.e., *Ca).

Type	Source	Variables, Acronym	Pre-Processing	Variable Type	Temporal
Proxies for potential forest AGB productivity	1 Arc second Shuttle Radar Topography Mission, SRTM, Google earth engine	Elevation, $ELEV$	Resampling bilinear interpolation	Co	Static
		Elevation categories, $CATELEV$	Division in 3 elevation categories of equal area	*Ca
		Slope, $SLP$	Computed from $ELEV$	Co
		Heat load index, $HTL$	Computed from $ELEV$	Co
	800 m resolution PRISM 30-year normals & Sun hours from SRTM	Paterson climate productivity index, $PCPI$	Resampling bilinear interpolation.Solar radiation ArcGIS tool	Co
		Categories Paterson climate productivity index, $CATPCPI$	Division in 3 categories of equal area	*Ca
	US Forest Service	Cleland’s level 3 ecoregions, $ECO$	Rasterization	Ca
Ownership	Bureau of land management, BLM	Ownership, $OWN$	Rasterization & reclassification	Ca
Proxies for disturbance	Monitoring trends in burn severity, MTBS.	Fire severity, $MAXFSE{V}_t$	Maximum fire severity for 10-year periods. Resampling nearest neighbors	Ca	Dynamic
	Landscape Change Monitoring System, LCMS.	Disturbances, $ACDIS{T}_t$	Computation of accumulated disturbances for 10-year periods	Co
	MTBS- LCMS	Disturbance-categories, $CATACDIS{T}_t$	Reclassification of $MAXFSE{V}_t$ and thresholds for $ACDIS{T}_t$	*Ca
	MRLC National Land Cover Database, NLCD.	Land cover change, $\mathrm{\Delta }NLC{D}_t$	Resampling nearest-neighbor. Reclassification and computation of change	Ca
Change in multi-year CMS AGB map	Fekety and Hudak, (2019) & Hudak et al., (2020) [8,9]	Independent prediction of $\mathrm{\Delta }AGB$ derived from Fekety and Hudak, (2019) [8] predictions of AGB for multiple years, $\mathrm{\Delta }CM{S}_t$	Resampling with bilinear interpolation. Computation of differences in predicted forest AGB between years.	Co
		Categories of $\mathrm{\Delta }AGB$ change derived from independent predictions of AGB for multiple years, $CAT\mathrm{\Delta }CM{S}_t$	Reclassification of $\mathrm{\Delta }CM{S}_t$ based on intervals defined from values reported by [16]	*Ca

.

2.3. Estimation Framework

2.3.1. Population

We considered a finite population approach for estimation purposes. For the 10-year period starting year $t$, the set of all 30 m × 30 m spatial pixels in the state of Oregon was considered as a finite population of interest, $U$. The base grid of the auxiliary information database is itself the sampling frame and does not change over time. Finally, pixels are considered population units and are indexed by $i=1,2,\dots N$, where $N$ is the total number of pixels for the state. For the 10-year period starting in year $t$, every population unit $i$ (i.e., pixel) has an associated q-dimensional vector of auxiliary variables denoted as ${\mathit{x}}_{\mathit{i}\mathit{t}}={\left(x_{i1t,}{x}_{i2t},\dots ,x_{iqt}\right)}^{\mathit{t}}$. This vector contains all auxiliary variables for the 10-year period starting in year $t$. Every unit $i$ also has associated values $\mathrm{\Delta }AGB$ for the period, but those changes are unknown except for the locations that were measured, according to the panel design, at the beginning and end of the period (see Section 2.3.3).

2.3.2. Target Parameter

PNW-FIA considers changes corresponding to 25 possible combinations of conditions at times $t$ and $t+10$. Still, some of these combinations correspond to unsampled areas at time one, time two, or both times or to sampled areas outside of the reporting scope of PNW-FIA. We focused our analysis on changes in areas that were forest at times $t$ and $t+10$. The response for a population unit $i$ in the decade starting in year $t$ is denoted in the remaining sections as $\mathrm{\Delta }AG{B}_{it}$ and is the annualized forest change of AGB per unit area.

$$\mathrm{\Delta }AG{B}_{it}=\frac{1}{10}(AG{B}_{it+10}-AG{B}_{it})$$

(1)

$AG{B}_{it}$ and $AG{B}_{it+10}$ represent the AGB for areas that are forest at years $t$ and $t+10$ in unit $i$ divided by the unit’s area. For non-forested units $AG{B}_{it}$, $AG{B}_{it+10}$ and $\mathrm{\Delta }AG{B}_{it}$ equal zero.

We defined the target parameter for the 10-year period starting on year $t$ as the mean of $\mathrm{\Delta }AG{B}_{it}$

$${\overline{\mathrm{\Delta }AGB}}_t=\frac{1}{10N}\left({\displaystyle \sum }_{i=1}^{N}AG{B}_{it+10}-{\displaystyle \sum }_{i=1}^{N}AG{B}_{it}\right)=\frac{1}{N}{\displaystyle \sum }_{i=1}^N\mathrm{\Delta }AG{B}_{it}.$$

(2)

We expressed the target parameter as a mean value per year and unit area change; however, converting this value to totals only requires multiplying ${\overline{\mathrm{\Delta }AGB}}_t$ by the state’s known area and by 10 years.

We also considered estimating running means covering the periods $t_0$ to $t_0+10$ through $t_0+a$ to $t_0+a+10$. The target parameter in the case of a running mean covering the periods $t_0$ to $t_0+10$ to $t_0+a$ to $t_0+a+10$ is

$${\overline{\mathrm{\Delta }AGB}}_{t0,t0+a}=\frac{1}{a}{\displaystyle \sum }_{t=t_0}^{t=t_0+a}\frac{1}{N}{\displaystyle \sum }_{i=1}^N\mathrm{\Delta }AG{B}_{it}$$

(3)

2.3.3. PNW-FIA Sampling Design and Sample

The sampling design used by PNW-FIA combines the general 2402.62 ha [17,18] FIA hexagonal grid with a 10-year rotating panel [15]. The result is a near systematic sampling design where one-tenth of the plots covering the entire state (i.e., a panel) are measured every year. In the PNW-FIA region, the sampling design was intensified on all National Forest Lands that were not designated Wilderness areas. Additional locations were selected for areas subject to intensification, which resulted in a field plot density 3.2 times higher than in non-intensified areas subject to the standard FIA design. Each of these new locations was assigned to one panel. The assignment to panels was done using a method that preserved the spatial balance of the sample within the intensified areas and provided an approximately constant number of intensification locations per panel [19].

All locations were initially inspected using orthophotos and GIS layers, and only those locations for which there was very strong evidence that they did not sustain any forest, i.e., urban areas, open farmland, or desertic areas, were labeled as non-forest [18]. All the remaining locations (including disturbed areas that currently had no live trees) were visited and measured by field crews that established plots and monumented the plot center to facilitate subsequent re-measurements (see [18] for a description of FIA measurement protocols).

For stratification purposes, three areas called estimation units, EU, are defined for each state in the PNW region. These areas are National Forest lands not designated as wilderness areas, NF, National Forest lands designated as wilderness areas, WL, and other lands, OL (Figure 1). Stratifications in the PNW region are developed in a way that ensures that a given stratum can only be present in a single EU. A layer with the delineation of the EUs was included in the auxiliary information database, and all EPS were developed respecting the constraint that no stratum can be present in more than one EU.

We considered that field plots for which less than 75% of the plot area was measured in both the first and second rounds of measurements were missing data. In addition, we also excluded plots for which the second measurement year was not 10 years after the first one (

Table 2. The total number of FIA plots used in model fitting was reported by time period and Estimation Units (EU). The number of excluded plots, which were classified as water in all NLCD maps, are also reported. Estimation units were: National Forest lands not designated as wilderness areas (NF), National Forest lands designated as wilderness areas (WL), and other lands (OL).

Period	Total Number of Plots	Number of Plots by EU			Excluded Plots
		NF	OL	WL
2001–2011	1310	675	592	28	15
2002–2012	1412	681	682	29	20
2003–2013	1402	687	656	29	30
2004–2014	1418	703	671	29	15
2005–2015	1420	704	662	33	21
2006–2016	1348	680	623	22	23
2007–2017	1331	650	645	18	18
2008–2018	1340	674	616	25	25

). Plots in which the first and second measurements were respectively taken in years $t$ and $t+10$ form the sample, $s_t$, for that 10-year period. Note that extreme cases include plots for which the first and second measurements can be close to 9 or 11 years apart (e.g., a lag between measurements close to 9 years occurs when the first measurement is taken at the end of year $t$ the second measurement at the beginning of year $t+10$). Moreover, note that, as indicated in Section 2.3.1, for many plots, the value of the response is zero because they are in non-forested areas. The sample for estimation of running means, Equation (3), is the union of all the samples for the periods considered in the running mean, i.e., $s_{t_0{t}_0+a}={{\displaystyle \mathrm{U}}}_{t=t_0}^{t=t_0+a}{s}_t$. Finally, we used the criteria in [4,20] and assigned each field plot measurement to the pixel that contained the center of the plot.

2.4. Development of Models for EPS

Estimators under consideration used assisting models for which the expected value of $\mathrm{\Delta }AG{B}_{it}$ equals

$$E\left(\mathrm{\Delta }AG{B}_{it}\right)={\beta }_t{}^t f\left(x_{it}\right)={\beta }_t{}^t z_{it},$$

(4)

where: $x_{it}$, is the $q$ dimensional vector that contains the auxiliary variables in the spatial database for unit $i$, $z_{it}$ results from applying to $x_{it}$ a function $f$ that maps vectors from the auxiliary variables space to a space that contains $p$-dimensional indicator vectors $z_{it}$ providing the stratum membership for each unit. All entries in $z_{it}$ are zero except for the one corresponding to the stratum where unit $i$ is assigned to $z_{itk}$ (i.e., the $k^{th}$ entry of $z_{it}$, equals 1 if unit $i$ belongs to the $k^{th}$ stratum during the period starting year $t$ and zero otherwise). The dimension $p$ provides the number of strata and varies depending on the EPS estimator. The term ${\beta }_t$ is a vector of coefficients, and super index $t$ is the transpose operator. Both $f$ and ${\beta }_t$ need to be determined. The function $f$ is linked to the type of EPS, and model coefficients ${\beta }_t$ need to be estimated using the sample data. Finally, for the assisting models, the variance of $\mathrm{\Delta }AG{B}_{it}$ is a constant for each stratum $k$

$$V\left(\mathrm{\Delta }AG{B}_{it}\right)={\sigma }_k^2, if i \in stratum k$$

(5)

2.4.1. General Model Selection Considerations

The number of variables in the auxiliary information database is large, and model selection steps were necessary to identify the function $f$ used in each model-assisted estimator. All model selection steps were carried out to ensure that the final strata did not cut across EUs. Plots in areas classified as open-water, NLCD code 11, in all the NLCD land-cover maps were excluded from the model selection stage (Table 2). Model selection steps differed between estimators, and specific details for each estimator are provided in Section 2.4.2, Section 2.4.3, Section 2.4.4. Common steps related to the model selection processes were:

First, all model selection steps were done with the datasets that combined the eight 10-year periods under consideration. This allowed identifying unique models for all 10-year periods. All 10-year periods were represented in approximately equal proportions in the combined datasets used for model selection, ensuring that no 10-year interval received more weight than the others.

Second, the model selection steps were only used to identify $f$ and the variables included in each assisting model. For a given 10-year period, final model coefficients were estimated using the field sample corresponding to the panel for the period under consideration, and for multiple 10-year periods, coefficients were estimated using the field sample of the corresponding panels.

2.4.2. GREG-EPS

In GREG-EPS, $f$ uses dummy variables resulting from the interaction of $ECO$, $OWN$, $CATPCPI$, $CATACDIS{T}_t$, $MAXFSE{V}_t$, and $\mathrm{\Delta }NLC{D}_t$. There were 4700 combinations of values of these variables, and most combinations had such a small areal representation that they could not be directly considered potential strata. These interactions were collapsed into groups with an average area for which the expected number of plots per year was four or more. For each EU, the original combinations of categorical variables were collapsed using an algorithm that:

• First, arranged all combinations on a tree that was constructed using the variables $CATACDIS{T}_t$, $ECO$ and $OWN$. Branches in the first level of the tree were defined using $CATACDIS{T}_t$, and branches for the second and third levels were based on $ECO$ and $OWN,$ respectively. Each combination of $ECO$,$OWN$,$CATPCPI$, $CATACDIS{T}_t$, $CATPCPI$, $MAXFSE{V}_t$, and $\mathrm{\Delta }NLC{D}_t$, resulted in a leaf that was placed in its corresponding branch depending on $CATACDIS{T}_t$, $ECO$ and $OWN$.
• Leaves with four plots per year or more were set as fixed leaves. The remaining leaves were merged with other leaves (fixed or not) in the same branch. The merging process ran separately in each branch and started with the leaf with a smaller area in the branch. The selected leaf was merged with the closest leaf in the same branch. The Gower distance [21], computed with all categorical variables, was used to determine which leaf was the closest to the selected leaf. This distance was selected because it allows treating differently categorical variables with an implicit ordering (i.e., all categorical variables derived from a continuous one) and categorical variables without an implicit ordering (e.g., $OWN)$. The two leaves were merged into a single leaf, and the expected number of plots per year was recomputed based on the area of the group resulting from the merge. If the expected number of plots per year of the resulting leaf was four or more, or if the small leaf was merged to a fixed leaf, the result was tagged as fixed, and it was not considered as a target for further merging steps.
• Step 2 was repeated until all the resulting leaves in a branch had an expected number of four plots per year or all leaves in one branch were merged into a single leaf.
• When merging all leaves in one branch did not yield an area from which to expect four plots per year, the merging process continued but considered merging groups from branches of the previous level of the tree (i.e., the algorithm continued with branches defined by $CATACDIS{T}_t$ and $ECO$ first, and then with branches defined by $CATACDIS{T}_t$ only).

This collapsing procedure is analogous to the one used by FIA for PS. Its result was a set of 223 groups, and each of those groups was a potential post-stratum. The number of potential post-strata was still large, and we used the elastic net implementation from the R package glmnet [22] to determine which ones could be considered as final strata. Groups with non-zero coefficients were considered as final strata, and groups with zero coefficients were assigned to selected strata by first tagging the former as fixed and then applying steps 1 to 4 to the groups with zero coefficients but setting an area threshold larger than the state’s area to ensure that groups with zero coefficients were merged to groups with non-zero coefficients. The resulting groups formed the final set of post-strata, and each pixel of the auxiliary information database was assigned to its corresponding group. This resulted in the final PS maps for GREG-EPS.

It is important to note that the combinations of auxiliary variables that give rise to strata are defined using dynamic variables that change over time. To compute the expected sample size per year, we used the average area of each combination of categorical variables over the 8 periods considered in this analysis. This implies that the possibility of having unsampled strata for a particular 10-year period is significantly reduced but not eliminated. After testing with different area thresholds, an average expected sample size of four plots per year was selected, as it maintained the area of the unsampled micro-strata below 2.5% of the total area.

To assess the improvement in the estimation of $\mathrm{\Delta }AGB$ when the auxiliary variables included the multiyear remote sensing CMS AGB map from [8], we repeated this process, adding $CAT\mathrm{\Delta }CM{S}_t$ to the pool of auxiliary variables. The branches of the tree for the collapsing process were defined with the same variables ($CATACDIS{T}_t$, $ECO$ and $OWN$). The elastic net model selection, final collapsing and strata mapping steps were repeated, and the result was an alternative EPS. We will refer to this EPS as GREG-EPS-CMS.

2.4.3. GL-EPS

The GL-EPS estimator proposed by Breidt and Opsomer [11] is developed based on a response that is selected as a PS variable. In our case, the PS variable is $\mathrm{\Delta }AG{B}_{it}$ itself. A set of intervals is defined a-priori by the modeler on the range of the PS variable, and each one of those intervals will result in a final post-stratum. For example, an over-simplified set of intervals for $\mathrm{\Delta }AG{B}_{it}$ could be $\left(-\infty ,\left.0\right]\right.$ and $\left(0,\infty \right)$, i.e., losses and gains of forest AGB. Once the intervals for the PS variable are defined, GL-EPS estimators are obtained assuming that the PS variable can be predicted using a function, $h$, of a linear combination, ${\mathsf{\gamma }}^t{x}_{ti}$, of the auxiliary information available

$${\widehat{\mathrm{\Delta }AGB}}_{it}=h\left({\widehat{\gamma }}^t{x}_{it}\right)$$

(6)

The function $h$ only needs to be a strictly monotone function [11] to ensure consistency of the GL-EPS estimator, which allows for very flexible specifications of $h$. In our study, we considered the simplest case and $h$ was the identity function. To form the linear predictor, we considered all the categorical variables that were not derived from a continuous one (e.g., $CATPCPI$ was not considered because it was derived from the continuous variable $PCPI$), all continuous variables, and all their pairwise interactions. The total number of variables and pairwise interactions present in the sample was 279. A model selection step was performed using the elastic-net regularization criteria [4] to reduce the number of coefficients in the linear predictor.

The intervals for $\mathrm{\Delta }AG{B}_{it}$ were defined using the values reported for this variable by [16] in Panther Creek, western Oregon. This study provides a good reference for large values of AGB growth because the study site is located in one of the most productive areas of the state. Thresholds were defined using the following criteria:

• The maximum value of $\mathrm{\Delta }AG{B}_{it}$ reported by [16] was 8.75 Mg ha⁻¹ year⁻¹. Based on this value, we defined the following five positive intervals $\left(0,\left.2\right],\right.$ $\left(2,\left.4\right],\right.\left(4,\left.6\right]\right.$ $\left(6,\left.8\right],\right.$ and $\left(8,\infty \right)$.
• For disturbances causing losses in forest AGB with magnitudes comparable to growth, we used the thresholds used for growth but with negative signs. The last interval $(-\infty ,\left.-8\right]$ accommodates large and negative values of $\mathrm{\Delta }AG{B}_{it}$ occurring after stand-replacing disturbances such as clear cuts.

Once models were selected, the non-zero coefficients ${\widehat{\gamma }}_t$ were used to obtain predictions $\widehat{\mathrm{\Delta }AG{B}_{it}}$ for each population unit. These predictions were reclassified according to the defined threshold to obtain the corresponding stratum membership vectors $z_{it}$ associated with each population unit and 10-year period. The model selection and mapping of strata was repeated, including $\mathrm{\Delta }CM{S}_t$ in the pool of auxiliary variables. We will refer to the EPS using $\mathrm{\Delta }CM{S}_t$ as GL-EPS-CMS.

2.4.4. TREE-EPS

For TREE-EPS, we applied the method proposed by McConville and Toth [14]. One partitioning tree was identified in the model selection stage for each EU. These trees were found using the rpms R package [23] that implements the partitioning algorithm described in [24]. All auxiliary variables, except $\mathrm{\Delta }CM{S}_t$ and those categorical variables derived from a continuous one, were used to find partitioning trees for each EU. Trees were used for each 10-year period to assign pixels and field plots to strata. This process was repeated including $\mathrm{\Delta }CM{S}_t$ in the pool of auxiliary variables to assess the improvement in the estimation of forest AGB derived from the use of this auxiliary variable. We will refer to the EPS based on tree models using $\mathrm{\Delta }CM{S}_t$ as TREE-EPS-CMS.

2.5. Estimators of $\Delta AGB$ and Variance Estimators

The three EPS variants were compared to the exogenous PS estimators currently used by PNW-FIA, FIA-PS, hereafter, and to the well-known HT estimator. As target parameters, we considered means and running means of forest AGB. Estimators are described in detail for the base case where the target parameter is the mean $\mathrm{\Delta }AGB$ for a single 10-year period. For running means, we present the modifications with respect to the base case that is needed. For both single 10-year periods or running means of more than one 10-year period, the possibility of having unsampled strata exists. This can cause differences in the areas used to make comparisons between methods. In order to overcome this problem, the presence of unsampled strata was addressed by merging all unsampled strata in a given EU to the sampled strata with the largest area in the EU.

2.5.1. Approximation to Sampling Design Weights and Point and Variance Estimators

All estimators under consideration are design-based estimators. Their expressions and associated variance estimators require knowing first and second-order inclusion probabilities. Due to many factors, such as denial of access or hazardous plot locations, the number of plots measured in the field does not match that planned in the design. In order to account for this fact, inclusion probabilities are empirically adjusted based on the sample size and the total number of pixels in the stratum of the plot. This adjustment implies the assumption that plots are missing at random within strata. For plot $i$ measured in the period $t$ to $t$ + 10 in stratum $k$, the adjusted inclusion probability used in the further analysis was

$${\pi }_{it}=\frac{n_{t k}}{N_k} ,$$

(7)

where $n_{tk}$ and $N_k$ are the number of plots and pixels in the stratum that contains $i$. When estimating running means, the sample size and number of pixels within each stratum were the sum of the sample sizes and the number of pixels for the periods considered in the running mean.

We assumed that the selection of FIA plots happened independently and computed pairwise inclusion probabilities for two elements, $i\ne j$ in the sample as

$${\pi }_{ijt}={\pi }_{it}{\pi }_{jt}$$

(8)

When $i=j$ then ${\pi }_{iit}={\pi }_{it}$. The total number of FIA plots, and the number of intensification plots for each 10-year period are provided in Table 2.

2.5.2. Point and Variance Estimators

PS and EPS estimators are obtained from the general formula for a GREG estimator

$${\widehat{\overline{\mathrm{\Delta }AGB}}}_t=\frac{1}{N}{\displaystyle \sum }_{i=1}^N{\widehat{\mathit{\beta }}}_{\mathit{t}}{}^t{\mathit{z}}_{it}+\frac{1}{N}{\displaystyle \sum }_{i\in s_t}\frac{\mathrm{\Delta }AG{B}_{it}-{\widehat{\mathit{\beta }}}_{\mathit{t}}{}^t {\mathit{z}}_{it}}{{\pi }_{it}} ,$$

(9)

The first summand in Equation (9) provides the mean value of $\mathrm{\Delta }AGB$ predicted by the assisting model and the second summand is a bias correction factor that accounts for possible misspecifications of the assisting model. For PS and EPS models verifying Equations (4) and (5), the second summand equals zero, and the estimator of $\mathrm{\Delta }AGB$ for the state can be simplified to

$${\widehat{\overline{\mathrm{\Delta }AGB}}}_t={\displaystyle \sum }_{k=1}^p{W}_{kt}{\widehat{\beta }}_{kt}$$

(10)

where $W_{kt}=\frac{N_{kt}}{N}$ is the $k^{th}$ stratum weight, with $N_{kt}$ the number of units in the population mapped by $f$ to stratum $k$ in the considered period. The coefficients ${\widehat{\beta }}_{kt}$ are the Hajek estimators of the stratum means of $\mathrm{\Delta }AGB$ for the decade extending from year $t$ to year $t+10$.

$${\widehat{\beta }}_{kt}=\frac{{\widehat{\mathrm{\Delta }AGB}}_{kt}}{{\widehat{N}}_{kt}}.$$

(11)

The numerator in Equation (11) is

$${\widehat{\mathrm{\Delta }AGB}}_{kt}={\displaystyle \sum }_{i\in s_t}\frac{\Delta AG{B}_{it}}{{\pi }_{it}}{z}_{ikt},$$

(12)

the HT estimator of the total $\mathrm{\Delta }AGB$ for the 10-year period in the $k^{th}$ stratum. The denominator of Equation (11) is

$${\widehat{N}}_{kt}={\displaystyle \sum }_{i\in s_t}\frac{1}{{\pi }_{it}}{z}_{ikt}$$

(13)

is the HT estimator of the stratum size.

A Taylor approximation for the variance of the estimator in Equation (10) is

$$V\left({\widehat{\overline{\mathrm{\Delta }AGB}}}_t\right)=\frac{1}{N^2}{\displaystyle \sum }_{i\in U}{\displaystyle \sum }_{j\in U}{\mathrm{\Delta }}_{ijt}\frac{e_{it}^0{e}_{jt}^0}{{\pi }_{it}{\pi }_{jt}}$$

(14)

where $e_{it}^0=\mathrm{\Delta }AG{B}_{it}-{\widehat{\mathit{\beta }}}_{0t}{}^t{\mathit{z}}_{it}$, and ${\widehat{\mathit{\beta }}}_{0t}{}^t{\mathit{z}}_{it}$ is the prediction of $\mathrm{\Delta }AG{B}_{it}$ if the assisting model was fit using the data for the entire population [7] p. 237. The variance in Equation (14) can be estimated using the variance estimator

$$\widehat{V}\left({\widehat{\overline{\mathrm{\Delta }AGB}}}_t\right)=\frac{1}{N^2}{\displaystyle \sum }_{iϵs_t}{\displaystyle \sum }_{jϵs_t}\left[\frac{{\pi }_{ijt}-{\pi }_{it}{\pi }_{jt}}{{\pi }_{ijt}}\left(\frac{g_{it}{e}_{it}}{{\pi }_{it}}\right)\left(\frac{g_{jt}{e}_{jt}}{{\pi }_{jt}}\right)\right]$$

(15)

where $e_{it}=\mathrm{\Delta }AG{B}_{it}-{\widehat{\mathit{\beta }}}_{\mathit{t}}{}^t{\mathit{z}}_{it}$ and $g_{it}=\frac{N_{gt}}{{\widehat{N}}_{gt}}$ [7] p. 264. From Equation (8), we obtain ${\pi }_{ijt}-{\pi }_{it}{\pi }_{jt}=0$, if $i\ne j$ and Equation (15) reduces to:

$$\widehat{V}\left({\widehat{\overline{\mathrm{\Delta }AGB}}}_t\right)=\frac{1}{N^2}{\displaystyle \sum }_{iϵs_t}\left[\left(1-{\pi }_{it}\right){\left(\frac{N_{kt}{e}_{it}}{{\widehat{N}}_{kt}{\pi }_{it}}\right)}^2\right]$$

(16)

Estimates of the running mean for the periods $t_0$ to $t_0+10$ through $t_0+a$ to $t_0+a+10$ and their associated variances were obtained by extending the sums in Equations (12), (13), (16), (17), and (18) to the elements of $s_{t_0{t}_0+a}$ and replacing $N$ by $aN$.

2.5.3. Comparison to Current PNW-FIA Estimators and Horvitz-Thompson Estimators

Improvements in precision with respect to the current PNW-FIA practices were assessed by comparing the EPS estimators developed in this study to those resulting from FIA-PS. PS and EPS estimators differ in how the post-strata are obtained, but once the post-strata are created, Formulas (9) though (18) are used for both types of stratifications. To analyze improvements in precision with respect to a situation in which no auxiliary information is available, we compared the developed EPS estimators and FIA-PS estimators to the HT estimator

$${\widehat{\overline{\mathrm{\Delta }AGB}}}_{\pi t}=\frac{1}{N}{\displaystyle \sum }_{i\in s_t}\frac{\mathrm{\Delta }AG{B}_{it}}{{\pi }_{it}}$$

(17)

In this case, no strata were used, and adjustments of inclusion probabilities were performed at the EU level. Thus, ${\pi }_{it}$ were computed as ${\pi }_{it}=\frac{n_{t EU}}{N_{EU}}$ with $N_{EU}$ and $n_{t EU}$ the number of pixels and sample size in the EU where unit $i$ is located. The variance of the HT estimator was computed as

$$\widehat{V}\left({\widehat{\overline{\mathrm{\Delta }AGB}}}_{\pi t}\right)=\frac{1}{N^2}{\displaystyle \sum }_{iϵs_t}\left[\left(1-{\pi }_{it}\right){\left(\frac{\mathrm{\Delta }AG{B}_{it}}{{\pi }_{it}}\right)}^2\right]$$

(18)

For a given PS or EPS estimator, $M$, where $M$ can be GREG-EPS, GREG-EPS-CMS, GL-EPS, GL-EPS-CMS, TREE-EPS, TREE-EPS-CMS and FIA-PS, the relative efficiency for the period starting year $t$ was computed as

$$\mathrm{\Delta }{\widehat{V}}_{Mt}=1-\frac{\widehat{V}\left({\widehat{\overline{\mathrm{\Delta }AGB}}}_{Mt}\right)}{\widehat{V}\left({\widehat{\overline{\mathrm{\Delta }AGB}}}_{\pi t}\right)}$$

(19)

This metric provides an estimate of the reduction in the variance that the estimator $M$ produces with respect to the HT estimator in terms relative to the variance of the latter. Comparisons with FIA-PS were analyzed by observing the differences between $\mathrm{\Delta }{\widehat{V}}_{Mt }$ and $\mathrm{\Delta }{\widehat{V}}_{FIA-PSt}$.

3. Results

3.1. EPS and PS Assisting Models and Summaries

The most important differences between EPS methods and FIA-PS were that EPS methods tended to explain a larger amount of the variance of $\mathrm{\Delta }AGB$, but created a smaller number of strata (Figure 2). Adjusted coefficients of determination, Adj-R², for the assisting models for EPS methods ranged between 35.93% to 40.19%, while R² for FIA-PS was 20.45%. Root mean squared errors, RMSE, for the assisting models ranged from 3.66 Mg ha⁻¹ year⁻¹ to 3.19 Mg ha⁻¹ year⁻¹ for EPS methods, while RMSE for FIA-PS was 4.19 Mg ha⁻¹ year⁻¹. The proportion of plots with observed and predicted values of $\mathrm{\Delta }AGB$ with the same sign (i.e., positive changes predicted as positive and negative changes predicted as negative changes) ranged from 84.25% to 84.94% for EPS methods. For FIA-PS, this percentage was 74.13%. These results indicate that the EPS variants analyzed in this study can improve the precision of the current estimates of $\mathrm{\Delta }AGB$. GL-EPS was the method that experienced the largest improvements in terms of Adj-R² when $\mathrm{\Delta }CM{S}_t$ was added to the stack of auxiliary variables, Adj-R² was 36.38% for GL-EPS and 40.19% for GL-EPS CMS. Adding $\mathrm{\Delta }CM{S}_t$ to the pool of auxiliary variables had a minor effect on GREG-EPS (i.e., Adj-R² was 37.15% for GREG-EPS and 38.38% for GREG-EPS-CMS). The model for TREE-EPS had an Adj-R² of 36.25% and explained slightly more variance than the model for TREE-EPS-CMS, with an Adj-R² of 35.93%.

Assisting models for EPS methods were substantially more parsimonious than the region’s current PS. The EPS model with the most strata was the model for GREG-EPS-CMS. This model created 97 strata, approximately one-half of the number of strata of FIA-PS (

Table 3. Total number of strata, number of sampled strata, and proportion of area sampled by 10-year periods and for all 10-year periods combined.

		2001–2011	2002–2012	2003–2013	2004–2014	2005–2015	2006–2016	2007–2017	2008–2018	All Periods
GREG-EPS	Total # of strata	84
	Sampled strata	84	84	83	84	83	84	84	84	84
	% area sampled	100.00	100.00	99.89	100.00	99.80	100.00	100.00	100.00	100.00
GREG-EPS-CMS	Total # of strata	97
	Sampled strata	95	97	97	97	97	97	94	95	97
	% area sampled	99.47	100.00	100.00	100.00	100.00	100.00	99.07	99.51	100.00
GL-EPS	Total # of strata	30
	Sampled strata	21	20	19	21	20	18	18	20	27
	% area sampled	99.75	99.87	99.80	99.85	99.74	99.83	99.56	99.78	100.00
GL-EPS-CMS	Total # of strata	30
	Sampled strata	27	23	24	24	23	21	23	25	30
	% area sampled	99.65	99.30	99.60	99.76	99.60	99.77	99.53	99.71	100.00
TREE-EPS	Total # of strata	44
	Sampled strata	44	44	44	44	44	44	44	44	44
	% area sampled	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00
TREE-EPS-CMS	Total # of strata	34
	Sampled strata	34	34	34	34	34	34	34	34	34
	% area sampled	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00
FIA-PS	Total # of strata	191
	Sampled strata	167	173	175	177	171	166	163	167	191
	% area sampled	96.79	97.68	97.64	97.43	97.50	96.61	95.40	96.92	100.00

). The most important similarities between the spatial representation of the EPS developed in this study were that: (1) strata close to the coast, and Willamette valley were the ones that had larger and positive stratum means, (2) strata with large negative changes appeared scattered on the western third of the state (west of the cascades mountains) where most intense industrial forest practices occur and (3) the eastern side of the state was characterized by stratum means that were close to zero, with the exception of the areas on the Ochoco National Forest and the Blue and Wallowa mountains located in northeast Oregon. The most important differences between EPS methods are derived from the number of strata they generated (Table 3). GREG-EPS and GREG-EPS-CMS were the methods that resulted in more strata, and associated strata maps provided a finer-grained representation of the variability of $\mathrm{\Delta }AGB$ (Figure 3). GL-EPS and GL-EPS-CMS were the most parsimonious EPS models, and the resulting strata maps captured large-scale trends but showed low variability at local scales (Figure 3). Finally, TREE-EPS and TREE-EPS-CMS resulted in stratifications with an intermediate level of detail (Figure 3).

EPS models used both static and dynamic auxiliary variables. The collapsing steps for GREG-EPS and GREG-EPS-CMS were driven by the variables $CATACDIS{T}_t$, $ECO$ and $OWN$. However, these variables alone do not justify the final number of strata which indicates that other variables play a role in the stratification. The model for GL-EPS used the variables $ELEV$ and $PCPI$, $ACDIS{T}_t$ and $\mathrm{\Delta }NLC{D}_t$, and the model for GL-EP-CMS used $ELEV$, $PCPI$, $ACDIS{T}_t$ and $\mathrm{\Delta }CM{S}_t.$ All variables except $MAXFSE{V}_t$, were used to form at least one node of the partitioning tree of TREE-EPS, but the variables $\mathrm{\Delta }NLC{D}_t$ and $ACDIS{T}_t$ determined the majority of the splits. For TREE-EPS-CMS, all variables except $ACDIS{T}_t$ and $MAXFSE{V}_t$ appeared in at least one split of the partitioning tree (see supplementary material Supplementary S2). $\mathrm{\Delta }CM{S}_t$ and $\mathrm{\Delta }NLC{D}_t$ determined most of the splits in the binary tree with $\mathrm{\Delta }CM{S}_t$ playing a role such as that played by $ACDIS{T}_t$ in TREE-EPS. The partitioning tree for EU WL was the same for TREE-EPS and TREE-EPS-CMS. The only variable that was never used in the models for TREE-EPS, TREE-EPS-CMS, GL-EPS or GL-EPS-CMS was $MAXFSE{V}_t$ which indicates that it provides information that is redundant with $ACDIS{T}_t$ and $\mathrm{\Delta }CM{S}_t,$ the other auxiliary variables derived from spectral data from Landsat images.

GREG-EPS, GREG-EPS-CMS, GL-EPS, GL-EPS-CMS, and FIA-PS resulted in micro-strata with a minimal presence in the territory and that were not sampled. This problem was most relevant for FIA-PS, where the proportion of area corresponding to unsampled strata ranged from 2.3% to 4.6% (Table 3). The collapsing process minimized but did not eliminate the problem of creating empty strata for GREG-EPS and GREG-EPS-CMS, where the area of the unsampled strata reached 0.2% and 0.9%, respectively (Table 3). For GL-EPS and GL-EPS-CMS, unsampled strata appeared in all EUs. Their area representation was less than 0.6% and 0.9% of the total area, respectively (Table 3). The recursive partitioning algorithm of TREE-EPS avoids creating empty strata and incorporates a tuning parameter that ensures that strata have a minimum sample size specified by the modeler. A minimal fraction of area corresponding with no-data values in some of the auxiliary variables used for TREE-EPS and TREE-EPS-CMS resulted in areas without stratum assignments. However, the proportion of those areas was less than 0.001% of the total area (Table 3). The areal extent of the unsampled strata decreased when considering multiple periods for running means. All strata were sampled for all methods except GL-EPS and GL-EPS-CMS when considering the eight periods combined. For GL-EPS and GL-EPS-CMS, the area that the unsampled strata represented in the territory was reduced to less than 0.05%.

3.2. Estimates of Changes for the State for Specific 10-Year Periods

All methods, including the HT estimator and the FIA-PS, estimated that total $\mathrm{\Delta }AGB$ (per unit area and year) for all periods under consideration was positive (Figure 4). This clearly indicates a steady increase in forest AGB in the region. Estimated change ranged from 0.33 Mg ha⁻¹ year⁻¹ to 0.88 Mg ha⁻¹ year⁻¹ depending on the 10-year period and estimator. Differences between estimated totals for different periods tended to be of larger magnitude than differences between methods for a given period (Figure 4).

Estimated improvements in precision with respect to the HT estimator, $\mathrm{\Delta }{\widehat{V}}_M$, varied across model-assisted estimators and also across periods; however, some general trends were observed. All model-assisted estimators, including the PS currently used by PNW-FIA, consistently improved the precision of the HT estimator (Figure 5). The only exception to this trend was 2007–2017 for all EPS variants except for GREG-EPS and the periods 2007–2017 and 2008–2018 for FIA-PS. The low value of $\mathrm{\Delta }{\widehat{V}}_M$ for 2007–2017 that was obtained for all EPS methods and FIA-PS is due to an unusually low variance of the HT estimator for this period.

The median improvement obtained by the EPS and PS estimators with respect to the HT estimator indicated that GL-EPS-CMS and GL-EPS produced larger gains in precision than other EPS methods. The median improvement for GL-EPS-CMS and GL-EPS was 48.02% and 47.99%. These methods were followed by GREG-EPS-CMS and GREG-EPS, TREE-EPS-CMS, and TREE-EPS, with improvements in the efficiency of around 40%. FIA-PS produced a median improvement with respect to the HT estimator of 28.24%. EPS variants were generally more efficient than FIA-PS (

Table 4. Median and mean improvement of EPS and PS estimators for the HT estimator.

Method	$\mathbf{Median} \mathbf{of} \mathrm{\Delta }{\hat{\mathit{V}}}_{\mathit{M}}$ for 10-Year Periods	$\mathbf{Mean} \mathbf{of} \mathrm{\Delta }{\hat{\mathit{V}}}_{\mathit{M}}$ for 10-Year Periods
GREG-EPS	37.95%	36.25%
GREG-EPS-CMS	40.13%	35.91%
GL-EPS	48.02%	38.36%
GL-EPS-CMS	47.99%	42.27%
TREE-EPS	41.05%	31.66%
TREE-EPS-CMS	40.27%	30.16%
FIA-PS	28.24%	20.68%

), but FIA-PS outperformed EPS methods in at least one of the 10-year periods (Figure 5).

3.3. Estimators of Running Means

When considering running means of more than one 10-year period, the patterns we observed were similar to those observed for single 10-year periods (Figure 6). All methods consistently estimated positive values of $\mathrm{\Delta }AGB$ and EPS methods and FIA-PS resulted in gains in performance with respect to the HT estimator. With only a few exceptions, EPS methods also resulted in gains in performance with respect to FIA-PS (Figure 6). As expected, when increasing the number of 10-year periods included in the running means, differences between the estimates and associated uncertainties provided by the tested methods decreased at a rate that was approximately inversely proportional to the square root of the combined sample size. Running means for the eight 10-year periods ranged between 0.58 Mg ha⁻¹ year⁻¹ and 0.65 Mg ha⁻¹ year⁻¹, and confidence intervals for the estimates provided by any method contained the estimates provided by the other methods (Figure 7). For the running means of eight 10-year periods, adding $\mathrm{\Delta }CM{S}_t$ to the pool of auxiliary variables consistently resulted in values of $\mathrm{\Delta }{\hat{V}}_M$ larger than those obtained by EPS models not using this auxiliary variable. Differences between GL-EPS and GL-EPS-CMS were largest at 6.7%, i.e., 38.3% vs. 45.0%, and differences between GREG-EPS and GREG-EPS-CMS and TREE-EPS and TREE-EPS-CMS were of small magnitude < 1%, i.e., 36.1% vs. 36.23% and 33.95 vs. 34.2%, respectively.

4. Discussion

4.1. Similarities between Model-Assisted Estimators

All methods estimated increases in forest AGB for the periods 2001–2011 to 2008–2018 (Figure 4). This indicates that despite increasing large-scale wildland fires and insect outbreaks (i.e., bark beetle) events over the last two decades, there is a net accumulation of forest AGB in the region. Depending on the method, the average estimated change over the eight 10--year periods analyzed in this study ranged between 0.58 Mg ha⁻¹ year⁻¹ and 0.65 Mg ha⁻¹ year⁻¹. These results indicate a consistent increase in forest AGB in the study area but do not inform about the type of changes that are occurring. The accumulation of forest AGB can result from many processes of a very different nature and with very different ecological implications. For example, both: (1) increases in the amount of land sustaining mature forest structures with larger trees and more fire resilient structures, and (2) increases in the amount of young and dense forest structures resulting from fire suppression efforts contribute to total forest AGB accumulation. The ecological implications of each type of increase mentioned above are very different and could be evaluated by a more detailed analysis of the field-measured plot data. However, PS and EPS allow constructing additive tables for non-overlapping categories in a straightforward manner. Thus, these methods can be used to derive estimates of change by typologies.

EPS estimators developed in this study use a method closely related to the direct approach proposed by McRoberts et al. (2015) because $\mathrm{\Delta }AGB$ is the modeled variable. However, small differences with the direct approach exist due to the inclusion of static auxiliary variables. Estimators developed in this study or using direct approaches [20,25,26,27] are typically preferred over indirect methods where assisting models are developed for $AGB$ and then estimates for time 2 and time 1 are differenced. Direct methods offer a more straightforward way of modeling change; however, because a different response is modeled and different auxiliary variables are used, the estimated changes from a direct estimator would not match the difference between the AGB totals for time 2 minus those for time 1 of an indirect estimator.

4.2. Model Selection

Models were selected with a dataset that combined all the 10-year periods under consideration. Once models were selected, they were re-fitted with the data for a specific 10-year interval or intervals in the case of running means. This implies trade-offs between consistency, performance, and operationality, which can be summarized in the following remarks. Identifying unique models for all 10-year periods combined is operationally simpler than developing models for each period separately, but the resulting models are not necessarily optimal for any particular 10-year interval. Although this lack of optimality can certainly occur, all model-assisted estimators outperformed the HT and the PS estimators currently used by PNW-FIA. When the model selection is developed with the combined sample of eight 10-year periods, estimating $\mathrm{\Delta }AGB$ for a specific 10-year interval is similar to obtaining estimates for a subpopulation of a larger one, the one resulting from the union of the populations of the eight 10-year periods. Differences in the relationships between auxiliary variables and $\mathrm{\Delta }AGB$ are not expected to vary substantially across periods, and obtaining stratum-specific means and variances for each 10-year period separately should minimize potential bias problems.

4.3. Differences between Model-Assisted Estimators

In terms of performance, EPS methods exceed the HT estimator and FIA-PS. The development of a model-assisted estimator for an NFI implies evaluating the estimator’s performance but also several factors related to their practical implementation. We found that differences in performance between EPS methods were small and decreased when considering running means. Thus, we envision that operational constraints will likely outweigh performance considerations when selecting an EPS method to produce official stratification tables. One factor that we found especially challenging for the development of EPS was that the chances of finding micro-strata that need to be collapsed increased very rapidly when adding auxiliary variables. Model selection methods such as elastic net facilitated EPS development and reduced the need for strata-collapsing processes. However, the application of these methods, except for TREE-EPS, was far from being completely automatic.

For GREG-EPS, directly performing a model selection step on a saturated model, including all categorical variables and their interactions, produced a very large number of strata, many of them with a very small or zero sample size. It was necessary to perform an initial pre-collapsing stage that grouped similar interactions. These groups were considered in the elastic-net selection procedure, and the resulting models required a final collapsing of groups discarded by the elastic-net method. The need for collapsing steps for this method can be considered its main disadvantage. The way in which categorial variables are defined (i.e., the number of categories per variable) and the number of auxiliary variables need to be carefully considered as the number of initial groups that need to be collapsed grows very quickly. This task involves decisions by the modeler and involves some degree of subjectivity. In addition, a direct consequence of the need for a collapsing step is that interpretability of the generated strata is less straightforward than for TREE-EPS or GL-EPS. After the collapsing process, strata can be interpreted as groups of pixels with similar auxiliary information. Still, these groups are likely to have some internal variability, and producing a full description of the strata requires backtracking a large number of initial combinations of categorical variables that were merged during the collapsing process.

The definition of strata for GL-EPS is relatively straightforward, and the strata themselves have a meaningful interpretation, which is an advantage of this method. An important difficulty of GL-EPS was finding meaningful thresholds of the variable of interest that ensured that no empty strata were created. This problem was discussed by [12] and is accentuated when strata are forced to belong to only one EU. We used broad categories of $\mathrm{\Delta }AGB$, and still some strata were unsampled for specific 10-year periods. The areal representation of unsampled strata was small, below 1% in all considered periods, but remained even when considering the running mean of all the 10-year periods (Table 3). Unsampled strata tended to correspond with the extreme intervals defined for $\mathrm{\Delta }AGB$ and could be merged to the neighboring strata. However, this solution implies that the final thresholds that give rise to the final strata are not defined a priori, which is the case considered in the definition of GL-EPS [11].

The simple GL-EPS version used in this study uses the identity function for $h$, which explains why observed improvements in performance when adding $\mathrm{\Delta }CM{S}_t$ were larger for this method than for the other EPS variants. Using the identity function for $h$ makes GL-EPS utterly reliant on the linear predictor ${\gamma }_t^t{x}_{ti}$. The auxiliary variable $\mathrm{\Delta }CM{S}_t$ was computed from validated AGB predictions and its association with $\mathrm{\Delta }AGB$ improves the linear predictor. Other EPS variants experienced small or negligible improvements when using $\mathrm{\Delta }CM{S}_t$. Because of their greater ability to model nonlinear relationships, they were able to find associations in the auxiliary information database that explained most of the variability explained by $\mathrm{\Delta }CM{S}_t$. The design consistency of EPS estimators, similar to GL-EPS in their construction, but using non-parametric models instead of a linear predictor ${\gamma }_t^t{x}_{ti}$ has been proven [12]. Further research must explore these non-parametric versions of GL-EPS estimators because more flexible model specifications will likely capture non-linear relationships between auxiliary information and $\mathrm{\Delta }AGB$. The simple specification used for GL-EPS is not necessary a disadvantage. The linear predictor provides pixel level predictions ${\gamma }_t^t{x}_{ti}$ of $\mathrm{\Delta }AGB$ and this intermediate output is interesting on its own. It provides a 30 m resolution map of $\mathrm{\Delta }AGB$ that can be used for multiple purposes, including the computation of GREG estimators similar to those in [6,20,28].

Finally, TREE-EPS clearly outperformed all other alternatives for PS when considering easiness of practical implementation. By design, the recursive partitioning algorithm from Toth and Eltinge (2011) ensured that no empty strata were created. The implementation of this method in the rpms R package [23] allows setting minimum sample sizes for the final strata, which allows for fine-tuning of the PS process and minimizes the chances of obtaining unsampled strata in specific 10-year intervals. Furthermore, TREE-EPS can model non-linear relationships, eliminating over-smoothing problems observed for GL-EPS. In addition, it easily accommodates large numbers of auxiliary variables of mixed type (i.e., continuous and categorical), and the assisting model itself is a binary classification tree with a straightforward interpretation (see supplementary material Supplementary S2).

4.4. Other Considerations of Practical Importance

4.4.1. Estimation of Variance and Estimation of Change for Periods Not Matching the PNW-FIA Panel Frequency

PNW-FIA, and most NFI globally, rely on systematic or near-systematic sampling designs that make it impossible to obtain unbiased variance estimators [29,30]. The standard procedure is to either treat systematic samples using formulas for simple random sampling or treat the elements in the sample as if they proceeded from independent draws, as we did in this study. Computing variances in any of these ways can lead to bias problems, especially when the response shows strong trends. A priori, this problem is less important for model-assisted estimators that explain a significant portion of the response variability because weaker trends should be present in the model residuals needed for variance estimation (e.g., Equation (15)). Nevertheless, the problem is not eliminated, and alternative variance estimators have been proposed to account for spatial correlations when estimating variances from systematic samples [31,32,33]. Recent studies have shown that using these estimators with FIA data can improve the efficiency when estimating variances of spatially correlated attributes, even after post-stratification [34].

Estimation of change for periods that do not match the panel cycle is another topic that deserves attention (e.g., estimation of changes for 3-year period). One feature of the rotating panel design of FIA is that panels do not overlap. This enforces the change estimation for lags that are not multiples of the rotating panel period to rely on differencing estimators of $AGB$ totals for state at times 1 and 2. If ${\widehat{\overline{AGB}}}_{\mathrm{t}1}$, and ${\widehat{\overline{AGB}}}_{\mathrm{t}2}$ are the estimators for the totals at time t1 and t2, with the difference between t1 and t2 not a multiple of the panel period, then the estimated change is ${\widehat{\overline{\mathrm{\Delta }AGB}}}_{\mathrm{t}1,\mathrm{t}2}={\widehat{\overline{AGB}}}_{\mathrm{t}2}-{\widehat{\overline{AGB}}}_{\mathrm{t}1}$, and its variance is

$$V\left({\widehat{\overline{\mathrm{\Delta }AGB}}}_{t1,t2}\right)=V\left({\widehat{\overline{AGB}}}_{t2}\right)+V\left({\widehat{\overline{AGB}}}_{t1}\right)-2Cov\left({\widehat{\overline{AGB}}}_{t2}, {\widehat{\overline{AGB}}}_{t1}\right)$$

(20)

From Equation (20), one can clearly see that a way to obtain more reliable estimators of change for a lag that is not a multiple of the rotating panel period is to minimize $V\left({\widehat{\overline{AGB}}}_{\mathrm{t}1}\right)$ and $V\left({\widehat{\overline{AGB}}}_{\mathrm{t}2}\right)$ [35]. This can be achieved by applying the techniques that we used for $\mathrm{\Delta }AGB$ to $AGB$: however, applying different models for different times can result in strata with cumbersome interpretations. Variance estimators for $V\left({\widehat{\overline{\mathrm{\Delta }AGB}}}_{\mathrm{t}1,\mathrm{t}2}\right)$ involve assumptions similar to those previously discussed. In particular, the estimation of $V\left({\widehat{\overline{AGB}}}_{\mathrm{t}1}\right)$ and $V\left({\widehat{\overline{AGB}}}_{\mathrm{t}2}\right)$ are instances of the same problem, and the same unknowns and potential solutions apply to them. Methods to provide spatially balanced samples with a time overlap have been proposed [36] and could be used to estimate $2Cov\left({\widehat{\overline{AGB}}}_{\mathrm{t}2}, {\widehat{\overline{AGB}}}_{\mathrm{t}1}\right)$. However, the implementation of those solutions would require prohibitively costly and deep modifications of the FIA sampling design.

4.4.2. Auxiliary Variables

Finally, one of the most time-consuming steps of this study was compiling a database with auxiliary variables potentially correlated with changes in forest AGB. This step was critical as model-assisted estimators’ performance improvements are directly related to the strength of the correlation between response and auxiliary variables. This study was developed using the state of Oregon as a relatively data-rich test area, but several questions regarding the development of a similar auxiliary information database need to be considered if similar estimators are to be developed in other regions. Static auxiliary variables such as climate normals or topographic indexes are not expected to change substantially over time, and suitable replacements exist for most of these datasets in other areas of the world. Dynamic auxiliary variables require more careful attention. The independence of $\mathrm{\Delta }CM{S}_t$ and $CAT\mathrm{\Delta }CM{S}_t$ of the FIA sampling design makes these datasets very appealing for FIA and can enhance estimates of biomass change or simplify model fitting processes. In general, dynamic auxiliary variables are expected to have increasing relevance in NFIs as new multitemporal datasets become available. Thus, the maintenance of these datasets will be a critical factor in ensuring that NFIs can derive official tables periodically.

5. Conclusions

The main conclusions of this study can be summarized as follows:

EPS methods allowed accommodating static and dynamic auxiliary variables and showed their potential to improve PS methods currently used in FIA for the estimation of $\mathrm{\Delta }AGB$.

All EPS methods had operational advantages and disadvantages, and issues related to implementation, generation of unsampled strata, or interpretability of the generated strata might outweigh performance considerations.

GL-EPS relied on the construction of a linear predictor and adding $\mathrm{\Delta }CM{S}_t$ to the pool of auxiliary variables produced substantial performance improvements. This indicates that $\mathrm{\Delta }CM{S}_t$ can help in the estimation of $\mathrm{\Delta }AGB$. TREE-EPS and GREG-EPS were able to model non-linear relationships that explained most of the variability explained by $\mathrm{\Delta }CM{S}_t$.

Changes in $\mathrm{\Delta }AGB$ in the study were always positive, even when considering running means of multiple 10-year periods. Estimated changes provide average values of change but do not inform about the nature and type of change. However, the generated strata correctly identified known features in the landscape, which indicates that these EPS can play an important role in studies disaggregating $\mathrm{\Delta }AGB$ by typologies.