Hierarchical models for small area estimation using zero-inflated forest inventory variables: comparison and implementation

The motivating data are from the USDA Forest Service FIA program and comprise inventory plot measurements of live aboveground tree biomass density (Mg/ha). These data were drawn from the most current sampled measurement for each plot in the FIA database downloaded on February 8, 2023 for the states of Washington and Nevada (Burrill et al.,, 2023). Washington was selected because it has large differences in biomass across counties, ranging from massive biomass densities on the Olympic Peninsula to near zero biomass in counties east of the Cascade Mountain Range. Nevada was selected because it has a fairly unique distribution of forest biomass, where much of the state’s arid environment has little to no biomass which is punctuated with sky islands where there is non-zero forest biomass. Each state has approximately ten years of FIA data, ending in year 2019, that were derived from a panel of plots measured annually across a systematic sample of hexagons approximately 2,500 hectare in size. Biomass values were from live trees only and include all trees 1.0 inch diameter and greater.

Estimators and simulations, described in Sections 2.2 and 2.6, respectively, were informed using five auxiliary variables: National Land Cover Dataset Analytical Tree Canopy Cover 2016 (hereafter tcc); LANDFIRE 2010 Digital Elevation Model (hereafter elev); US Geological Survey Terrain Ruggedness Index (hereafter tri); PRISM mean annual precipitation, 30yr normals (1991-2020) (hereafter ppt); and LANDFIRE 2014 tree/non-tree lifeform mask (hereafter tnt) (Yang et al.,, 2018; U.S. Geological Survey,, 2019; Daly et al.,, 2002; Rollins,, 2009; Picotte et al.,, 2019). The tcc variable is a measure of average tree canopy cover in a given pixel, the elev variable gives the elevation at a given pixel, the tri variable gives the terrain ruggedness at a given pixel, the ppt variable is a measure of average precipitation at a given pixel over 30 years, and the tnt variable is a binary variable distinguishing between pixels with and without trees. These auxiliary variables were resampled to 90 meter resolution and available wall-to-wall in both states. At locations with FIA plots these variables are matched with the corresponding plot and then used as predictors in the estimators’ regression component and to inform simulated population generation. The use of these variables in different model types and simulated population generation is shown in Table 1.

Nine candidate model-based estimators are used to estimate the average biomass at the county level across Nevada and Washington. The first estimator follows a frequentist mode of inference and is constructed using a two-stage regression. The eight additional estimators use a Bayesian mode of inference and are constructed using one- and two-stage regressions. A brief description of the candidate estimators is provided in Table 2.

Models are fit at the unit level, which sets plot-level biomass (Mg/ha) as the response variable. To better meet subsequent models’ assumption of normally distributed residuals and to ensure positive support for predictions, models are fit to a transformed response variable. Often in such settings the logarithm is a natural choice; however, we found the logarithm to be too strong, resulting in a skewed response distribution once zeros were accommodated. We prefer root function transformations that were chosen to be state specific, with stronger roots for high biomass states and weaker roots for low to moderate biomass states. Specifically, we used a fourth and square root transformation for Washington and Nevada, respectively.

Despite our models being fit at the unit level, our inferential goal is to estimate average biomass per hectare at the county level, which we denote $\mu_{j}$ where $j$ indexes county within a given state. In subsequent sections we present the candidate unit-level models and how each is used to estimate the small area parameters of interest.

Here we develop a frequentist two-stage approach to estimation (F ZI CVI). This estimator and its MSE estimator were introduced by Chandra and Sud, (2012) and later applied to forest inventory data by White et al., (2025) with promising results in the state of Nevada. At generic spatial location $\ell$ the transformed non-zero forest biomass is modeled as

where $\beta_{0}$ is the intercept, $\tilde{\beta}_{0}\left(\ell\right)$ is the county specific random effect with $\tilde{\beta}_{0}\left(\ell\right)=\tilde{\beta}_{0,j}$ when $\ell$ is in the $j$th county and with $\tilde{\beta}_{0,j}\overset{\text{iid}}{\sim}\mathcal{N}(0,\sigma^{2}_{\tilde{\beta}_{0}})$, ${\boldsymbol{x}}(\ell)$ is a $p\times 1$ vector of predictor variables, ${\boldsymbol{\beta}}$ is a $p\times 1$ vector of regression coefficients, and $\varepsilon\left(\ell\right)$ is a residual error term with $\varepsilon\left(\ell\right)\overset{\text{iid}}{\sim}\mathcal{N}(0,\tau^{2})$.

Biomass presence and absence is modeled using a Bernoulli mixed model with a logit link function defined as

where $p(\ell)$ denotes the probability of non-zero response value at location $\ell$, $\alpha_{0}$ is the intercept, $\tilde{\alpha}_{0}\left(\ell\right)$ is the county specific random effect with $\tilde{\alpha}_{0}\left(\ell\right)=\tilde{\alpha}_{0,j}$ when $\ell$ is in the $j$th county and with $\tilde{\alpha}_{0,j}\overset{\text{iid}}{\sim}\mathcal{N}(0,\sigma^{2}_{\tilde{\alpha}_{0}})$, ${\boldsymbol{v}}(\ell)$ is a $q\times 1$ vector of predictor variables, ${\boldsymbol{\alpha}}$ is a $q\times 1$ vector of regression coefficients.

Here we consider a class of hierarchical Bayesian models. These models are primarily two-stage estimators similar to the frequentest two-stage estimator developed in Section 2.3. For comparison, we also include some simpler, single-stage, models that do not explicitly accommodate excess zero values in the response.

The two-stage hierarchical model used in Finley et al., (2011) and adopted here is

The first level of hierarchy is the model for $z(\ell)$. In our case, we only consider two options for estimation of $z(\ell)$: first, setting $z(\ell)$ to 1; and second, using a linear mixed model with logit link function. In the second case, for simplicity, we only use one model for estimation of $z(\ell)$ across all Bayesian estimators. Here, a realization of $z(\ell)$ indicates whether or not the location $\ell$ is predicted to have biomass (1) or not (0). Then, we pass the realization of $z(\ell)$ into the second level of hierarchy where it determines the expression of the mean $m(\ell)$ and the associated variance term. In the subsequent development, $m(\ell)$ is a given linear mixed model with county-varying intercept (CVI), county-varying coefficients (CVC), or space-varying intercept (SVI) components (Table 2). We implement and compare four different model forms for $m(\ell)$. The residual variance is estimated through a parameter $\tau^{2}_{1}$ when $z(\ell)$ is 1, and set to a small value via a constant, $\tau^{2}_{2}$, when $z(\ell)$ is 0 (see Finley et al.,, 2011, Section 3 for details). For some candidate models we allow for county-specific residual variance (CRV) terms via county-specific $\tau^{2}_{1}$s (Table 2).

We now turn to the particular models that use this hierarchical structure. The simplest model we consider is the county-varying intercept (B CVI) model, which is the Bayesian equivalent to Eq. (1) and is defined as

with parameter and hyperparameter distributions defined as follows, $\beta_{0}\sim\mathcal{N}(0,\sigma^{2}_{\beta_{0}})$, $\tilde{\beta}_{0}(\ell)$ is the county specific random effect, i.e., $\tilde{\beta}_{0}\left(\ell\right)=\tilde{\beta}_{0,j}$ and $\tilde{\beta}_{0,j}\left(\ell\right)\overset{\text{iid}}{\sim}\mathcal{N}(0,\sigma^{2}_{\tilde{\beta}_{0}})$ when $\ell$ is in the $j$th county, ${\boldsymbol{\beta}}\sim\mathcal{N}({\boldsymbol{0}},\sigma^{2}_{\beta}{\boldsymbol{I}})$ with ${\boldsymbol{I}}$ being the $p$-dimentional identify matrix, and $\varepsilon(\ell)\overset{\text{iid}}{\sim}\mathcal{N}(0,\tau^{2})$, $\sigma^{2}_{\tilde{\beta}_{0}}\sim\mathcal{IG}(a_{\sigma^{2}_{\tilde{\beta}_{0}}},b_{\sigma^{2}_{\tilde{\beta}_{0}}})$, and $\tau^{2}\sim\mathcal{IG}(a_{\tau^{2}},b_{\tau^{2}})$. All hyperparameters were set to induce non-informative prior distributions.

Next, we consider the county-varying coefficient (B CVC) model that allows regression coefficients to vary by county. This model is defined as

where $\tilde{\boldsymbol{\beta}}(\ell)=(\tilde{\beta}_{1,j},\tilde{\beta}_{2,j},\ldots,\tilde{\beta}_{p,j})^{\top}$ with $\tilde{\beta}_{k,j}\overset{\text{iid}}{\sim}\mathcal{N}(0,\sigma^{2}_{\tilde{\beta}_{k,j}})$ and $\sigma^{2}_{\tilde{\beta}_{k,j}}\sim\mathcal{IG}(a_{\sigma^{2}_{\tilde{\beta}_{k,j}}},b_{\sigma^{2}_{\tilde{\beta}_{k,j}}})$ for $k=1,2,\ldots,p$ when $\ell$ is in the $j$th county.

The B CVI and B CVC models defined above are used in our analyses both as a single-stage model (by setting $z(\ell)=1$ in the hierarchical model) and in the two-stage setting. We are most interested in these simpler, single-stage, models for comparison with two-stage models as laid out in Table 2.

Turning now to the two-stage models. Each two-stage model introduced in this section uses the same first stage model; however, the framework we defined in Eq. (3) lends itself to a variety of different first stage models. We consider a range of second stage models.

The first stage model used for all two-stage models is a generalized linear mixed model with Bernoulli response and a county-varying intercept defined as

where $p(\ell)$ denotes the probability of non-zero response, and parameter and hyperparameter distributions defined as follows, $\alpha_{0}\sim\mathcal{N}(0,\sigma^{2}_{\alpha_{0}})$, $\tilde{\alpha}_{0}(\ell)$ is the county specific random effect, i.e., $\tilde{\alpha}_{0}\left(\ell\right)=\tilde{\alpha}_{0,j}$ when $\ell$ is in the $j$th county, ${\boldsymbol{\alpha}}\sim\mathcal{N}({\boldsymbol{0}},\sigma^{2}_{\tilde{\alpha}_{0}}{\boldsymbol{I}})$, and $\sigma^{2}_{\tilde{\alpha}_{0}}\sim\mathcal{IG}(a_{\sigma^{2}_{\tilde{\alpha}_{0}}},b_{\sigma^{2}_{\tilde{\alpha}_{0}}})$. All hyperparameters were set to induce non-informative prior distributions.

Now, we can combine Eq. (6) with Eq. (4) and specify the first two-stage hierarchical model (B ZI CVI) as

where $\varepsilon_{1}\left(\ell\right)\overset{\text{iid}}{\sim}\mathcal{N}(0,\tau_{1}^{2})$ and $\varepsilon_{2}\left(\ell\right)\overset{\text{iid}}{\sim}\mathcal{N}(0,\tau_{2}^{2})$ with $\tau_{1}^{2}\sim\mathcal{IG}(a_{\tau_{1}^{2}},b_{\tau_{1}^{2}})$ and $\tau^{2}_{2}=0.00001$.

The hierarchical model specified in Eq. (7) is analogous to the frequentist two-stage model specified in Section 2.3. Notably, this hierarchical structure combines the Bernoulli model with the B CVI model to produce estimates that account for zero inflation.

We next extend Eq. (7) to include county-varying coefficients (B ZI CVC) defined as

In practice, it is common for residual variance to increase with increasing biomass, i.e., heteroscedasticity. Given disparity in forest density at the county-level, we might expect residual variance to vary across counties. For example, western counties in Washington have much more forest biomass than counties in eastern Washington, and if this disparity in biomass is not entirely captured by the predictors and random effects, then we would expect quite different residual variances. A county-specific residual variance term can help accommodate heteroscedasticity and improve county-level estimates; hence, we extend the B ZI CVI Eq. (7) and B ZI CVC Eq. (8) models with county-specific residual variance parameters. Specifically, the zero-inflated county-varying intercept model with county-specific residual variance (B ZI CVI CRV) and the corresponding county-varying coefficient model (B ZI CVC CRV) are defined analogous to Eq. (7) and Eq. (8) but with $\varepsilon_{1}\left(\ell\right)\overset{\text{iid}}{\sim}\mathcal{N}(0,\tau_{1,j}^{2})$ with $\tau_{1,j}^{2}\sim\mathcal{IG}(a_{\tau_{1}^{2}},b_{\tau_{1}^{2}})$ when $\ell$ is in the $j$th county.

In addition to county scale differences in mean biomass and predictor variables’ relationships with biomass, captured through county-varying intercepts and county-varying coefficients, respectively, we might expect to see smoothly-varying spatially structured changes in mean biomass caused by disturbance history, climate impacts, species composition, or any spatially dependent factors not captured by predictors variables. Such spatial changes in mean biomass can be accommodated via a space-varying intercept random effect. The models below extend the B ZI CVI CRV and B ZI CVC CRV models to include such a space-varying intercept. Specifically, the B ZI CVI SVI CRV model is defined as

where $w(\ell)$ is a spatial random effect that adjusts the intercept based on residual spatial dependence. Here we estimate $w(\ell)$ using a GP approximation called the Nearest Neighbor Gaussian Process (NNGP; Datta et al., 2016; Finley et al., 2019) that provides substantial improvements in run time, with negligible differences in inference and prediction, compared to a model that uses a full Gaussian Process. In brief, for this specification, the vector of random effects collected over $n$ locations ${\boldsymbol{w}}=\left(w(\ell_{1}),w(\ell_{2}),\ldots,w(\ell_{n})\right)^{\top}$ is distributed multivariate normal with mean zero and covariance matrix that captures the spatial dependence among random effects, i.e., ${\boldsymbol{w}}\sim\mathcal{MVN}({\boldsymbol{0}},\sigma^{2}_{w}{\boldsymbol{R}}(\phi))$, where $\sigma^{2}_{w}$ is the spatial variance and ${\boldsymbol{R}}(\phi)$ is the NNGP-derived correlation matrix that depends on a spatial correlation function, which in our case is an exponential, and decay parameter $\phi$ used to estimate the strength of correlation between any two locations. As with other variance parameters, we assume $\sigma^{2}_{w}\sim\mathcal{IG}(a_{\sigma^{2}_{w}},b_{\sigma^{2}_{w}})$ with hyperparameters set to induce a non-informative prior distribution. The spatial decay parameter is assumed to follow a uniform distribution with broad, non-informative, spatial support.

The last candidate model ZI CVC SVI CRV extends B ZI CVI SVI CRV Eq. (9) to include county-varying coefficients.

A simulation study was used to assess qualities of the estimators introduced in Section 2.2. Simulation studies are particularly useful in SAE research as they allow for us to assess estimators against simulated true parameter values, letting us gain accurate insight into estimator performance. In order to assess estimators in a manner that is fair and realistic, we utilize methodology introduced in White et al., 2024a for our simulation study. Briefly, the approach uses a $k$ nearest neighbors ($k$NN) algorithm on auxiliary data weighted with bootstrap inclusion probabilities to impute forest inventory attributes at each population unit. Formally, the algorithm to generate the population for a given state is defined in White et al., 2024a Algorithm 1.

It is important to note that we generate simulated populations separately in Washington and Nevada. Further, as discussed in Section 2.1, we use tcc, elev, and ppt as auxiliary data in the $k$NN matching in Washington, and tcc, elev, and tri in Nevada. We also use tnt as a stratification variable for generating the simulated population. These auxiliary variables are centered and scaled before the matching occurs. Without the centering and scaling step, variables of different magnitudes would hold different weight in the $k$NN search. We implement the generation of the simulated population via the kbaabb R package (White et al., 2024b, ).

After generating the simulated population we take design-based samples from the simulated population for the purposes of estimator assessment. In order to create one sample from the simulated population, we take a simple random sample from each county of the same size as the number of FIA plots in that county. Therefore, both the county and overall sample sizes remain constant.

We fit and evaluated the nine estimators introduced in Section 2 across all samples taken from the simulated population generated in Section 2.6. These estimators are fit separately between the two states considered in this study, Nevada and Washington. We first turn to evaluating performance metrics in Washington, and then to Nevada. Generally, estimator performance differs between the two states, likely due to the substantial differences in the ecological landscape of the states, with Nevada’s forests primarily occupying high elevation sky islands, and Washington’s forests primarily on the west side of the Cascade mountain range.

Figure 1 displays values for the four performance metrics evaluated in Nevada, with each point on the plot represented a county-performance metric-estimator combination. Turning our attention to Figure 1a, we see that the single stage estimators exhibit the most bias. Further, we see that estimators that do not allow for county-specific variances tend to exhibit more bias than those that allow for county-specific variances. The least bias estimators are the B ZI CVI CRV and the B ZI CVI SVI CRV, both performing quite similar to each other. The CVC models, when compared to their CVI counterparts, tend to exhibit more bias after CRV and/or SVI terms are added. When turning to Figure 1b, the negative effect of the CVC effects in Nevada becomes even clearer, as we see these models exhibit much higher RMSE than those without the CVC terms. All two-stage estimators with the CVI effect perform comparably in terms of RMSE, with the others performing comparably to each other. Figure 1c displays the bias of $\widehat{\mbox{RMSE}}$, indicating if an estimator is over- or under-estimating its variability. The B CVI and B CVC tend to substantially under estimate their variability, likely due to the poorly specified underlying model. Somewhat surprisingly, the B ZI CVI estimator also exhibits some negative bias. The F ZI CVI does quite a good job of estimating its RMSE, but has more variability in those estimates than the B ZI CVC, B ZI CVI CRV, B ZI CVC CRV, B ZI CVI SVI CRV, and B ZI CVC SVI CRV, all of which estimate the RMSE quite well, but have a few negative outliers. Figure 1d displays the 95% uncertainty interval coverage, where we can see the B CVI and B CVC estimators performing very poorly. The F ZI CVI performs better than its Bayesian counterpart in this case, but both have some counties with very low coverage rates. The remaining estimators all perform quite well and similar to each other, with the B ZI CVI SVI CRV performing the best on median.

Figure 2 displays values for the four performance metrics evaluated in Washington, with each point on the plot represented a county-performance metric-estimator combination. Figure 2a displays the bias of each estimator, and we see in Washington that the F ZI CVI estimator exhibits the most bias, substantially more than its Bayesian counterpart. Initially this may come as a surprise due to the estimators similar structure, but one must consider the differences in the back-transformation of the response variable between these estimators. While the frequentist estimator back-transforms the predictions from each model, the Bayesian estimators back-transforms each sample from the posterior predictive distribution, potentially leading to more sensible estimates under certain conditions. The B CVI and B CVC estimators tend to exhibit greater magnitude bias than the two-stage Bayesian estimators, all of which exhibit similar amounts of bias to each other. One notable observation is the B ZI CVC CRV and B ZI CVC SVI CRV estimators exhibit a slight amount more bias than their CVI counterparts, and when we turn to Figure 2b we can see their RMSE is higher than their CVI counterparts. The B ZI CVI CRV and B ZI CVI SVI CRV estimators have the lowest RMSE, while the F ZI CVI estimator has the highest RMSE. Other two-stage Bayesian estimators exhibit very similar RMSE to each other, and the single-stage estimators exhibit a bit higher RMSE. Turning to Figure 2c, we see the bias of the RMSE estimator. Notably, the F ZI CVI and B CVI estimators exhibit the most bias here, while the other estimators exhibit very low bias. Once CRVs are accounted for, we see almost no bias. Finally, turning to Figure 2d, we examine the 95% uncertainty interval coverage rate. The B CVI and B CVC estimators have very poor coverage, while the other estimators, all of which are two-stage, have decent coverage among themselves. Notably, the estimators which account for CRVs and have CVC effects exhibit the best coverage rates.

Stepping back from the details of Figures 1 and 2, we can gain some broader insights about the estimators and the simulation study. First of all, it is striking that the estimators that include a SVI perform almost identically to their non-spatial counterparts. This is not to say that spatial effects are unuseful for estimating biomass or other forest attributes, in fact we know from other studies these effects are quite useful (Finley et al.,, 2024). Further, from the cross-validation carried out in Section 3.2 we see more substantial improvements with the spatial models. It is the case though, that the simulated population we generated does not necessarily reflect the spatial structure of the true population, and in fact upon inspection it has very little spatial structure. Incorporating some spatial smoothing into the population generation methodology would help improve the simulated population’s utility for assessing these spatial models. Another broad observation we can make from the simulation study is that the CVCs generally did not improve estimation, and in some cases they worsened estimation. We see this particularly in Nevada, where county-level biomass is more homogeneous across the state. In Washington, the estimators with CVCs did have better 95% coverage rates than their CVI counterparts, but at the cost of higher bias and RMSE. The use of CVCs may be more useful for models fit across larger spatial domains (such as the entire United States), where we would expect the effect of predictors on the response to vary substantially. States with more diverse forest types may also benefit from the use of CVCs.

We applied the estimators discussed to actual plot data collected by FIA in both Nevada and Washington. We first compare the estimators using 10-fold cross-validation, and then turn to discussing estimates produced by the B ZI CVI SVI CRV estimator, an estimator that showed favorability in the simulation study and in cross-validation.

We performed 10-fold cross-validation at the unit-level to help assess the estimators performance, with the unit-level predictions as a proxy for estimation of county means. Results from the cross-validation are shown in Table 3. To compute 95% uncertainty interval coverage rates (“Coverage” in Table 3), we used the quantiles of the posterior predictive distribution of $y(\ell)$ for the Bayesian estimators. We omit 95% confidence interval coverage for the frequentist estimators as the task of computing closed form or bootstrap prediction intervals at the unit-level for this two-stage model is beyond the scope of this article. Regarding the results, we can see similar patterns as we saw in Figures 1 and 2, but now we can better see the performance improvement associated with the spatial models when considering the root mean square prediction error (RMSPE) and bias. We see that across states bias and coverage rates were favorable for all estimators that included a county-specific variance term, and in Nevada for all estimators that account for zero-inflation. Examining the results of this cross-validation is helpful for model evaluation and in this case provides similar insights to the simulation study, and in Section 4 we discuss the tradeoffs of the two model evaluation methods.

In both Nevada and Washington, the B ZI CVI SVI CRV estimator performs quite well when considering all computed metrics in the simulation study and cross-validation, so we chose to use this estimator as an example for producing pixel-level predictions and county-level estimates in Washington and Nevada. Figure 3a, b show pixel-level estimates of biomass probability and Figure 3c, d show estimated biomass (Mg/ha) in both Washington and Nevada produced by the B ZI CVI SVI CRV estimator. While we are not particularly interested in the pixel-level predictions for the study at hand, these sorts of maps can be very useful for management at very small spatial scales, such as the stand level. In this Bayesian framework, not only do we have pixel-level predictions, but we also have pixel-level uncertainty predictions. Further, we can easily summarize these pixel level predictions to get estimates and uncertainty estimates for any area that may be of interest. For the study at hand, we only produce estimates for counties, but we could have just as easily produced estimates for any custom small area of interest (e.g., ecoregions, watersheds, fires, stands, etc.) with any of the Bayesian estimators.

Figure 4 shows the estimated biomass (Mg/ha) in both Nevada and Washington produced by the B ZI CVI SVI CRV estimator. These sorts of estimates are what we are most interested for SAE, and for demonstration we only produce estimates at the county-level. County-level estimates of biomass show the stark change in biomass across the cascades in Washington state, and the relatively constant county-level biomass across the state of Nevada.

The authors thank the Forest Inventory & Analysis Program for the data.

This work was supported by the USDA Forest Service Forest Inventory and Analysis Program, Region 9, Forest Health Protection, Michigan State University AgBioResearch, National Science Foundation DEB-2213565 and DEB-1946007, and National Aeronautics and Space Administration CMS Hayes-2023. The findings and conclusions in this publication are those of the authors and should not be construed to represent any official US Department of Agriculture or US Government determination or policy.

The data include confidential plot data, which can not be shared publicly. FIA data can be accessed through the FIA DataMart (https://apps.fs.usda.gov/fia/datamart/datamart.html). Requests for data used here or other requests including confidential data should be directed to FIA’s Spatial Data Services (https://www.fs.usda.gov/research/programs/fia/sds).