Addressing nonresponse bias in forest inventory change estimation using response homogeneity classifications

1.   Introduction

Although assessments of changes in forest resources based on time series of remotely sensed data are becoming more technically advanced and commonplace (Cohen et al., 2010; Nguyen et al., 2019; Schleeweis et al., 2020; Woodcock et al., 2020), estimates from forest inventory field plots still often serve as the primary information source for official statistics and many analytical purposes (Zhao et al., 2019; Oswalt et al., 2019). A key aspect of many national forest inventory (NFI) programs is the ability to monitor change in numerous forest resource attributes. For example, change in total forest area (possibly by forest type or stand size class) as well as biomass and carbon attributes often rely on NFI data (Charru et al., 2010; Gschwantner et al., 2016). Accurate estimation of change statistics is often vital for making effective management and policy decisions affecting forest resources. Further, establishment of sound change accounting practices is necessary for participation in programs aimed at reducing deforestation and forest degradation (Park et al., 2013). As NFI often provide the basis for these information needs, statistical rigor in the sampling design is paramount for achieving accurate statistical estimates. Fortunately, there exists considerable literature in this regard due to longstanding implementation of forest inventories in various countries (Bechtold and Patterson, 2005; Tomppo et al., 2011; Zeng et al., 2015). However, it is nearly impossible to avoid nonresponse in sample based NFIs, where no information is obtained for some of the selected sample units primarily due to hazardous conditions or lack of access permission by private landowners (Patterson et al., 2012). In the U.S., there is some evidence for increasing incidence of denied access by private landowners (Westfall et al., 2022a, b), with similar increasing nonresponse by landowners to natural resource-related mail surveys (Connelly and Brown, 2003). Although some research of nonresponse mitigation in forest inventory has been conducted (Fischer et al., 2011; Corona et al., 2014; Fattorini et al., 2013), there is often little operational effort to explicitly address potential biases in the sample due to nonresponse and consequently some bias in the estimates may be present (Goeking and Patterson, 2013; Magnussen et al., 2017). Rather, nonresponse is often assumed to have minor effects on estimation bias due to the infrequent occurrence but increasing rates of nonresponse may negate that assumption.

A key area of interest in change estimation is to quantify effects of disturbance on forest resources. Disturbance events occur across a wide range of intensity and spatial scale (Stueve et al., 2011; Wilson et al., 2019; Fahey et al., 2020). While all disturbances have the potential to alter forest ecosystem structure and function, it is often that relatively large catastrophic events are of most concern due to the impact potential (Edgar et al., 2019; Senf and Seidl, 2021). In a forest inventory context, these severe disturbances can be problematic in that resulting hazardous conditions can preclude field crew visits to plots and thus nonresponse is more prevalent in the disturbance zone (Nelson et al., 2009). As most forest inventories employ equal probability sampling designs, the result is that information from the disturbance zone is underrepresented in the sample outcome. Failure to acknowledge and address this phenomenon and proceeding to apply standard estimation procedures will likely produce inaccurate population estimates of the disturbance effect (Nelson et al., 2009).

In this study, the effect of nonresponse outcomes on estimates of change is examined in the context of a typical NFI post-stratified estimation method and an alternative approach that explicitly incorporates knowledge of factors leading to nonresponse. The analysis accounts for two scenarios: 1) typical NFI nonresponse outcomes in a general population, and 2) a specific population in which a large disturbance occurs. Thus, the objectives are to: 1) increase awareness of nonresponse issues in the context of repeated plot measurements, 2) assess the potential magnitude of nonresponse bias that may be occurring in change estimates derived from forest inventory data, and 3) present a potential framework for using response homogeneity classes in a change estimation context when nonresponse does not occur uniformly at randomly in the population.

2.   Methods

2.1   Data

Initial testing of methods was conducted using a synthetic population where the values for a generic forest attribute response variable (mimicking tree biomass in tonnes per hectare) were known for each population element at two points in time. For each population element at each time point, classifications were assigned as either being forest or nonforest land use status. To facilitate use of post-stratified estimation methods used to improve the precision of estimates, three post-strata were created and each population element was assigned to a post-stratum. In practice, maps constructed from classified remote sensing data are often used to create the post-strata and determine the post-stratum weights (McRoberts et al., 2002). Post-stratum weights for the synthetic population were determined from the proportions of population elements in each post-stratum. The response variable for each forested population element at T1 was generated based on the post-stratum assignment (h) and a N(μ_h, σ_h) normal distribution where the distributions were 1) μ₁ ​= ​25 and σ₁ ​= ​5, 2) μ₂ ​= ​100 and σ₂ ​= ​12.5), and 3) μ₃ ​= ​200 and σ₃ ​= ​25 for post-strata h ​= ​1 to 3 (Table 1). The population elements having no forest area had the response value set to zero.

Table  1.  Description of synthetic NFI population construction by post-stratum and forest/nonforest population element status at both points in time.

Post-stratum	Time 1				Time 2
	Population element	Population (%)	Observation distribution		Population element	Population (%)	Observation distribution
1	Forest	10.5	N (25, 5)		Forest	9.4	T1 + (5%–10%)
					Nonforest	1.1	0
1	Nonforest	24.5	0		Forest	1.2	N (50, 10)
					Nonforest	23.3	0
2	Forest	40.0	N (100, 12.5)		Forest	36.0	T1 + (5%–10%)
					Nonforest	4.0	0
2	Nonforest	10.0	0		Forest	4.0	N (50, 10)
					Nonforest	6.0	0
3	Forest	9.6	N (200, 25)		Forest	8.6	T1 + (5%–10%)
					Nonforest	1.0	0
3	Nonforest	5.4	0		Forest	1.1	N (50, 10)
					Nonforest	4.3	0

 | Show Table

DownLoad: CSV

Population element response values at T2 were generated using the following methods regardless of post-stratum: 1) values of the response variable at T2 for forested elements that were nonforest at T1 were obtained randomly from a N(50, 10) distribution as areas reverting to forest over the measurement interval would likely have relatively small values, 2) a 'growth' factor approach was used for elements that were forest at both T1 and T2, where the T1 value was increased by a randomly selected value between 5% and 10% to obtain the T2 amount, and 3) nonforest elements were always zero-valued by definition (Table 1). The elements forested at T1 that became nonforest at T2 had a negative change equal to the T1 value. As such, these elements that changed status between T1 and T2 tended to have relatively large negative change amounts (Tables 1 and 2). Overall, there was approximately a 2% reduction in the population total between T1 and T2.

Table  2.  Synthetic NFI population post-stratum RHG weights ($W_{hg}$), number of population elements ($N_{hg}$), total change (T2 – T1), and percentage of office plots (NF) at T1 and T2. Overall values for the post-strata are also given for weight ($W_h$) and population size ($N_h$), total change (T2 – T1), and percentage of office plots (NF) at T1 and T2.

Post-stratum	RHG	$W_{hg}$a	$N_{hg}$	T2 – T1	Office plots (%)
					T1	T2
1	NF–NF	0.665	46, 561	0	100.0	100.0
1	NF–F/F–NF	0.066	4603	‒21, 946.8	53.8	46.2
1	F–F	0.269	18, 836	35, 370.3	0	0
1		0.350	70, 000	13, 423.6	70.1	69.6
2	NF–NF	0.120	11, 990	0	100.0	100.0
2	NF–F/F–NF	0.160	15, 988	‒701, 744.9	49.7	50.3
2	F–F	0.720	72, 022	540, 503.1	0	0
2		0.500	100, 000	‒161, 241.8	19.9	20.0
3	NF–NF	0.289	8657	0	100.0	100.0
3	NF–F/F–NF	0.138	4136	‒368, 979.1	52.2	47.8
3	F–F	0.574	17, 207	258, 423.6	0	0
3		0.150	30, 000	‒110, 555.5	36.1	35.4
a The RHG weights shown here reflect the proportion of observations in the stratum population. In practice, the actual RHG weights within strata are estimated from the selected sample.

 | Show Table

DownLoad: CSV

A sample of the population elements is considered to be analogous to plots that may be selected in a forest inventory sample. Once the sample is drawn, complete remeasurement sampling is assumed where each sample plot is remeasured periodically (Schreuder et al., 2004). Selected plots were designated for a field visit when forestland was present or are considered as an 'office' plot where no field visit was required due to no forest area being present. This corresponds to sampling methodology of the United States NFI, wherein sample plots having zero probability of forest condition (e.g., open water, cropland, developed impervious surface) as determined by office staff via photointerpretation of high resolution aerial imagery, are not assigned for field measurement. As office-determined nonforest plots are not subject to nonresponse, the only plots used in estimation that had zero probability of nonresponse were those that were nonforest office plots at both time points (Note that occasionally a field plot may be nonforest, but it is still subject to potentially being nonresponse). Thus, only the field plots were subject to either partial or complete nonresponse.

To account for different nonresponse probabilities, response homogeneity groups (RHG) were created within each post-stratum (Table 2). In the context of change estimation, plots can be classified in four ways based on their forest/nonforest status at each time point: 1) forested at T1 and nonforest at T2 (F–NF), 2) nonforest at T1 and forest at T2 (NF–F), 3) forest at both times (F–F), and 4) nonforest at both times (NF–NF). Only plots in the latter classification (NF–NF) are assured to have no nonresponse, as these plots would have been nonforest office plots for both measurements where both T1 and T2 forest attribute values are zero and thus the change is also zero. The other three categories would have been subject to nonresponse at either T1 (F–NF), T2 (NF–F), or both (F–F), with the latter being considered the most likely to incur nonresponse due to access requirements at both measurements. Thus, three RHGs are created to account for differing nonresponse probabilities (Särndal et al., 1992; section 15.6.2): 1) no nonresponse possible (NF–NF), 2) nonresponse possible at only one of two measurements (F–NF, NF–F), and 3) potential nonresponse at both measurements (F–F).

A second synthetic population was created to illustrate a specific case where forest disturbance affects nonresponse. In this situation, it is considered that the disturbance creates hazardous conditions such that plots within the disturbance zone are more prone to nonresponse than those elsewhere. In this case, consider a smaller population with 50, 000 elements that is composed entirely of forest land at both T1 and T2, e.g., a large publicly owned forest where access to field plots is always granted. Three post-strata are created and the generic response values for population elements at T1 are generated from N(μ_h, σ_h) distributions where the distributions were N(75, 7.5), N(125, 10), and N(200, 25) for post-strata h ​= ​1 to 3, respectively. As done previously, T2 values were obtained from the T1 values using a randomly selected increase between 5% and 10%. Additionally, an extreme disturbance occurring between T1 and T2 caused losses in the generic response variable on 10, 000 population elements that define the disturbance zone. It was assumed the T2 value of the generic response variable was randomly reduced by 80%–100% (i.e., 0–20% remaining) for each population element subject to the disturbance (Table 3), simulating severe damage levels typical of high intensity disturbances such as wildfires, derechos, or hurricanes (Boose et al., 1994). Overall, this structure results in an approximately 12% reduction in the population total from T1 to T2. Under this scenario, RHGs within post-strata are created to differentiate between population elements within and external to the disturbance zone (Table 4). In practice, this differentiation would likely be accomplished via map information that delineates the disturbance area.

Table  3.  Description of synthetic disturbance population construction by post-stratum and disturbance zone at both points in time.

Time 1					Time 2
Post-stratum	Zone	Population (%)	Plot observation		Zone	Population (%)	Plot observation
1	Undisturbed	20.0	N (75, 7.5)		Disturbed	4.1	T1 + (5%–10%)
					Undisturbed	15.9	[T1 + (5%–10%)] × (0–20%)
2	Undisturbed	30.0	N (125, 10)		Disturbed	6.0	T1 + (5%–10%)
					Undisturbed	24.0	[T1 + (5%–10%)] × (0–20%)
3	Undisturbed	50.0	N (200, 25)		Disturbed	10.0	T1 + (5%–10%)
					Undisturbed	40.0	[T1 + (5%–10%)] × (0–20%)

 | Show Table

DownLoad: CSV

Table  4.  Disturbance population post-stratum RHG weights ($W_{hg}$), number of population elements ($N_{hg}$), and total change (T2 – T1) values. Overall values for the post-strata are also given for weight ($N_h$) and population size ($N_{hg}$), and total change (T2 – T1).

Post-stratum	RHG	$W_{hg}$a	$N_{hg}$	T2 – T1
1	DIST	0.203	2026	‒135, 728
1	UNDIST	0.797	7974	44, 878
1		0.200	10, 000	‒90, 850
2	DIST	0.120	2988	‒334, 203
2	UNDIST	0.720	12, 012	112, 775
2		0.500	15, 000	‒221, 428
3	DIST	0.289	4986	‒885, 119
3	UNDIST	0.574	20, 014	299, 910
3		0.150	25, 000	‒585, 209
a The RHG weights shown here reflect the proportion of observations in the stratum population. In practice, the actual RHG weights within strata are estimated from the selected sample.

 | Show Table

DownLoad: CSV

2.2   Analysis

Analytical results were obtained for a typical implementation of post-stratified estimation (Scott et al., 2005) as well as for implementation of RHGs within the post-strata (Goeking and Patterson, 2013; Westfall, 2022). First, the standard FIA post-stratified procedures are outlined, with subsequent extension to accommodate the RHG approach. Initially, plot values to be used in estimation are obtained by adjusting the observed value to account for plots having partial nonresponse within each post-stratum h.

where $a_{hi}$ a is the proportion of measured area for plot i in post-stratum h, $m_h$ is the number of plots (excluding completely nonsampled plots) in post-stratum h, ${\overline{p}}_h$ is the post-stratum h adjustment factor, $y_{hi}$ is the original observed value for plot i in post-stratum h, and $y_{hi}^{\prime }$ is the adjusted value for plot i in post-stratum h. The mean (${\overline{y}}_h^{\prime }$) and estimated variance of the mean ($v\left({\overline{y}}_h^{\prime }\right)$) for each post-stratum h are obtained from simple random sample estimators:

In this implementation of post-stratification, the mean of plots that are completely nonsampled are, in effect, assumed to be equal to the post-stratum mean calculated from the observed plots. Finally, estimates of the population total and corresponding variance are calculated by combining across post-strata. Equation (5) is used to obtain the estimated population total (${\hat{y}}^{\prime }$) using post-stratum mean ${\overline{y}}_h^{\prime }$ values, post-stratum weights ($W_h$), and the population size (N). The estimated variance of the population total ($v\left({\hat{y}}^{\prime }\right)$) includes both sampling uncertainty and the variability due to random post-stratum sample sizes (Scott et al., 2005).

where $m=\sum _{h=1}^H{m}_h$ is the total number of at least partially observed plots.

Implementation of an RHG estimation framework is accomplished by using the RHGs created within each post-stratum, wherein the nonresponse occurrence is considered to be random within each RHG. From an estimation perspective, Westfall et al., 2022a, b suggest the use of ratio-to-size estimators within each RHG to avoid potential bias in estimates of variance due to the use of the adjustment factor method described above:

where ${\overline{y}}_{hg}$ is the mean response for RHG g in post-stratum h, $y_{hgi}$ is the observed value from the ith plot in RHG g for post-stratum h, $a_{hgi}$ is the proportion of the plot that was accessible, and $s_{hg}^2$ is the sample variance in RHG g for post-stratum h. Subsequently, RHG need to be combined to estimate means and variances for each post-stratum. This requires weights to be determined for each RHG, which are based on the proportion of plots selected into the sample for each RHG relative to the number of plots selected for post-stratum h (Goeking and Patterson, 2013).

where $w_{hg}$ is the estimated weight for RHG g in post-stratum h, $n_{hg}$ is the number of plots in RHG g within post-stratum h, and $n_h$ is the number of plots in post-stratum h. It is important here to recognize that n refers to the sample sizes that also include completely nonresponse plots; whereas m notation used earlier excludes plots that were entirely inaccessible. In the context of change estimation, some assumptions must be made in practice to assign nonresponse plots to an RHG for calculation of $w_{hg}$ - acknowledging that nonresponse occurs when the plot is unobserved at either T1, T2, or both. In the case of the U.S. NFI and office/field plot considerations, it may be reasonable to assume that plots observed at one of the two time periods had the same status when unobserved such that these plots would be classified as either F–F or NF–NF accordingly. When a plot is nonresponse at both T1 and T2, it was probably a field plot both times and is likely to be F–F. Of course, exceptions will occur and it must be accepted that RHG assignments will not always be entirely accurate.

Additionally, unlike $W_h$, the $w_{hg}$ are estimated from the sample and thus RHG weights are an additional source of variability. Estimation of the $w_{hg}$ can be considered as a special case of double sampling for post-stratification where the first phase sample consists only of points associated with the population elements selected into the sample. Further, Särndal et al. (1992) propose that variation in nonresponse that would occur among different samples should be included in estimates of uncertainty. In consideration of these factors, it is proposed the estimators presented in Westfall (2022) be used to obtain post-stratum means and variances:

Hereafter, the remainder of the estimation process proceeds as described above in equations (5), (6) by substituting ${\overline{y}}_h$ and $v\left({\overline{y}}_h\right)$ in place of ${\overline{y}}_h^{\prime }$ and $v\left({\overline{y}}_h^{\prime }\right)$, respectively, to obtain RHG-based estimates of the population total $\hat{y}$ and its variance $v\left(\hat{y}\right)$. Construction of 95% confidence intervals was accomplished using ${\hat{y}}^{\prime }$ ​± ​$1.96\sqrt{v\left({\hat{y}}^{\prime }\right)}$ for the typical post-stratified estimates and $\hat{y}$ ​± ​$1.96\sqrt{v\left(\hat{y}\right)}$ for RHG estimates.

Comparisons between estimations approaches were based on Monte Carlo simulations. For the generalized change estimation scenario, a sample of 2000 population elements was drawn randomly at each iteration without regard to post-stratum or RHG assignment. This approach produces samples consistent with a post-stratified design where sample sizes within post-strata are random and proportions of sample plots within post-strata often differ from post-strata weights ($W_h$). For the primary change estimation population, mean nonresponse rates for completely inaccessible sample units were 10% (F–F), 5% (F–NF, NF–F), and 0 (NF–NF). Similarly, partial nonresponse was applied on average to 4% of F–F plots, 2% of F–NF/NF–F, and 0 if NF–NF plots where the nonresponse proportion for each plot was randomly determined from a U(0, 1) uniform distribution. This nonresponse scenario produces roughly 6% of plots selected into the sample being complete nonresponse and 2% of plots having partial nonresponse Sample units subject these nonresponse rates were chosen randomly at each iteration. The mean of R ​= ​20, 000 population total estimates for the typical post-stratified estimates and those where RHGs were employed were compared to the known population total change of Y ​= ​−258, 373.7. Similarly, the means of the estimated variances for both estimation approaches ($v\left(\hat{y}\right)$ and $v\left({\hat{y}}^{\prime }\right)$) were assessed against the variance of the R ​= ​20, 000 population total estimates.

For the disturbance change scenario, a sample of 500 population elements was obtained at each iteration in the same manner as described above. It was assumed that all sample units outside the disturbance zone were accessible. Four nonresponse outcomes were examined within the disturbance zone, where complete/partial nonresponse rates of 20/2, 40/3, 60/4 and 80/4 percent were examined. Sample units subject to either complete or partial nonresponse were selected randomly at each iteration. The nonresponse is applied to the sample units in the disturbance zone to reflect the hazardous conditions that directly result in nonresponse that otherwise would not be present. As before, comparisons were made between mean estimates of change from the R ​= ​20, 000 iterations and the known population change of Y ​= ​−897, 487.0. As with the NFI scenario, formula-based variances were compared with the simulation results.

3.   Results

In the standard NFI change estimation scenario, implementation of the typical post-stratified estimator that assumes nonresponse sample unit observations are on average equal to the post-stratum mean produced an estimated change that was −63182.2 (−24.4%) different than the known population change (Table 5). In comparison, the RHG estimator resulted in an estimate that differed from the known value by only −461.3 (0.18%). The formula-based estimates of variance were consistent with the variance of the simulations estimates for both estimators.

Table  5.  Estimates, percent bias, and 95% confidence interval coverage for RHG and standard post-stratified estimation methods for the NFI and disturbance populations.

Population (CNR/PNR)	RHG method							Standard post-stratification
	Estimate			Bias (%)		95% CI coverage (%)		Estimate			Bias (%)		95% CI coverage (%)
	$\hat{y}$	$v\left(\hat{y}\right)$		$\hat{y}$	$v\left(\hat{y}\right)$			${\hat{y}}^{\prime }$	$v\left({\hat{y}}^{\prime }\right)$		${\hat{y}}^{\prime }$	$v\left({\hat{y}}^{\prime }\right)$
Typical NFI (6/2)	‒257, 912.4	1.72E + 10		0.18	‒0.16	94.6		‒321, 555.9	1.91E + 10		‒24.45	‒0.55	93.7
Disturbance (20/2)	‒896, 952.7	1.97E + 10		0.06	0.36	94.9		‒639, 875.3	1.75E + 10		28.70	0.27	49.9
Disturbance (40/3)	‒896, 777.0	1.99E + 10		0.08	1.30	94.9		‒365, 360.2	1.47E + 10		59.29	0.97	1.6
Disturbance (60/4)	‒898, 699.9	2.04E + 10		‒0.14	0.81	94.9		‒65, 387.3	1.09E + 10		92.71	1.70	0
Disturbance (80/4)	‒894, 251.2	2.17E + 10		0.36	0.03	94.7		255, 548.7	5.85E + 09		128.47	‒2.13	0

 | Show Table

DownLoad: CSV

For change estimation under a large disturbance scenario where nonresponse is caused by the disturbance, standard post-stratified estimation produced biases in estimates that were dependent on the amounts of nonresponse present. When 20% of disturbance plots were complete nonresponse and 2% of disturbance plots had partial nonresponse (denoted 20/2), the estimated change was 28.7% smaller than the known population total (Table 5). Under the 40/3 nonresponse scenario, the estimate was 59.3% too small, the 60/4 nonresponse rate produced an underestimation of 92.7%, and 80/4 nonresponse was 128.5% too small. In contrast, the differences between RHG estimates and known population totals were approximately 0.1% for nonresponse scenarios 20/2, 40/3, and 60/4, respectively, and nearly 0.4% for the 80/4 nonresponse (Table 5).

From a more practical perspective, users often rely on a confidence interval to provide an indication of estimate reliability. Typical interpretation of a 95% confidence interval would be that in repeated sampling of the population, 95% of the confidence intervals would contain the true population value. However, biased estimates can cause disruption to the theoretical construct due to the shift in confidence interval central location. As the estimates from the RHG approach for the disturbance population were largely unbiased, the confidence interval coverage of the true value was consistently near 95%. However, post-stratified estimates of the disturbance suffered poor coverage, with only about 50% coverage for the 20/2 nonresponse situation and near zero coverage for more intense nonresponse outcomes (Table 5).

4.   Discussion

The results for the NFI population indicated the RHG method produced estimates of change in the response variable that were essentially unbiased. However, there was a tendency for nonresponse to cause overestimation (more change was estimated than occurred in the population) of change in the NFI population when using post-stratified estimation. In this case, the preponderance of nonresponse on F/F plots that generally have positive change values shifts the results more towards NF/NF (zero change) and the group of status change (NF/F and F/NF) plots which had an overall tendency to negative change (Table 2). Other population change scenarios may provide different results. For example, if the population change was positively valued and status change plots had smaller change than the F/F plots then the change estimate would be too small (underestimation). Outcomes would also be dependent on the proportions of plot types within the RHG and on both the post-stratum and RHG weights. Temporal correlations in plot nonresponse may also interact with population change dynamics to affect estimation results. Thus, the population used here may be considered one of numerous potential NFI scenarios and inventory practitioners are encouraged to make evaluations in the context of their specific inventory circumstances. Nonetheless, potentially adverse effects of nonrandom nonresponse on population change estimates using the typical post-stratified estimator were demonstrated. This outcome is troubling as inaccurate estimates of change that are too large are more likely to be judged as being of consequential practical/statistical importance with subsequent implementation or continuation of activities and policies that may actually not be achieving desired outcomes. Similarly, underestimation of change may provide indications of status quo in the population and thus no responses to alter trajectories are needed (Field et al., 2004; Legg and Nagy, 2006).

Further, it may be worthwhile for those desiring to estimate disturbance impacts to consider the effects disturbance can have on nonresponse patterns and amounts (Nelson et al., 2009; Smith and Gray, 2021). As with the NFI population, estimates based on the RHG approach for the disturbance population were essentially unbiased at all levels of nonresponse studied (Table 5). In standard post-stratified estimation, differences in response probabilities between disturbed and undisturbed areas produced substantial underestimation of disturbance amounts. The results pertaining to the population studied here suggest that the disturbance effect is underestimated by nearly 30% under a 20/2 nonresponse scenario, with even further underestimation associated with higher levels of nonresponse. In extreme cases such as 80/4 nonresponse occurrence, the bias was severe enough to produce a positive estimate of change when the true population value was negative. Thus, catastrophic disturbances that cause increased nonresponse in the disturbed area can produce very misleading results under a typical implementation of post-stratified estimation. This suggests management and policy reactions to disturbance outcomes may be inadequate as the severity of impacts is underestimated in proportion to the amount of nonresponse. Failure to recognize the full extent of disturbance impacts may have long-term detrimental socio-economic consequences as well as inhibiting effective land management practices (De Grandpré et al., 2018; Montagné-Huck and Brunette, 2018).

The primary issue that affects the standard post-stratified estimation approach is a combination of post-stratum heterogeneity and the assumption that completely inaccessible sample units are reflective of the mean based on observed plots within the post-stratum. For example, in post-stratum 1 under a typical NFI change scenario (Table 2), there are nearly 90% of population elements that are NF–NF and therefore: 1) have a change value of zero, and 2) will never be subjected to nonresponse. Thus, the post-stratum total (or mean) change will tend to be relatively small due to the large number of zeros. In contrast, other sample elements in the post-stratum may have relatively large change values, e.g., F–NF, and also are potential nonresponse candidates. When the post-stratum mean is small but sample units with relatively large change values are nonresponse, it is incongruous that the post-stratum mean and nonresponse sample unit values will be of similar magnitude. For estimates of current status, it is often the case that post-stratum means are underestimated in the presence of nonresponse (Westfall, 2022). In the context of change estimation, the phenomenon is more complex due to the presence of both positive and negative values which may result in either overestimation or underestimation of post-stratum means depending on how nonresponse transpires in a given sample. Nonetheless, bias in the sample often arises from nonresponse and mitigation via RHG or other methods (Magnussen et al., 2017; Westfall and Wilson, 2022) can reduce negative impacts on population change estimates.

Estimated variances for both estimation approaches performed well with respect to the simulation variance, however the RHG approach generally had smaller differences than the standard post-stratified method. Astute readers may have noticed the estimated variance was larger using the post-stratified estimators in comparison to RHG estimation for the NFI population but was smaller in the analysis of the disturbance population. Examination of the results using the NFI population suggested that creation of post-stratum variances from RHG variances (Eq. (11)) produced smaller values than directly estimating post-stratum variances (Eq. (4)). Conceptually, this is intuitive as RHG essentially post-stratify the original post-strata into more homogenous assemblages – which is the underlying premise of variance reduction via (post-) stratified estimation (Cochran, 1977). This effect translated into overall variances being smaller for RHG estimation when combining variances across post-strata (Eq. (6)).

Conversely, for the disturbance population the means of the DIST and UNDIST RHG (Table 4) are substantially different than the overall post-stratum mean. When the post-stratum variance is calculated from RHG variances, the third term of Eq. (11) becomes relatively large and contributes to larger values than are realized using the standard post-stratified approach (Eq. (4)). In addition, the 4th term of Eq. (11) makes a larger contribution as the amount of nonresponse increases, e.g., the 4th term contributes an additional 8% to the total variance under the 80/4 nonresponse scenario. Another phenomenon that occurs with the disturbance population is that as nonresponse amounts increase in the disturbed zone, the sample becomes more homogeneous (mostly undisturbed plots) and this situation favors smaller variances using post-stratified estimators; whereas the RHG approach maintains the DIST/UNDIST distinction mentioned previously that produces higher variances. Särndal and Lundström (2005) provide a useful perspective on the situation by stating that ' … if an estimator is greatly biased, it is poor consolation that its variance is low'.

Within each post-stratum, various factors (e.g., ownership) can cause nonrandom nonresponse among sample units. The purpose of the RHG method is to further partition sample units such that the assumption of nonresponse occurring randomly and with equal probability is better satisfied. The simulations applied nonresponse under this assumption and thus robust agreement was obtained between population estimates and known values. It should be acknowledged that, while the RHGs used in this study are intuitive with regards to basic nonresponse tenets, in practice the construction of RHGs for nonresponse in forest inventories such as NFIs may include various and possibly interrelated factors that are difficult to identify. Thus, formulation of RHG where plots therein truly have (approximately) equal response probabilities may be a challenging endeavor. First, identification of the underlying nonresponse mechanisms may prove elusive without further focused investigation (Gao et al., 2020). For some NFIs, further refinements that account for different nonresponse probabilities between public and private land ownership or other known accessibility factors may be appropriate if not already included in development of the post-strata. Also, methodological differences in forest (F) and nonforest (NF) plot status delineation and plot vs. stratum definition (e.g., land use vs. land cover) may alter the distribution of plots requiring a field visit or introduce class mismatches between plots and their assigned RHG. Another factor not explicitly addressed within this study is RHG construction to minimize bias that may occur if nonresponse plots tend to cluster spatially within a certain areal subdomain of the population. Subsequently, a large number of combinations of relevant factors may preclude full accommodation into an RHG structure due to small sample sizes in certain assemblages and may result in corresponding sub-optimal remediation of nonresponse bias effects on population estimates.

The RHG methodology presented in the paper represents one potential way of addressing nonresponse bias in a forest inventory context. Other methods include ignoring missing values or employing different types of weighting schemes or imputation approaches, for which Schreuder et al. (1993) provide helpful coverage in forest inventory applications. Schreuder et al. (1993) also reference the Kalton and Kasprzyk (1986) publication which provides an excellent overview of weighting and imputation methods. In the context of this paper, an approach that considers weighting based on plot nonresponse probabilities may be considered. However, a remaining research need is to assess these different approaches in practical forest inventory settings. It is likely that optimal solutions may vary due to a number of factors, including the ability to mitigate nonresponse bias, effects on the standard errors of estimates, operational cost and feasibility, and complexity of implementation. Nonetheless, a broader understanding of the advantages and disadvantages of various missing data compensation methods would provide practitioners a framework for identification and implementation of a suitable solution pertinent to their specific circumstances.

5.   Conclusion

Accounting for potential nonresponse bias in estimates of change is more complex than for current status evaluations, primarily due to the need to assess nonresponse outcomes and underlying mechanisms at two points in time. In this study, nonresponse was based on plots requiring a field visit at either measurement time with commonly reported mechanisms of hazardous conditions or lack of access permission being assumed. Depending on the structure of change in the population and nonresponse severity, substantial differences between the post-stratified estimate and the true value can occur. For the forest inventory change scenarios presented, this study demonstrated that simply relying on post-stratum means to represent the average value of nonresponse plots is often an inadequate strategy for mitigating nonresponse bias. As an alternative, identification of nonresponse causal factors allows sample plots to be assigned to RHG that provide better representation of nonresponse plot values via group means from observed plots. In most cases, the plot information needed for placing plots into RHG based on factors such as office vs field plot and public vs private ownership already exists in the inventory database or can be otherwise easily obtained, e.g., via GIS map layers and plot location information. Accounting for more complex causal factors (e.g., ecological) may be considerably more difficult. In that context, nonresponse will remain nonrandom with respect to any underlying factors not specifically addressed in RGH formation. Nonetheless, the RHG change estimation framework described herein provides a potential pathway for forest inventory practitioners to reduce nonresponse effects on forest resource change estimates.

Author contributions

James A. Westfall: Conceptualization, Formal Analysis, Writing- Original draft preparation. Mark D. Nelson: Conceptualization, Writing- Reviewing and Editing.

Data availability

The data used in this study are available from the author upon request.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.