Models to support forest inventory and small area estimation using sparsely
sampled LiDAR: A case study involving G-LiHT LiDAR in Tanana, Alaska

Andrew O. Finley1,2, Hans-Erik Andersen3, Chad Babcock4, Bruce D. Cook5,

Douglas C. Morton5, and Sudipto Banerjee6

1Department of Forestry, Michigan State University, East Lansing, MI, USA.
2Department of Statistics and Probability, Michigan State University, East
Lansing, MI, USA.
3USDA Forest Service, Pacific Northwest Research Station, Seattle, WA, USA.
4Department of Forest Resources, University of Minnesota, St Paul, MN, USA
5National Aeronautics and Space Administration, Goddard Space Flight Center,
Greenbelt, MD, USA.
6Department of Biostatistics, University of California, Los Angeles, CA, USA.

Corresponding Author: Andrew O. Finley<EMAIL_ADDRESS>
## Abstract

A two-stage hierarchical Bayesian model is developed and implemented to
estimate forest biomass density and total given sparsely sampled LiDAR and
georeferenced forest inventory plot measurements. The model is motivated by
the United States Department of Agriculture (USDA) Forest Service Forest
Inventory and Analysis (FIA) objective to provide biomass estimates for the
remote Tanana Inventory Unit (TIU) in interior Alaska. The proposed model
yields stratum-level biomass estimates for arbitrarily sized areas. Model-
based estimates are compared with the TIU FIA design-based post-stratified
estimates. Model-based small area estimates (SAEs) for two experimental
forests within the TIU are compared with each forest’s design-based estimates
generated using a dense network of independent inventory plots. Model
parameter estimates and biomass predictions are informed using FIA plot
measurements, LiDAR data that are spatially aligned with a subset of the FIA
plots, and complete coverage remotely detected data used to define
landuse/landcover stratum and percent forest canopy cover. Results support a
model-based approach to estimating forest variables when inventory data are
sparse or resources limit collection of enough data to achieve desired
accuracy and precision using design-based methods.

##  1 Introduction

Large scale forest monitoring programs have traditionally used design-based
inference that uses probability sampling and associated estimators to deliver
forest variable estimates. For example, the United States Department of
Agriculture (USDA) Forest Service Forest Inventory and Analysis (FIA) program
conducts the US national forest inventory (NFI), collecting data describing
the condition of forest ecosystems on a large network of permanent inventory
plots distributed across all lands in the nation (Smith,, 2002). These data
offer a unique and powerful resource for determining the extent, magnitude,
and causes of long-term changes in forest health, timber resources, and forest
landowner characteristics across large regions in the US (Wurtzebach et al.,,
2020). The FIA program uses design-based post-stratified estimators to improve
precision of point and change estimates (Westfall et al.,, 2011; Bechtold and
Patterson,, 2005). Depending on the desired level of estimate precision, such
approaches often require costly measurements over a relatively dense network
of inventory plots; hence, from a cost efficiency standpoint, there is
interest in methods that can deliver comparable inference using fewer
inventory plots. At the same time, like other NFIs (Breidenbach and Astrup,,
2012; Köhl et al.,, 2006), FIA has experienced increased demand for estimates
within smaller spatial, temporal, and biophysical extents than design-based
inference can reasonably deliver (e.g., annual or stand-level estimates).
Developing estimation methods that support inference on small areas—referred
to as small area estimation (SAE) methods—using FIA data is an active area of
research, with considerable progress made in the last several years (Hou et
al.,, 2021; Coulston et al.,, 2021; Schroeder et al.,, 2014; Lister et al.,,
2020). SAE methods are numerous and diverse, though most seek to improve
inference on small areas by making use of statistical models and auxiliary
information that is correlated with outcome variables (Rao and Molina,, 2015).

The FIA and similar NFIs continue to explore the efficacy of model-assisted
and model-based modes of inference to reduce cost and improve precision for
both large and small area estimates. These approaches often leverage rich
information content of satellite and airborne remote sensing to augment
information gleaned from the inventory plot network. Model-assisted approaches
employing complete coverage auxiliary (e.g., remote sensing) information to
improve the precision of inventory estimates within the design-based
inferential paradigm have been presented (Breidt and Opsomer,, 2017;
McConville et al.,, 2017; Strunk et al.,, 2019; Magnussen et al.,, 2018;
Ekström and Nilsson,, 2021), while model-based techniques have been developed
to provide estimates in cases where useful models are available to relate
remote sensing metrics to inventory measurements (Ståhl et al., 2011a, ;
McRoberts,, 2010; Finley et al.,, 2011; Saarela et al.,, 2016; Babcock et
al.,, 2018; May et al.,, 2023).

Model-based SAE methods offer a valuable alternative to the design-based post-
stratified estimators implemented by FIA. Model-based SAE methods seek to
borrow information from outside the area of interest (e.g., from neighboring
areal-units or point-referenced observations) and auxiliary data (e.g., remote
sensing data) to improve precision of estimated quantities. SAE methods can
generally be classified into two groups: unit-level and area-level models.
Unit-level models are constructed at the level of population units, where a
population unit is defined as the minimal unit that can be sampled from a
population. With respect to FIA’s survey design, field plot centers represent
population units. Unit-level models typically relate outcome variable
measurements on sampled population units to auxiliary data that is available
for all population units (e.g., complete coverage remote sensing data).
Prediction for a small area is achieved by aggregating unit-level predictions
within the given areal extent (Rao and Molina,, 2015). In contrast, area-level
models are constructed across areal units where relationships are built
between area-specific outcome variable direct estimates (e.g., generated using
design-based estimators) and auxiliary data (Rao and Molina,, 2015). Hence,
area-level models effectively “adjust” direct estimates given auxiliary
information.

The models developed and implemented in this paper are motivated by FIA’s
objective to improve biomass estimates in interior Alaska. FIA plot
measurements in interior Alaska’s Tanana Inventory Unit (TIU) were carried out
over four years starting in 2014. Implementing FIA’s standard sampling
intensity is cost-prohibitive in this remote region, due to challenging
logistics and high transportation costs—lack of roads requires that virtually
every plot is accessed via helicopter. For this reason, FIA has implemented a
modified sampling design in interior Alaska using a reduced sampling intensity
for field plots (1 plot per 12 140 ha), supplemented with high-resolution
airborne imagery acquired by NASA Goddard’s Lidar, Hyperspectral, and Thermal
(G-LiHT) Airborne Imager (Cook et al.,, 2013) in a sampling mode (Cahoon and
Baer,, 2022). Given the large geographic expanse of interior Alaska (approx.
46 million ha of forestland), a strip sampling mode is considered to be the
most economical and spatially balanced acquisition strategy for the airborne
data. The use of high-resolution airborne LiDAR strip sampling to support
forest inventories has been the focus of several recent studies in Europe and
North America where a variety of estimation approaches have been developed and
evaluated, including model-assisted (Andersen et al., 2009a, ; Gregoire et
al.,, 2011; Strunk et al.,, 2019), model-based (Saarela et al.,, 2018; Ståhl
et al., 2011b, ), and hybrid estimation (Ståhl et al.,, 2016). A comprehensive
account of model-based inference for survey data and its richness over design-
based estimators is offered by Little, (2004), which includes specific
pitfalls of the Horwitz-Thompson estimator. Recent explorations into Bayesian
survey sampling from spatially correlated frameworks include Chan-Golston et
al., (2020) and Chan-Golston et al., (2022) who show the benefits of modeling
finite populations as realizations of a spatial process with the latter also
accommodating modeling spatial associations in sampling indicators (also see
Finley et al.,, 2011).

Here, we apply a model-based approach to estimate forest biomass within the
TIU. The devised model yields stratum specific biomass estimates for
arbitrarily sized areas. We generate biomass estimates for the entire TIU and
two small areas. Model-based estimates are compared with the TIU FIA design-
based post-stratified estimates. Model-based SAEs for two experimental forests
within the TIU are compared with each forest’s design-based estimates
generated using a dense network of independent inventory plots (i.e.,
independent, meaning the data were not used to inform the model-based
estimates). Model parameter estimates and biomass predictions are informed
using FIA plot measurements, LiDAR data that is spatially aligned with a
subset of the FIA plots, and complete coverage remotely sensed data used to
define four landuse/landcover stratum and percent forest canopy cover. Model
data input and subsequent biomass predictions are point-referenced (i.e.,
indexed by spatial locations) hence, from a SAE standpoint, we pursue a unit-
level approach.

The remainder of the paper is as follows. Section 2 provides a description of
the data and exploratory analysis. The modeling framework and assessment
criteria are described in Section 3. Design-based estimators used to generate
estimates to compare with model-based estimates are provided in Section 4. A
description of the TIU and small area analyses are provided in Section 5.
Model selection results and biomass estimates for the TIU and small areas are
presented in Section 6. Results and description of possible next steps are
discussed in Section 7.

##  2 Data

We describe TIU data used to inform the proposed biomass models presented in
Section 3. The 13.533 million (ha) TIU is shown in Figure 1a along with data
locations and extents.

(a)

(b)

(c)

(d)

Figure 1: LABEL:sub@fig:tiu-data Tanana Inventory Unit with data collection
locations (FIA locations are perturbed) and two small areas of interest (SAI).
LABEL:sub@fig:tiu-data-zoom Zoom-in that shows all three location types (FIA
locations are perturbed). The dense grid of prediction locations exist across
the entire TIU. LABEL:sub@fig:tiu-strata Strata used in subsequent modeling.
LABEL:sub@fig:tiu-hansen Percent forest cover.

###  2.1 FIA plots

The standard FIA plot comprises four 1/60-th (ha) fixed-area, circular
subplots spaced 36.6 (m) apart (see Cahoon and Baer,, 2022, for details on FIA
plot design and measurements). Subplot coordinates were obtained using a
GLONASS-enabled Trimble GeoXH mapping-grade GNSS receivers that provide $<$ 2
(m) geolocation error (McGaughey et al.,, 2017; Andersen et al., 2009b, ;
Andersen et al.,, 2022). Following FIA procedures, individual tree dry biomass
was estimated, summarized to the plot-level, and expressed on a per ha basis
(see FIA DRYBIOT variable in Burrill et al.,, 2021). This plot-level
summarized and expanded DRYBIOT variable is $y$ (Mg/ha) in subsequent model
development.

As noted previously, TIU FIA plot measurements were collected over four years
(2014, 2016-2018). Over this period, 1 091 FIA plots were sampled. The subset
of 880 plots that spatially align with G-LiHT data (described in the Section
2.2) are depicted in Figure 1a. This subset of plots is used in subsequent
model development, whereas all sampled plots are used for subsequent design-
based estimates. To protect the plot integrity and ownership information, FIA
plot location is proprietary. While actual locations were used for modeling,
figures presented here depict spatially perturbed (i.e., fuzzed) locations.

###  2.2 G-LiHT LiDAR

To augment sparsely sampled FIA plots, linear swaths of high-resolution
airborne remote sensing measurements, placed approximately 9.2 (km) apart and
spatially aligned with most of the FIA plots, were acquired in 2014 using
G-LiHT mounted on a fixed-wing aircraft platform (Cook et al.,, 2013). G-LiHT
data specifications are provided in Table 1.

Specification |  | Specification |
---|---|---|---
Instrument | Riegl VQ-480 | Footprint size | 10 cm
Laser wavelength | 1550 nm | Half-scanning angle | 15 degrees
Flying height | 3 35 m (AGL) | Average pulse density | 3 pulses/m2
Beam divergence | 0.3 mrad | Swath width | 400 m
Table 1: G-LiHT LiDAR specifications for TIU.

A 1-by-1 (m) canopy height model (CHM) was computed as the difference between
G-LiHT derived digital terrain model (DTM) and a canopy surface model (CSM).
The average canopy height (m) was computed from the CHM over the 880 spatially
coinciding FIA plot footprints and a series of 61 029 “LiDAR plots” each with
area equal to the FIA plot and spaced 200 (m) apart along the linear flight
lights. The average CHM height variable is $x_{CH}$ in subsequent analyses.

G-LiHT LiDAR plot locations are shown in Figures 1a and 1b. We could have
tessellated the G-LiHT swath into 10s of millions of FIA plot sized LiDAR
plots; however, results from Finley et al., (2019), Shirota et al., (2022),
and Peruzzi et al., (2021) show there is strong spatial dependence among
G-LiHT canopy height metrics in TIU and hence little additional information
gain at the expense of massive model fitting computational cost.

###  2.3 Stratification and forest canopy cover

A complete coverage (i.e., covering the entire TIU) 30-by-30 (m) resolution
stratification variable was formed using forest and non-forest National Land
Cover Database (NLCD) (Homer et al.,, 2015) with stratum “Deciduous” (Class
41), “Conifer” (Class 42), “Mixed” (Class 43), and “Other” (all non-forest
classes) (Figure 1c). Using a model-assisted approach to analyze the TIU data
described above, Andersen et al., (2024) found this NLCD stratifying variable
explained a substantial portion of variability in biomass. Hence, we too use
this stratification variable to facilitate comparison with Andersen et al.,
(2024).

In addition to the NLCD stratification variable, a complete coverage 30-by-30
(m) resolution percent tree cover variable was formed using Hansen et al.,
(2013) fractional tree cover, updated to account for forest cover loss prior
to 2014 (Figure 1d). This tree cover variable is $v_{TC}$ in subsequent
analyses.

###  2.4 Exploratory data analysis

In Section 3, we define a hierarchical Bayesian regression model used to
predict biomass ($y$) using NLCD strata and G-LiHT $x_{CH}$. As with all
modeling efforts, we begin with exploratory data analysis (EDA) to build
intuition about the data and relationships of interest.

(a)

(b)

Figure 2: LABEL:sub@fig:boxplots-bio distribution of biomass and
LABEL:sub@fig:boxplots-rs G-LiHT mean canopy height $x_{CH}$ by NLCD strata.
Horizontal boxplot lines indicate the distribution’s quantiles and points
identify extreme measurements.

Figure 2a shows the distribution of FIA plot biomass density measurements is
different for each stratum. Conifer forests have lower median biomass than
Deciduous or Mixed forest types, but also several plots with very high biomass
(i.e., right skewed distribution). Among the forested strata, Deciduous has
the largest median biomass. As the stratum name suggests, biomass distribution
within Mixed appears to be a mix of Conifer and Deciduous distributions. As
described in Section 2.3, the Other stratum is a catch-all for non-forest
classes. Across the TIU, non-forest landcover is dominated by barren,
herbaceous, and wetlands, hence the forest biomass distribution in Other is
concentrated just above zero. The substantial right skew in Other is due to
forested FIA plots that fall within non-forest NLCD classes. Figure 2b shows
distribution of mean forest canopy height, summarized using $x_{CH}$, is also
different for each stratum. As observed in numerous studies, forest biomass
and forest canopy structure metrics, like $x_{CH}$, have a positive
correlation, therefore it is not surprising the patterns seen in Figure 2a are
reflected in Figure 2b.

Next we explore the relationship between biomass and mean forest canopy height
by stratum using the simple linear regression

$y=\beta_{0}+\beta_{CH}x_{CH}+\epsilon,$ (1)

where $\beta_{0}$ and $\beta_{CH}$ are the intercept and slope coefficient,
respectively, and $\epsilon\overset{\mathrm{iid}}{\sim}N(0,\tau^{2})$. This
model was fit to stratum specific data. These data and resulting regression
lines are shown in Figure 3. Sample size, $n$, and parameter estimates are
given in Table 2. Results suggest a strong linear relationship between biomass
and $x_{CH}$ within each stratum. Among the forested strata, the regression
slope coefficient for Conifer is larger than Deciduous and Mixed. The slope
coefficient for Other is substantially smaller than forested strata
coefficients. These differences suggest subsequent model development could
benefit from stratum-varying coefficients.

Figure 3: Biomass and G-LiHT derived mean canopy height $x_{CH}$ observed at
FIA plot locations. Model (1) regression lines correspond to parameter
estimates given in Table 2.

The last row in Table 2 shows parameter estimates for the pooled model (i.e.,
not broken out by stratum). Given stratum specific differences in $\tau^{2}$
estimates and magnitude of the pooled model’s $\tau^{2}$ estimate compared
with that of the Other’s model (i.e., $\sim$251 vs. $\sim$85), we might expect
subsequent model development could benefit from stratum specific residual
variance parameters.

|  | $\beta_{0}$ | $\beta_{CH}$ | $\tau^{2}$
---|---|---|---|---
Stratum | $n$ | ${}_{\text{(L. 95\%)}}\;\;\text{Mean}\;\;_{\text{(U. 95\%)}}$ | ${}_{\text{(L. 95\%)}}\;\;\text{Mean}\;\;_{\text{(U. 95\%)}}$ | ${}_{\text{(L. 95\%)}}\;\;\text{Mean}\;\;_{\text{(U. 95\%)}}$
Conifer | 265 | ${}_{(\text{-}5.206)}\;\;\text{-}1.616\;\;_{(1.986)}$ | ${}_{(12.059)}\;\;13.008\;\;_{(13.923)}$ | ${}_{(288.871)}\;\;341.878\;\;_{(403.191)}$
Deciduous | 29 | ${}_{(\text{-}15.817)}\;\;\text{-}4.501\;\;_{(6.070)}$ | ${}_{(7.973)}\;\;9.141\;\;_{(10.333)}$ | ${}_{(283.412)}\;\;413.825\;\;_{(609.904)}$
Mixed | 56 | ${}_{(\text{-}12.645)}\;\;\text{-}2.328\;\;_{(7.459)}$ | ${}_{(7.619)}\;\;8.952\;\;_{(10.435)}$ | ${}_{(141.354)}\;\;226.322\;\;_{(396.705)}$
Other | 530 | ${}_{(\text{-}2.765)}\;\;\text{-}1.906\;\;_{(\text{-}0.948)}$ | ${}_{(7.161)}\;\;7.596\;\;_{(8.022)}$ | ${}_{(76.062)}\;\;85.039\;\;_{(96.779)}$
Pooled | 880 | ${}_{(\text{-}2.464)}\;\;\text{-}1.166\;\;_{(0.132)}$ | ${}_{(9.518)}\;\;9.866\;\;_{(10.212)}$ | ${}_{(227.651)}\;\;251.259\;\;_{(275.621)}$
Table 2: Sample size $n$ and parameter estimates for stratum specific model
(1) fit using biomass and G-LiHT derived mean canopy height $x_{CH}$ observed
at FIA plot locations. Data and regression lines are shown in Figure 3.

Regression model diagnostic analyses identified a few potentially high
leverage observations and moderate residual heteroskedasticity. Additionally,
semivariogram analysis suggested negligible spatial dependence among model
residuals (i.e., spatial dependence among biomass measurements was, for the
most part, captured using $x_{CH}$). Previous TIU modeling efforts conducted
using FIA subplots, which allow for better estimate of the semivariogram
nugget due to closer spatial proximity, and denser plot networks (see, e.g.,
Babcock et al.,, 2018; Taylor-Rodriguez et al.,, 2019; Shirota et al.,, 2022),
showed that even after accounting for LiDAR derived forest canopy variables,
residual spatial dependence in biomass was present, and model fit and
predictive performance was improved via addition of spatial random effects.

##  3 Models

We propose a model to predict biomass density at any location(s) within the
TIU. The model accommodates stratum-varying regression intercept, slope, and
variance parameters, as well as residual spatial dependence. Biomass
prediction for an area (e.g., entire TIU or small area within the TIU) is
approximated by summarizing an appropriately dense grid of unit-level
predictions within the area. For example, biomass density and total estimates
for the TIU will be based on posterior predictive distributions estimated at
each location within a 250-by-250 (m) grid that extends across the TIU—a
portion of which is illustrated in Figure 1b. In our setting, a key issue with
this approach is that biomass is to be conditioned on G-LiHT derived canopy
structure metrics, e.g., $x_{CH}$, and these metrics are observed only along
the flight lines. This misalignment is remedied using a two-stage approach.
The first stage is a biomass process model that is conditioned on canopy
structure metrics, and the second stage comprises a separate process model for
each canopy structure metric. For the first stage model, we assume all biomass
measurements have spatially coinciding canopy structure measurements. When
cast within a hierarchical Bayesian framework, this two-stage approach allows
uncertainty in canopy structure metric predictions to be propagated through to
biomass predictions.

Let us consider a first stage biomass model at generic location $\ell$ within
the TIU. We model biomass $y(\ell)$ using a set of $p$ canopy structure
variables ${\bm{x}}(\ell)$ measured using G-LiHT. Based on the exploratory
analysis, we allow the relationship between these canopy structure variables
and biomass to vary by stratum via $p$ stratum-varying coefficients
$\tilde{{\bm{\beta}}}_{j}$ for $j$ in $1,2,\ldots,q$, where $q$ is the number
of strata. Additionally, an underlying latent spatial process $w(\ell)$
accounts for local changes in biomass attributed to unobserved, but smoothly
varying, environmental influences. Note that a spatial location $\ell$ resides
in one specific stratum, which implies that the stratum index $j$ is
completely determined by $\ell$. The biomass $y(\ell)$ for a location $\ell$
situated inside the $j$-th stratum is modeled as

$y(\ell)=\beta_{0}+\tilde{\beta}_{0}(\ell)+{\bm{x}}(\ell)^{\top}({\bm{\beta}}+\tilde{{\bm{\beta}}}(\ell))+w(\ell)+\epsilon(\ell),$
(2)

where $\beta_{0}$ is the intercept,
$\tilde{\beta}_{0}(\ell)=\tilde{\beta}_{0,j}$ for all $\ell$ in stratum $j$,
${\bm{x}}(\ell)$ is the $p\times 1$ vector of canopy structure variables with
associated global regression coefficients ${\bm{\beta}}$ and stratum-specific
coefficients $\tilde{{\bm{\beta}}}(\ell)=\tilde{{\bm{\beta}}}_{j}$ for all
$\ell$ in stratum $j$, $w(\ell)$ is a spatial random effect that adds local
adjustment with spatial dependence, and the residual is modeled as
$\epsilon(\ell)\overset{ind.}{\sim}N(0,\tau^{2}(\ell))$ with
$\tau^{2}(\ell)=\tau^{2}_{j}$ for all $\ell$ in stratum $j$. Stratum effects
are modeled as
$\tilde{\beta}_{0,j}\overset{\mathrm{iid}}{\sim}N(0,\sigma^{2}_{0})$ and
$\tilde{{\bm{\beta}}}_{j}=(\tilde{\beta}_{1,j},\tilde{\beta}_{2,j},\ldots,\tilde{\beta}_{p,j})^{\top}$
with $\tilde{\beta}_{k,j}\overset{\mathrm{iid}}{\sim}N(0,\sigma^{2}_{k})$ for
$k$ in $1,2,\ldots,p$. The spatial random effect is modeled as a Nearest
Neighbor Gaussian Process (NNGP) $w(\ell)\sim
NNGP(0,\sigma^{2}_{w}\rho(\cdot,\cdot;\phi_{w}))$, where $\sigma^{2}_{w}$ is
the variance, $\rho$ is a spatial correlation function defined for pairs of
locations within the domain, and $\phi_{w}$ captures correlation between
locations based on spatial separation (see, e.g., Datta et al.,, 2016;
Banerjee,, 2017). Briefly, the NNGP specification implies the spatial random
effect vector over $n$ locations
${\bm{w}}=(w(\ell_{1}),w(\ell_{2}),\ldots,w(\ell_{n}))^{\top}$ has probability
distribution $N({\bf 0},{\bm{K}}_{w})$, where ${\bm{K}}_{w}$ is an $n\times n$
NNGP covariance matrix derived from the covariance function
$K_{w}(\ell,\ell^{\prime};{\bm{\theta}}_{w})=\sigma_{w}^{2}\rho(\ell,\ell^{\prime};\phi_{w})=\sigma_{w}^{2}\exp(-\phi_{w}||\ell-\ell^{\prime}||)$
where $\exp(\cdot)$ is an exponential spatial correlation function,
${\bm{\theta}}_{w}=\\{\sigma_{w}^{2},\phi_{w}\\}$, and
$||\ell-\ell^{\prime}||$ is the Euclidean distance between possibly different
locations $\ell$ and $\ell^{\prime}$.

To allow prediction of $y(\ell)$ at locations where G-LiHT is not observed,
each canopy structure variable is modeled similar to biomass and includes both
stratum-varying regression coefficients and a continuous spatial process.
Specifically, the second stage model for the $k$-th canopy structure variable
at location $\ell$ situated inside the $j$-th stratum is

$x_{k}(\ell)=\alpha_{k,0}+\tilde{\alpha}_{k,0}(\ell)+{\bm{v}}_{k}(\ell)^{\top}({\bm{\alpha}}_{k}+\tilde{{\bm{\alpha}}}_{k}(\ell))+u_{k}(\ell)+\eta_{k}(\ell),$
(3)

where $\alpha_{k,0}$ is the intercept,
$\tilde{\alpha}_{k,0}(\ell)=\tilde{\alpha}_{k,0,j}$ for all $\ell$ in stratum
$j$, ${\bm{v}}_{k}(\ell)$ is a $r\times 1$ vector of predictors with
associated global regression coefficients ${\bm{\alpha}}_{k}$ and stratum-
specific effects $\tilde{{\bm{\alpha}}}_{k}(\ell)=\tilde{{\bm{\alpha}}}_{k,j}$
for all $\ell$ in stratum $j$, $u_{k}(\ell)$ is a spatial random effect, and
$\eta_{k}(\ell)\overset{ind.}{\sim}N(0,\gamma^{2}_{k}(\ell))$ with
$\gamma^{2}_{k}(\ell)=\gamma^{2}_{k,j}$ for all $\ell$ in stratum $j$. Stratum
effects are modeled as
$\tilde{\alpha}_{k,0,j}\overset{\mathrm{iid}}{\sim}N(0,\nu^{2}_{k,0}$) and
$\tilde{{\bm{\alpha}}}_{k,j}=(\tilde{\alpha}_{k,1,j},\tilde{\alpha}_{k,2,j},\ldots,\tilde{\alpha}_{k,r,j})^{\top}$
with $\tilde{\alpha}_{k,i,j}\overset{\mathrm{iid}}{\sim}N(0,\nu^{2}_{k,i})$
for $i$ in $1,2,\ldots,r$. The spatial random effect is modeled as
$u_{k}(\ell)\sim NNGP(0,\nu^{2}_{k,u}\rho(\cdot,\cdot;\phi_{k,u}))$, where
$\nu^{2}_{k,u}$ is the variance, $\rho$ is again taken as an exponential
spatial correlation function, and $\phi_{k,u}$ is the spatial decay parameter.
We again collect the spatial process parameters into a vector
${\bm{\theta}}_{k,u}=\\{\nu^{2}_{k,u},\phi_{k,u}\\}$.

For the set of $n$ locations $\mathcal{L}=(\ell_{1},\ell_{2},\ldots,\ell_{n})$
where both biomass and G-LiHT canopy structure data are observed, we define
the $n\times 1$ vector
${\bm{y}}=(y(\ell_{1}),y(\ell_{2}),\ldots,y(\ell_{n}))^{\top}$, $n\times 1$
vector of ones ${\bf 1}$, $n\times q$ matrix $\mathds{1}$ with the $(i,j)$-th
element equaling 1 if $\ell_{i}$ falls within the $j$-th stratum and zero
otherwise, $n\times p$ matrix ${\bm{X}}$ with row $i$ equal to
${\bm{x}}(\ell_{i})^{\top}$, $n\times pq$ matrix $\tilde{{\bm{X}}}$ with row
$i$ and columns $(j-1)p+1$ through $(j-1)p+p$ equal to
${\bm{x}}(\ell_{i})^{\top}$ when location $\ell_{i}$ falls within the $j$-th
stratum and zero otherwise,
${\bm{w}}=(w(\ell_{1}),w(\ell_{2}),\ldots,w(\ell_{n}))^{\top}$, and
${\bm{\epsilon}}=(\epsilon(\ell_{1}),\epsilon(\ell_{2}),\ldots,\epsilon(\ell_{n}))^{\top}$.
We then write model (2) as

${\bm{y}}=\beta_{0}{\bf
1}+\mathds{1}\tilde{{\bm{\beta}}}_{0}+{\bm{X}}{\bm{\beta}}+\tilde{{\bm{X}}}\tilde{{\bm{\beta}}}+{\bm{w}}+{\bm{\epsilon}},$
(4)

where $\beta_{0}$ and ${\bm{\beta}}$ are as defined earlier,
$\tilde{{\bm{\beta}}}_{0}=(\tilde{\beta}_{0,1},\tilde{\beta}_{0,2},\ldots,\tilde{\beta}_{0,q})^{\top}$,
and
$\tilde{{\bm{\beta}}}=(\tilde{{\bm{\beta}}}_{1}^{\top},\tilde{{\bm{\beta}}}_{2}^{\top},\ldots,\tilde{{\bm{\beta}}}_{q}^{\top})^{\top}$.
Let
${\bm{\Omega}}_{y}=\\{\beta_{0},\tilde{{\bm{\beta}}}_{0},{\bm{\beta}},\tilde{{\bm{\beta}}},\\{\sigma^{2}_{i}\\}_{i=0}^{p},{\bm{w}},{\bm{\theta}}_{w},\\{\tau_{j}^{2}\\}_{j=1}^{q}\\}$
denote the collection of parameters in the above model including those
specifying the prior distributions and the spatial process as described below
(2) to complete the hierarchical model.

Similarly, for the $n_{s}$ locations
$\mathcal{L}_{s}=(\ell_{1},\ell_{2},\ldots,\ell_{n_{s}})$ where G-LiHT is
observed, we define the $n_{s}\times 1$ vector for the $k$-th canopy structure
variable
${\bm{x}}_{k}=(x_{k}(\ell_{1}),x_{k}(\ell_{2}),\ldots,x_{k}(\ell_{n_{s}}))^{\top}$,
$n_{s}\times 1$ vector of ones ${\bf 1}_{k}$, $n_{s}\times q$ matrix
$\mathds{1}_{k}$ with row $i$ and column $j$ equal 1 when $\ell_{i}$ falls
within the $j$-th stratum and zero otherwise, $n_{s}\times r$ matrix
${\bm{V}}$ with row $i$ equal to ${\bm{v}}(\ell_{i})^{\top}$, $n_{s}\times rq$
matrix $\tilde{{\bm{V}}}$ with row $i$ and columns $(j-1)r+1$ through
$(j-1)r+r$ equal to ${\bm{v}}(\ell_{i})^{\top}$ when location $\ell_{i}$ falls
within the $j$-th stratum and zero otherwise,
${\bm{u}}_{k}=(u_{k}(\ell_{1}),u_{k}(\ell_{2}),\ldots,u_{k}(\ell_{n_{s}}))^{\top}$,
and
${\bm{\eta}}_{k}=(\eta_{k}(\ell_{1}),\eta_{k}(\ell_{2}),\ldots,\eta_{k}(\ell_{n_{s}}))^{\top}$.
We then write model (3) as

${\bm{x}}_{k}=\alpha_{k,0}{\bf
1}_{k}+\mathds{1}_{k}\tilde{{\bm{\alpha}}}_{k,0}+{\bm{V}}{\bm{\alpha}}_{k}+\tilde{{\bm{V}}}\tilde{{\bm{\alpha}}}_{k}+{\bm{u}}_{k}+{\bm{\eta}}_{k}.$
(5)

where
$\tilde{{\bm{\alpha}}}_{k,0}=(\tilde{\alpha}_{k,0,1},\tilde{\alpha}_{k,0,2},\ldots,\tilde{\alpha}_{k,0,q})^{\top}$,
and
$\tilde{{\bm{\alpha}}}_{k}=(\tilde{{\bm{\alpha}}}_{k,1}^{\top},\tilde{{\bm{\alpha}}}_{k,2}^{\top},\ldots,\tilde{{\bm{\alpha}}}_{k,q}^{\top})^{\top}$.
The spatial random effect vector ${\bm{u}}_{k}\sim NNGP({\bf
0},{\bm{K}}_{k,u})$, where ${\bm{K}}_{k,u}$ is an $n_{s}\times n_{s}$ NNGP
covariance matrix defined analogously to ${\bm{K}}_{w}$. Model parameters are
collected in the vector
${\bm{\Omega}}_{x_{k}}=(\alpha_{k,0},\tilde{{\bm{\alpha}}}_{k,0},{\bm{\alpha}}_{k},\tilde{{\bm{\alpha}}}_{k},\\{\nu^{2}_{k,i}\\}_{i=0}^{r},{\bm{u}}_{k},{\bm{\theta}}_{k,u},\\{\gamma^{2}_{k,j}\\}_{j=1}^{q})$,
which includes parameters specifying prior distributions in spatial process as
described below (3).

To complete the Bayesian specification of these models, a prior distribution
is assigned to each parameter. For model (2) and (3) we assign: flat priors to
$\beta_{0}$, $\alpha_{k,0}$ and elements in ${\bm{\beta}}$ and
${\bm{\alpha}}_{k}$; implicit Normal distributions for elements in
$\tilde{{\bm{\beta}}}_{0}$, $\tilde{{\bm{\beta}}}$,
$\tilde{{\bm{\alpha}}}_{k,0}$, and $\tilde{{\bm{\alpha}}}_{k}$; NNGP
distribution for ${\bm{w}}$ and ${\bm{u}}_{k}$; weakly informative inverse-
Gamma distributions with hyperparameter shape equal to 2 and scale value
guided by EDA results for all variance parameters; Uniform distributions with
support between 1 to 500 (km) for decay parameters $\phi_{w}$ and $\phi_{u}$.

Following the Bayesian mode of inference, we employ a Markov chain Monte Carlo
(MCMC) algorithm to generate samples from model (2) parameters’ joint
posterior distribution,

$[{\bm{\Omega}}_{y}\,|\,{\bm{y}},\mathds{1},{\bm{X}}]\propto[{\bm{\Omega}}_{y}]\times[{\bm{y}}\,|\,{\bm{\Omega}}_{y},\mathds{1},{\bm{X}}],$
(6)

where $[{\bm{\Omega}}_{y}\,|\,{\bm{y}},\mathds{1},{\bm{X}}]$ represents the
joint posterior distribution of ${\bm{\Omega}}_{y}$ conditioned on the data,
$[{\bm{\Omega}}_{y}]$ is the joint prior distribution, and
$[{\bm{y}}\,|\,{\bm{\Omega}}_{y},\mathds{1},{\bm{X}}]$ represents the
likelihood. Parameter inference is obtained by sampling from (6) using MCMC.
After diagnosing convergence, we collect $L$ samples from (6), which are
denoted as
$({\bm{\Omega}}_{y}^{(1)},{\bm{\Omega}}_{y}^{(2)},\ldots,{\bm{\Omega}}_{y}^{(L)})$.
Similarly, we obtain estimates for (3) by sampling from the joint posterior
distribution,

$[{\bm{\Omega}}_{x_{k}}\,|\,{\bm{x}}_{k},\mathds{1}_{k},{\bm{V}}]\propto[{\bm{\Omega}}_{x_{k}}]\times[{\bm{x}}_{k}\,|\,{\bm{\Omega}}_{x_{k}},\mathds{1}_{k},{\bm{V}}],$
(7)

with $L$ post-convergence MCMC samples collected as
$({\bm{\Omega}}_{x_{k}}^{(1)},{\bm{\Omega}}_{x_{k}}^{(2)},\ldots,{\bm{\Omega}}_{x_{k}}^{(L)})$.

###  3.1 Posterior predictive inference and areal estimates

As described in Section 1, our inferential objective is to estimate biomass
density (Mg/ha) and total (Mg) for any user-defined area within the TIU. This
inference uses model (2) to estimate the biomass density posterior predictive
distribution at each of $n^{\ast}$ prediction locations,
$\mathcal{L}_{0}=(\ell_{0,1},\ell_{0,2},\ldots,\ell_{0,n^{\ast}})$, laid in a
dense grid over the area of interest. These predictions are then used to
estimate the desired area summaries of biomass density and total. However,
because G-LiHT canopy structure variables are not observed at prediction
locations, we must condition biomass predictions from model (2) on predictions
of canopy structure variables from model (3). Ideally, this is done in a way
that propagates the uncertainty in canopy structure variable predictions
through to biomass predictions.

For inference at unobserved locations, we extend model (3) from observed
locations to any arbitrary location $\ell_{0}$. If
${\bm{x}}^{\ast}_{k}=(x^{\ast}_{k}(\ell_{0,1}),x^{\ast}_{k}(\ell_{0,2}),\ldots,x^{\ast}_{k}(\ell_{0,n^{\ast}}))^{\top}$
is the vector of unknown measurements and ${\bm{u}}_{k}^{\ast}$ is the
analogously defined vector of spatial random effects corresponding to the
$k$-th canopy structure variable, the NNGP extends model (3) such that if
$\ell_{0,m}$ falls within the $j$-th stratum then

$x^{\ast}_{k}(\ell_{0,m})=\alpha_{k,0}+\tilde{\alpha}_{k,0,j}+{\bm{v}}_{k}(\ell_{0,m})^{\top}({\bm{\alpha}}_{k}+\tilde{{\bm{\alpha}}}_{k,j})+u_{k}(\ell_{0,m})+\eta_{k}(\ell_{0,m})\,\;\quad
m=1,2,\ldots,n^{*}\;,$ (8)

with $\eta_{k}(\ell_{0,m})\overset{ind}{\sim}N(0,\gamma_{k,j}^{2})$ and an
NNGP predictive model for $u_{k}(\ell_{0,m})$. We briefly explain below and
refer the reader to further details in Banerjee, (2017) Section 3.2 and their
Equation (19).

It will be convenient, in what follows, to denote
${\bm{K}}(\mathcal{A},\mathcal{B})$ for finite sets of locations $\mathcal{A}$
and $\mathcal{B}$ to be the matrix whose $(i,j)$-th element is evaluated from
the covariance function $K(\ell_{i},\ell_{j})$, where $\ell_{i}$ and
$\ell_{j}$ are the $i$-th entry in $\mathcal{A}$ and $j$-th entry in
$\mathcal{B}$, respectively. We build a sequence of neighbor sets
$N(\ell_{0,m})$ for each $m=1,2,\ldots,n^{\ast}$ that consist of a fixed
number of nearest neighbors of $\ell_{0,m}$ from the “past,” where “past”
refers to the $n_{s}$ locations already accounted for in ${\cal L}_{s}$ and
the set of $\ell_{0,i}$’s for $i<m$. The distribution of $u_{k}(\ell_{0,m})$
is specified as

$u_{k}(\ell_{0,m})=\sum_{\ell^{\prime}\in
N(\ell_{0,m})}a(\ell_{0,m},\ell^{\prime})u_{k}(\ell^{\prime})+\omega_{k}(\ell_{0,m})\,\;\quad
m=1,2,\ldots,n^{*}\;,$ (9)

where solving
${\bm{K}}(N(\ell_{0,m}),N(\ell_{0,m})){\bm{a}}(\ell_{0,m})={\bm{K}}(N(\ell_{0,m}),\ell_{0,m})$
for the vector ${\bm{a}}(\ell_{0,m})$ renders the values of
$a(\ell_{0,m},\ell^{\prime})$ in (9), while $\omega_{k}(\ell_{0,m})\sim
N(0,\delta^{2}_{0,m})$ with

$\delta^{2}_{0,m}=K(\ell_{0,m},\ell_{0,m})-{\bm{K}}(\ell_{0,m},N(\ell_{0,m})){\bm{K}}(N(\ell_{0,m}),N(\ell_{0,m}))^{-1}{\bm{K}}(N(\ell_{0,m}),\ell_{0,m})\;.$

Predictive inference is achieved by sampling from
$[x_{k}(\ell_{0,m})\,|\,{\bm{x}}_{k}]$, where ${\bm{x}}_{k}$ is the collection
of observations for the $k$-th canopy structure variable. Since the joint
posterior distribution of the hierarchical model extended over ${\cal L}_{0}$
is
$[{\bm{\Omega}}_{x_{k}},{\bm{x}}_{k}^{\ast},{\bm{u}}_{k}^{\ast}\,|\,{\bm{x}}_{k}]\propto[{\bm{\Omega}}_{x_{k}}\,|\,{\bm{x}}_{k}]\times[{\bm{u}}_{k}^{\ast}\,|\,{\bm{\Omega}}_{x_{k}}]\times[{\bm{x}}_{k}^{*}\,|\,{\bm{u}}_{k}^{\ast},{\bm{\Omega}}_{x_{k}}]$,
we can use the posterior samples from
${\bm{\Omega}}_{x_{k}}^{(l)}\sim[{\bm{\Omega}}_{x_{k}}\,|\,{\bm{x}}_{k}]$ for
each $l=1,2,\ldots,L$ and each $\ell_{0,m}\in{\cal L}_{0}$ to draw
$u_{k}^{(l)}(\ell_{0,m})\sim[u_{k}(\ell_{0,m})\,|\,{\bm{\Omega}}_{x_{k}}^{(l)}]$
and
$x_{k}^{(l)}(\ell_{0,m})\sim[x_{k}(\ell_{0,m})\,|\,u_{k}^{(l)}(\ell_{0,m}),{\bm{\Omega}}_{x_{k}}^{(l)}]$.
The resulting $x_{k}^{(l)}(\ell_{0,m})$ are the desired posterior predictive
samples.

Predictive inference for $y(\ell)$ over ${\cal L}_{0}$ proceeds analogously.
We draw $y^{(l)}(\ell_{0,m})$ from
$[y(\ell_{0,m})\,|\,{\bm{\Omega}}_{y}^{(l)},{\bm{x}}^{(l)}(\ell_{0,m})]$ by
generating a value from (2) for each sampled ${\bm{\Omega}}_{y}^{(l)}$ and
${\bm{x}}^{(l)}(\ell_{0,m})$, where the latter is generated from its posterior
predictive distribution as above. We model biomass over an area $D$ as a
density, i.e., on a per unit area basis Mg/ha, which we denote by $y(D)$. The
posterior distribution for $y(D)$ is evaluated by calculating
$y^{(l)}(D)\approx\sum^{n^{\ast}}_{m=1}y^{(l)}(\ell_{0,m})/n^{\ast}$ for each
sampled $y^{(l)}(\ell_{0,m})$. Similarly, posterior samples of total biomass
are obtained from
$y_{Tot.}^{(l)}(D)=|D|\sum^{n^{\ast}}_{m=1}y^{(l)}(\ell_{0,m})/n^{\ast}$,
where $|D|$ is the area in ha. Samples from the posterior distribution of
$y(D)$ and $y_{Tot.}(D)$ are then summarized to described any quantity of
interest, e.g., posterior distribution mean, standard deviation, and 95%
credible interval.

###  3.2 Implementation

Methods developed in the preceding sections were programmed in C++ and used
openBLAS Zhang, (2016) and Linear Algebra Package (LAPACK;
www.netlib.org/lapack) for efficient matrix computations. openBLAS is an
implementation of Basic Linear Algebra Subprograms (BLAS; www.netlib.org/blas)
capable of exploiting multiple processors. Computing details about efficiently
updating NNGP spatial random effects are provided in Finley et al., (2019).
MCMC sampler code for (2) and (3) are provided in the supplemental material
along with simulated data used for testing the code and estimation procedures.

The computer used for subsequent analyses was running a linux operating system
with a AMD Ryzen Threadripper 3990X 64-Core Processor (128 threads) with 264
GB of RAM.

###  3.3 Submodels and model selection

Connecting Sections 2 and 3, for model (2) our analysis uses the NLCD strata
to form $\mathds{1}$ and single G-LiHT canopy structure variable $x_{CH}$ to
form ${\bm{X}}$ and $\tilde{{\bm{X}}}$ (i.e., $q=4$ and $p=1$). Model (3) for
$x_{CH}$ is similarly informed using the NLCD strata and a single variable,
$v_{TC}$, to form ${\bm{V}}$ and $\tilde{{\bm{V}}}$ (i.e., $r=1$).

We consider the full model (2) and four submodels. All submodels follow model
(2) notation and indexing unless noted otherwise.

Submodel 1:
$\displaystyle\;\;y(\ell)=\beta_{0}+x_{CH}(\ell)\beta_{CH}+\epsilon(\ell)\text{,
where }\epsilon(\ell)\overset{\mathrm{iid}}{\sim}N(0,\tau^{2}),$ Submodel 2:
$\displaystyle\;\;y(\ell)=\beta_{0}+x_{CH}(\ell)\beta_{CH}+\epsilon(\ell)\text{,
where }\epsilon(\ell)\overset{ind}{\sim}N(0,\tau^{2}_{j}),$ Submodel 3:
$\displaystyle\;\;y(\ell)=\beta_{0}+\tilde{\beta}_{0,j}+x_{CH}(\ell)(\beta_{CH}+\tilde{\beta}_{CH,j})+\epsilon(\ell)\text{,
where }\epsilon(\ell)\overset{\mathrm{iid}}{\sim}N(0,\tau^{2}),$ Submodel 4:
$\displaystyle\;\;y(\ell)=\beta_{0}+\tilde{\beta}_{0,j}+x_{CH}(\ell)(\beta_{CH}+\tilde{\beta}_{CH,j})+\epsilon(\ell)\text{,
where }\epsilon(\ell)\overset{ind}{\sim}N(0,\tau^{2}_{j}),$ Full model:
$\displaystyle\;\;y(\ell)=\beta_{0}+\tilde{\beta}_{0,j}+x_{CH}(\ell)(\beta_{CH}+\tilde{\beta}_{CH,j})+w(\ell)+\epsilon(\ell)\text{,
where }\epsilon(\ell)\overset{ind}{\sim}N(0,\tau^{2}_{j}).$

To assess the contribution of G-LiHT derived canopy structure information, we
also consider the full model without $x_{CH}$, i.e.,
$y(\ell)=\beta_{0}+\tilde{\beta}_{0,j}+w(\ell)+\epsilon(\ell)$.

We consider several criterion for selecting the “best” model. The deviance
information criterion (DIC; Spiegelhalter et al.,, 2002) and widely applicable
information criterion (WAIC; Watanabe,, 2010) model fit criterion were
computed for each candidate model. DIC equals $-2(\text{L}-\text{p}_{D})$
where L is goodness of fit and pD is a model penalty term viewed as the
effective number of parameters. Two WAIC criteria were computed based on the
log pointwise predictive density (LPPD) with
$\text{WAIC}_{1}=-2(\text{LLPD}-\text{p}_{1})$ and
$\text{WAIC}_{2}=-2(\text{LLPD}-\text{p}_{2})$ where penalty terms p1 and p2
are defined in Gelman et al., (2014) just prior to, and in, their Equation
(11). Models with lower DIC, WAIC1, and WAIC2 values have better fit to the
observed data and should yield better out-of-sample prediction, see Gelman et
al., (2014), Vehtari et al., (2017), or Green et al., (2020) for more details.
Lastly, we compute root mean squared error (RMSE) between the observed and
model fitted values.

##  4 Design-based and model-based estimates

Design- and model-based approaches follow different theories of inference.
Both are well developed and compared in statistical and forestry literature
(see, e.g., Särndal et al.,, 1978, 2003; Gregoire,, 1998; McRoberts,, 2010).
Little, (2004) offers an excellent review from diverse perspectives. The
design-based approach assumes a fixed finite population that can be accessible
(in principle without error) through a census if all population units were
observed. Randomness is incorporated via the selection of population units
into a sample according to a randomized sampling design. A sampling design
assigns a probability of selection to each sample. This is often effective
when the variability and dependence across the population units can be
adequately captured by the sampling design.

A comprehensive review and assessment of design-based and model-based
inference for survey data is beyond the scope of this paper. Nevertheless, an
increasingly prevalent view seems to be that if the units of the population
exhibit associations or dependencies that are too complex to be accounted for
by a sampling design, then a model-based approach to inference is preferable
(see, e.g., the developments in Little,, 2004; Ghosh,, 2012; Banerjee,, 2023,
for model-based as well as fully Bayesian perspectives to inference for finite
populations). Furthermore, it is well-known that popular design-based
estimators such as the Horvitz-Thompson (Horvitz and Thompson,, 1952) emerge
as special cases of certain hierarchical models (see, e.g., Section 2 in
Banerjee,, 2023). While design-based inference has been developed in SAE
problems (see, e.g., the text by Rao and Molina,, 2015), in applications such
as here the nature of the spatial structure in the data and complex
dependencies suggest that the population is a realization from a data-
generating stochastic process. As argued in Banerjee, (2023), introducing the
complex dependencies arising from a spatial process is difficult in design-
based weights. Post-stratification issues for spatial data are also
complicated. This renders the model-based approach more appropriate, where
randomness is incorporated through probability laws endowed by a spatial
process.

Design-based estimates presented in Section 6 are generated using the FIA
post-stratified estimators detailed in Bechtold and Patterson, (2005). Here,
the mean biomass within the $j$-th stratum is estimated as

$\bar{y}_{j}=\frac{\sum^{n_{j}}_{i=1}y_{j,i}}{n_{j}}$ (10)

with associated standard error of the mean estimated as

$s_{\bar{y}_{j}}=\sqrt{\frac{\sum^{n_{j}}_{i=1}\left(y_{i,j}-\bar{y}_{j}\right)^{2}}{n_{j}\left(n_{j}-1\right)}},$
(11)

where $n_{j}$ is the number of plots within the $j$-th stratum. The post-
stratified estimator combines stratum estimates using area-based weights. The
$j$-th stratum’s weight is $W_{j}=A_{j}/A$, where $A_{j}$ is the $j$-th
stratum area and $A$ is the area of all strata (i.e., $\sum^{q}_{j=1}A_{j}$).
The post-stratified mean is estimated as

$\bar{y}=\sum^{q}_{j=1}W_{j}\bar{y}_{j}$ (12)

with associated standard error of the mean estimated as

$s_{\bar{y}}=\sqrt{\frac{1}{n}\left(\sum^{q}_{j=1}W_{j}n_{j}s^{2}_{\bar{y}_{j}}+\sum^{q}_{j=1}\left(1-W_{j}\right)\frac{n_{h}}{n}s^{2}_{\bar{y}_{j}}\right)}.$
(13)

The different uses of randomization distributions between design- and model-
based inference yield different population parameter estimates, particularly
for uncertainty summaries of these estimates. In our setting, for example, a
design-based 95% confidence interval for the TIU’s biomass total is
interpreted as “if a large number of independent and equally sized samples
were collected according to the sampling design, 95% of the intervals computed
using these samples would include the population total.” In contrast, a
Bayesian model-based 95% credible interval is interpreted as “there is a 95%
chance or probability the interval includes the population total, given the
posited model and observed data.”

While the interpretation of these intervals is different, it is worth pointing
out two important benefits of the model-based approach in our application.
First, we are able to predict the biomass at arbitrary locations while
accounting for its dependence on canopy structure and using the spatial random
field associated with biomass. Design-based estimators, on the other hand,
will not offer inference beyond the finite population units. Second, our data
analysis reveals that exploiting the spatial structure and the use of a model
for canopy structure delivers improved uncertainty quantification in the
biomass for different strata. We defer to our specific discussion surrounding
the results in Tables 5 and 6.

##  5 Analyses

In Section 6, we present the full model estimates of biomass density and total
for the TIU. We also present the design-based post-stratified estimates
generated using the complete sample of 1 091 TIU plots. In addition to the
entire TIU, we consider two illustrative small areas—Caribou-Poker Creeks
Research Watershed (CPC) and Bonanza Creek Experimental Forest (BCEF). CPC and
BCEF are located on Alaska state land and run as long-term ecological research
sites by the University of Alaska, Fairbanks. The location and extent of these
small areas are shown in Figure 1a.

As illustrated in Figures 6a and 7a, CPC and BCEF have a continuous forest
inventory (CFI) plot network where each plot follows the FIA layout and
measurement protocol. These inventory data were used to generate design-based
post-stratified estimates that we compare with our model-based estimates.
Importantly, the TIU data inform the model-based SAEs, not the CPC and BCEF
CFI data. In this way, design-based CPC and BCEF estimates provide an
independent assessment using data separate from those used to inform the
model-based estimates.

The candidate models were fit using $n$=880 spatially coinciding FIA and
G-LiHT locations (Figure 1a). Model (3) parameter estimates and subsequent
$x_{CH}$ predictions were informed using the $n_{s}$=61 029 G-LiHT only
locations (Figure 1a). Estimates for TIU biomass density and total were based
on samples from posterior predictive distributions at $n^{\ast}$=2 165 220
locations (see, e.g., Figure 1b) laid out on a 250-by-250 (m) grid (one
prediction location per 6.25 ha). This chosen grid is sufficiently dense to
yield stable areal estimates and provide detailed maps (see Figure S1 for the
relationship between estimate stability and prediction grid density over the
TIU). Given CPC’s and BCEF’s small spatial extent and desire to provide more
detailed prediction maps, estimates were based on samples from posterior
predictive distributions at locations on a 50-by-50 (m) grid.

Three chains of each model’s MCMC sampler were run for 25 000 iterations. The
full model’s runtime for one chain was $\sim$2 hours, with all submodels
having substantially shorter runtimes. Posterior inference for all model
parameters and predictions was based on 1 000 post-burn and thinned MCMC
samples from each of three chains. Runtime to generate the 3 000 posterior
predictive samples over each of the TIU’s $n^{\ast}$ locations was $\sim$10
hours.

##  6 Results

###  6.1 Model selection

Candidate model fit criteria scores are given in Table 3. As suggested by
Section 2.4 EDA, Submodel 2’s improved fit over Submodel 1 and 3 supports
stratum specific residual variances (i.e., $\tau^{2}_{j}$ parameters).
Submodel 4, which includes stratum-varying coefficients and residual
variances, provides the best fit among the submodels. The full model, with its
spatially varying intercept, provides marginally better fit over that of
Submodel 4. Given the strong relationships between biomass and $x_{CH}$ seen
in the EDA, it is not surprising to see substantially degraded fit when this
canopy structure variable is not included in the full model (i.e., comparing
full model fit criteria with and without $x_{CH}$). Given the full model (2)
with $x_{CH}$ provides the best fit, it is used for all subsequent inference.

|  Submodels (with $x_{CH}$) | Full model
---|---|---
Criteria | 1 | 2 | 3 | 4 |  With $x_{CH}$ |  Without $x_{CH}$
WAIC1 | 7369.82 | 7100.65 | 7112.79 | 6921.25 | 6894.00 | 8105.62
WAIC2 | 7369.94 | 7102.80 | 7115.22 | 6924.03 | 6914.94 | 8047.47
p1 | 10.30 | 18.47 | 20.18 | 30.98 | 152.05 | 71.48
p2 | 10.36 | 19.54 | 21.39 | 32.37 | 162.52 | 42.40
LPPD | -3674.61 | -3531.86 | -3536.22 | -3429.64 | -3294.95 | -3981.33
DIC | 7362.51 | 7088.07 | 7099.19 | 6899.82 | 6864.26 | 8056.36
pD | 2.99 | 5.89 | 6.57 | 9.56 | 122.31 | 22.21
RMSE | 16.84 | 15.97 | 14.41 | 13.68 | 12.98 | 26.90
Table 3: Fit criterion for candidate models. Scores that indicate the “best”
model are bolded.

###  6.2 Tanana Inventory Unit estimates

Estimates for the full model’s stratum-varying regression coefficients and
residual variance parameters are given in Table 4. The remainder of this
model’s parameter estimates are given in Table S1. As expected, coefficients’
sign and magnitude generally follow EDA stratum specific regression
coefficient estimates in Table 2 and show $x_{CH}$ explains a substantial
portion of variability in biomass (i.e., 95% credible interval bound for each
slope coefficient does not include zero). Compared with EDA models, the full
model’s information pooling, via stratum and spatial random effects, yields
smaller and more precise residual variance estimates (i.e., $\tau^{2}_{j}$
parameters). Spatial process parameter estimates in Table S1 show spatial
random effects capture a relatively small amount of residual variation (i.e.,
$\sigma^{2}_{w}$ to $\tau^{2}_{j}$s ratios are small). Spatial decay parameter
estimates suggest mean effective spatial range is 37.92 (km) with lower and
upper 95% credible intervals 25.39 and 115.22 (km), respectively (we define
the “effective spatial range” as the distance at which the spatial correlation
drops below 0.05, which equals $-\log(0.05)/\phi_{w}$).

| $\beta_{0}+\tilde{\beta}_{0,j}$ | $\beta_{CH}+\tilde{\beta}_{CH,j}$ | $\tau^{2}_{j}$
---|---|---|---
Stratum | ${}_{\text{(L. 95\%)}}\;\;\text{Mean}\;\;_{\text{(U. 95\%)}}$ | ${}_{\text{(L. 95\%)}}\;\;\text{Mean}\;\;_{\text{(U. 95\%)}}$ | ${}_{\text{(L. 95\%)}}\;\;\text{Mean}\;\;_{\text{(U. 95\%)}}$
Conifer | ${}_{(\text{-}4.213)}\;\;\text{-}1.735\;\;_{(0.879)}$ | ${}_{(12.219)}\;\;12.945\;\;_{(13.699)}$ | ${}_{(257.461)}\;\;314.385\;\;_{(382.148)}$
Deciduous | ${}_{(\text{-}5.392)}\;\;\text{-}2.116\;\;_{(1.118)}$ | ${}_{(8.288)}\;\;8.896\;\;_{(9.493)}$ | ${}_{(246.176)}\;\;352.329\;\;_{(547.472)}$
Mixed | ${}_{(\text{-}4.824)}\;\;\text{-}1.945\;\;_{(0.784)}$ | ${}_{(8.221)}\;\;9.009\;\;_{(9.770)}$ | ${}_{(97.547)}\;\;167.053\;\;_{(294.068)}$
Other | ${}_{(\text{-}3.194)}\;\;\text{-}1.943\;\;_{(\text{-}0.849)}$ | ${}_{(7.291)}\;\;7.721\;\;_{(8.145)}$ | ${}_{(46.736)}\;\;69.811\;\;_{(84.513)}$
Table 4: Biomass model (2) parameter estimates, with $j$ indexing stratum.
Associated process parameter estimates are given in Table S1.

Parameter estimates for $x_{CH}$’s model (3) are given in Tables S2 and S3.
Here, estimates of stratum-varying regression coefficients suggest that both
strata and $v_{TC}$ explain a substantial portion of variability in $x_{CH}$.
The model’s spatial random effect explains a portion of residual variation
(i.e., $\nu^{2}_{CH,u}$ to $\gamma^{2}_{j}$s ratios are larger); however, the
effective spatial range for the spatial process is short, extending only about
5 (km). This is not too surprising given $v_{TC}$ explains a substantial
portion of $x_{CH}$’s variability.

(a)

(b)

Figure 4: LABEL:sub@fig:tiu-xmu and LABEL:sub@fig:tiu-xsd mean and standard
deviation, respectively, of the posterior predictive distribution for the
G-LiHT mean canopy height variable estimated using (3) and associated
predictive model.

The predictive model for $x_{CH}$ (8) yields a posterior predictive
distribution at the TIU’s $n^{*}$ prediction locations. The posterior
predictive distribution mean and standard deviation for $x_{CH}$ at each
prediction location is given in Figures 4a and 4b, respectively. Given the
short effective spatial range, approximately 5 (km), the G-LiHT flight lines
where $x_{CH}$ was measured appear as increased fidelity and precision stripes
in the maps—predictions are improved along and adjacent to flight lines.

(a)

(b)

Figure 5: LABEL:sub@fig:tiu-ymu and LABEL:sub@fig:tiu-ysd mean and standard
deviation, respectively, of the posterior predictive distribution for biomass
estimated using (2) and associated predictive model.

The two stage modeling approach developed in Section 3 is designed to use
$x_{CH}$’s posterior predictive information to inform biomass predictions.
This information, propagated across model components, is seen as increased
fidelity and precision along flight lines in biomass posterior predictive
distribution summary maps Figures 5a and 5b, respectively.

Turning now to area biomass estimates. Table 5 gives the stratum area, sample
size, and design-based biomass density and total estimates. Corresponding
model-based estimates derived from the $n^{*}$ posterior predictive
distributions are given in Table 6. The estimators yield similar biomass
density point estimates. The Mixed stratum shows the largest disparity, with a
design-based estimate of 52.834 (Mg/ha) and model-based estimate of 40.288
(Mg/ha).

The model-based point estimate for TIU biomass total is 9 831.88 (1 000 Mg)
more than the design-based point estimate (i.e.,
$261\,441.680-251\,609.800=9\,831.88$). The Other stratum has the largest land
area in the TIU. Hence, the seemingly small difference in biomass density
between the design- and model-based estimates (i.e., 4.624 and 7.525,
respectively) results in a large difference in biomass total (i.e., 39 748.36
and 64 689.31, respectively). Difference between the design- and model-based
estimate for the Other stratum comprises about 62% of the difference between
the estimators’ TIU biomass total estimate, i.e., 62% of the 9 831.88 (1 000
Mg) is due to different Other density estimates (contribution of the remaining
strata are Conifer 18%, Deciduous 7%, and Mixed 12%).

Comparing stratum specific and TIU-wide design-based standard error to
corresponding model-based posterior predictive distribution standard deviation
in Tables 5 and 6 shows the model-based approach consistently yields more
precise uncertainty quantification.

|  |  | Biomass (Mg/ha) | Biomass (1000 Mg)
---|---|---|---|---
Stratum | Area (1000 ha) | $n$ | Mean | SE | Total | SE
Conifer | 3659.234 | 278 | 36.183 | 2.240 | 132402.390 | 8197.627
Deciduous | 890.847 | 58 | 66.259 | 6.171 | 59026.840 | 5497.304
Mixed | 386.728 | 31 | 52.834 | 7.025 | 20432.210 | 2716.658
Other | 8596.150 | 724 | 4.624 | 0.542 | 39748.360 | 4655.796
TIU | 13532.960 | 1091 | 18.592 | 0.804 | 251609.800 | 10882.700
Table 5: Tanana Inventory Unit strata area, number of FIA plots $n$, and design-based estimates for biomass density and total with associated standard error (SE). | Biomass (Mg/ha) | Biomass (1000 Mg)
---|---|---
Stratum | ${}_{\text{(Lower 95\%)}}\;\;\text{Mean}\;\;_{\text{(Upper 95\%)}}$ | SD | ${}_{\text{(Lower 95\%)}}\;\;\text{Mean}\;\;_{\text{(Upper 95\%)}}$ | SD
Conifer | ${}_{(32.312)}\;\;34.188\;\;_{(36.502)}$ | 1.040 | ${}_{(118238.780)}\;\;125102.840\;\;_{(133568.720)}$ | 3806.029
Deciduous | ${}_{(58.907)}\;\;62.939\;\;_{(67.118)}$ | 2.179 | ${}_{(52476.700)}\;\;56069.070\;\;_{(59791.490)}$ | 1941.281
Mixed | ${}_{(35.632)}\;\;40.288\;\;_{(44.415)}$ | 2.163 | ${}_{(13779.930)}\;\;15580.460\;\;_{(17176.670)}$ | 836.452
Other | ${}_{(6.629)}\;\;7.525\;\;_{(8.571)}$ | 0.529 | ${}_{(56981.520)}\;\;64689.310\;\;_{(73676.780)}$ | 4544.735
TIU | ${}_{(18.312)}\;\;19.319\;\;_{(20.239)}$ | 0.513 | ${}_{(247810.560)}\;\;261441.680\;\;_{(273898.360)}$ | 6945.950
Table 6: Tanana Inventory Unit model (2) biomass density and total posterior
predictive distribution mean and standard deviation (SD).

###  6.3 Small area estimates

Next we consider the CPC and BCEF analysis results. As noted at the end of
Section 3.3, posterior inference comes directly from the TIU data and model;
however, each small area uses a denser prediction location grid for posterior
predictive inference and hence area estimates (i.e., one prediction every 0.25
ha for the small areas vs. 6.25 ha for the TIU). The denser grid yields more
detailed maps and, given the fairly large size of each small area, has
negligible effect on subsequent biomass density and total estimates (see,
e.g., Figure S1).

We begin with the CPC, which is the northernmost small area shown in Figure
1a. Figure 6a shows the CPC boundary, locations with G-LiHT derived $x_{CH}$,
single FIA plot, and CFI plot network. The strata are shown in Figure 6b. The
posterior predictive distribution mean and standard deviation for $x_{CH}$ at
each prediction location is given in Figures 6c and 6d, respectively. The
increased precision seen in the broader TIU $x_{CH}$ standard deviation map
Figure 4b is less pronounced in Figure 6d but still visible where the G-LiHT
only observations coincide with the Other stratum (the large ratio of _signal_
variance $\nu^{2}_{CH,u}$ to _noise_ variance $\gamma^{2}_{j}$, where $j$
equals the Other stratum index, allows the random effect influence to be more
apparent). Variability in $x_{CH}$ within a given stratum is due primarily to
variability in tree cover $v_{TC}$ values. Biomass posterior predictive
distribution mean and standard deviation maps are given in Figures 6e and 6f,
respectively. These maps make apparent the stratum-varying biomass density and
variability, and propagated $x_{CH}$ uncertainty.

a) (a)  |  b) (b)
---|---
c) (c)  |  d) (d)
e) (e)  |  f) (f)
Figure 6: Caribou-Poker Creek Research Watershed (CPC) data and analysis
results. LABEL:sub@fig:cpcrw-data locations of G-LiHT, FIA, and CPC’s
continuous forest inventory (CFI) data. LABEL:sub@fig:cpcrw-strata strata use
for design- and model-based biomass estimates. LABEL:sub@fig:cpcrw-mux and
LABEL:sub@fig:cpcrw-sdx mean and standard deviation, respectively, of the
posterior predictive distribution for the G-LiHT mean canopy height variable
estimated using (3) and associated predictive model (8). LABEL:sub@fig:cpcrw-
muy and LABEL:sub@fig:cpcrw-sdy mean and standard deviation, respectively, of
the posterior predictive distribution for biomass estimated using (2) and
associated predictive model.

Table 7 gives the CPC’s stratum areas, CFI sample size, and design-based post-
stratified biomass density and total estimates. Comparing the CPC and TIU
design-based estimates (i.e., Tables 5 and 7), we see the CPC’s stratum
specific biomass density differs substantially from that of the broader TIU.
For example, CPC’s Conifer and Mixed strata have about half the density of the
TIU, i.e., 15.261 vs. 36.183 (Mg/ha) and 28.434 vs. 52.834 (Mg/ha),
respectively. Whereas the CPC’s Deciduous and Other strata have a greater
density than the TIU, i.e., 73.195 vs. 66.259 (Mg/ha) and 12.263 vs. 4.625
(Mg/ha). CPC model-based estimates derived from the $n^{*}$ posterior
predictive distributions are given in Table 8. Relative to the CPC’s design-
based biomass density point estimates, the model-based estimates look more
similar to the broader TIU. This is not surprising, because the model draws
information from the entire TIU dataset to inform CPC biomass estimates.
Despite the differences between the design- and model-based stratum specific
density and total estimates seen when comparing Tables 7 and 8, the CPC-wide
densities and totals are quite similar, i.e., 32.284 vs. 37.695 (Mg/ha) and
342.341 vs. 399.725 (1 000 Mg), respectively.

Comparing the CPC design-based standard errors in Table 7 with model-based
standard deviations in Table 8 shows the model delivers slightly larger but
comparable stratum specific and CPC-wide uncertainty quantification despite
having only one FIA plot located within the CPC.

|  |  | Biomass (Mg/ha) | Biomass (1000 Mg)
---|---|---|---|---
Stratum | Area (ha) | $n$ | Mean | SE | Total | SE
Conifer | 3793.290 | 13 | 15.261 | 5.378 | 57.888 | 20.400
Deciduous | 2922.567 | 10 | 73.195 | 6.398 | 213.918 | 18.698
Mixed | 1413.225 | 6 | 28.434 | 7.290 | 40.183 | 10.303
Other | 2475.102 | 6 | 12.263 | 5.669 | 30.352 | 14.032
CPC | 10604.180 | 35 | 32.284 | 3.217 | 342.341 | 34.118
Table 7: Caribou-Poker Creek Research Watershed strata area, number of inventory plots $n$, and design-based estimates for biomass density and total with associated standard error (SE). | Biomass (Mg/ha) | Biomass (1000 Mg)
---|---|---
Stratum | ${}_{\text{(Lower 95\%)}}\;\;\text{Mean}\;\;_{\text{(Upper 95\%)}}$ | SD | ${}_{\text{(Lower 95\%)}}\;\;\text{Mean}\;\;_{\text{(Upper 95\%)}}$ | SD
Conifer | ${}_{(22.810)}\;\;36.338\;\;_{(49.723)}$ | 7.125 | ${}_{(86.524)}\;\;137.839\;\;_{(188.614)}$ | 27.026
Deciduous | ${}_{(52.516)}\;\;64.087\;\;_{(75.373)}$ | 5.983 | ${}_{(153.481)}\;\;187.300\;\;_{(220.281)}$ | 17.486
Mixed | ${}_{(24.555)}\;\;36.604\;\;_{(48.815)}$ | 6.240 | ${}_{(34.702)}\;\;51.730\;\;_{(68.986)}$ | 8.818
Other | ${}_{(1.907)}\;\;9.247\;\;_{(16.340)}$ | 3.645 | ${}_{(4.721)}\;\;22.886\;\;_{(40.444)}$ | 9.021
CPC | ${}_{(27.620)}\;\;37.695\;\;_{(47.904)}$ | 5.350 | ${}_{(292.884)}\;\;399.725\;\;_{(507.978)}$ | 56.731
Table 8: Caribou-Poker Creek Research Watershed model (2) biomass density and
total posterior predictive distribution mean and standard deviation (SD).

Next we consider BCEF analysis results. The BCEF is the southernmost small
area shown in Figure 1a. Figure 7a shows the BCEF boundary, locations with
G-LiHT derived $x_{CH}$, two FIA plots, and CFI plot network. The strata are
shown in Figure 7b. The posterior predictive distribution mean and standard
deviation for $x_{CH}$ at each prediction location is given in Figures 7c and
7d, respectively. Like the CPC, increased precision in $x_{CH}$ prediction
(seen in Figure 7d) is most apparent in the Other stratum (due to the larger
signal to noise ratio, i.e., $\nu^{2}_{CH,u}/\gamma^{2}_{j}$, where $j$ equals
the Other stratum index). Biomass posterior predictive distribution mean and
standard deviation maps are given in Figures 7e and 7f, respectively.

a) (a)  |  b) (b)
---|---
c) (c)  |  d) (d)
e) (e)  |  f) (f)
Figure 7: Bonanza Creek Experimental Forest (BCEF) data and analysis results.
LABEL:sub@fig:bnz-data locations of G-LiHT, FIA, and BCEF’s continuous forest
inventory (CFI) data. LABEL:sub@fig:bnz-strata strata use for design- and
model-based biomass estimates. LABEL:sub@fig:bnz-mux and LABEL:sub@fig:bnz-sdx
mean and standard deviation, respectively, of the posterior predictive
distribution for the G-LiHT mean canopy height variable estimated using (3)
and associated predictive model (8). LABEL:sub@fig:bnz-muy and
LABEL:sub@fig:bnz-sdy mean and standard deviation, respectively, of the
posterior predictive distribution for biomass estimated using (2) and
associated predictive model.

Table 9 gives the BCEF’s stratum areas, CFI sample size, and design-based
post-stratified biomass density and total estimates. The BCEF and TIU design-
based biomass density estimates (i.e., Tables 5 and 9) are similar, with the
exception of the Other stratum. The Other point estimate for the BCEF is
18.070 vs. TIU’s 4.624 (Mg/ha). BCEF’s model-based estimates derived from the
$n^{*}$ posterior predictive distributions are given in Table 10. The BCEF’s
model-based biomass density point estimates for Conifer and Deciduous stratum
are a bit higher than the design-based estimates. Despite differences between
design- and model-based stratum specific density and total estimates seen when
comparing Tables 9 and 10, BCEF-wide densities and totals are quite similar,
i.e., 41.397 vs. 49.391 (Mg/ha) and 867.503 vs. 1 198.903 (1 000 Mg),
respectively.

Comparing the BCEF design-based standard errors in Table 9 with model-based
standard deviations in Table 10 shows the model delivers consistently more
precises uncertainty quantification across stratum specific and BCEF-wide
estimates.

|  |  | Biomass (Mg/ha) | Biomass (1000 Mg)
---|---|---|---|---
Stratum | Area (1000 ha) | $n$ | Mean | SE | Total | SE
Conifer | 7701.783 | 21 | 39.139 | 6.872 | 301.439 | 52.926
Deciduous | 5615.655 | 30 | 60.951 | 7.394 | 342.278 | 41.520
Mixed | 2124.105 | 8 | 58.444 | 17.532 | 124.140 | 37.239
Other | 5514.323 | 17 | 18.070 | 6.677 | 99.645 | 36.820
BCEF | 20955.870 | 76 | 41.397 | 4.157 | 867.503 | 87.106
Table 9: Bonanza Creek Experimental Forest strata area, number of inventory plots $n$, and design-based estimates for biomass density and total with associated standard error (SE). | Biomass (Mg/ha) | Biomass (1000 Mg)
---|---|---
Stratum | ${}_{\text{(Lower 95\%)}}\;\;\text{Mean}\;\;_{\text{(Upper 95\%)}}$ | SD | ${}_{\text{(Lower 95\%)}}\;\;\text{Mean}\;\;_{\text{(Upper 95\%)}}$ | SD
Conifer | ${}_{(43.211)}\;\;54.079\;\;_{(64.846)}$ | 5.396 | ${}_{(332.800)}\;\;416.503\;\;_{(499.429)}$ | 41.557
Deciduous | ${}_{(63.168)}\;\;71.335\;\;_{(80.474)}$ | 4.440 | ${}_{(354.728)}\;\;400.591\;\;_{(451.914)}$ | 24.934
Mixed | ${}_{(44.372)}\;\;54.328\;\;_{(63.964)}$ | 5.275 | ${}_{(94.252)}\;\;115.397\;\;_{(135.866)}$ | 11.204
Other | ${}_{(12.230)}\;\;18.600\;\;_{(25.275)}$ | 3.326 | ${}_{(67.440)}\;\;102.569\;\;_{(139.372)}$ | 18.338
BCEF | ${}_{(41.741)}\;\;49.391\;\;_{(57.211)}$ | 3.936 | ${}_{(874.711)}\;\;1035.022\;\;_{(1198.903)}$ | 82.483
Table 10: Bonanza Creek Experimental Forest model (2) biomass density and
total posterior predictive distribution mean and standard deviation (SD).

##  7 Discussion and future work

We obtain stratum specific biomass estimates for arbitrarily sized areas,
e.g., large area TIU and small area CPC and BCEF. From a SAE standpoint, the
unit-level model we propose also provides maps for the given areas and at a
user-defined spatial resolution. Under the Bayesian paradigm used here, these
maps (and resulting areal estimates) can summarize a posterior predictive
distribution characteristic, e.g., mean, median, standard deviation, or
credible interval. Further, in the same way we condition biomass prediction on
$x_{CH}$, access to biomass posterior predictive samples facilitates
uncertainty propagation through other functions or models that take biomass as
input (e.g., economic or ecological models).

Model assessment metrics used here suggest substantive information gain using
G-LiHT derived mean canopy height, stratum-varying effects, a space-varying
intercept, and stratum specific residual variance parameters. These features
were identified using EDA. One might also consider allowing for stratum-
varying spatial random effects (see, e.g., Chan-Golston et al.,, 2020);
however, given the paucity of spatial structure after accounting for the mean
canopy height, it is unlikely stratum-varying spatial processes are warranted.

Section 3.3 reveals that our model-based approach offers more precise
uncertainty quantification for biomass than design-based inference (compare
standard errors in Tables 5 with posterior standard deviations in Table 6) by
exploiting the spatial structure in the data and incorporating information
from canopy structure. With the acknowledgment that the true population
parameter is unknown (and perhaps unknowable), we compare population point
estimates. The design-based post-stratified point estimate, using 1 091 FIA
plots, puts TIU’s total biomass at 251 609.8 (1 000 Mg). The model-based point
estimate, using 880 FIA plots and auxiliary information, puts TIU’s total
biomass at a slightly higher 261 441.68 (1 000 Mg). Much of the difference in
these estimates is due to the model placing more biomass in the Other stratum.
Looking at the observed distribution of biomass in Other (Figure 2a) and
canopy cover map (Figure 1d), future work should consider a refined set of
strata that partition Other into non-vegetation (e.g., water, barren) and
vegetation. Or, following Finley et al., (2011) and May et al., (2024), one
might introduce an additional hierarchical level to differentiate between
locations with and without biomass.

We demonstrate the model-based approach is particularly useful for SAE. Here,
paucity of FIA plots in CPC and BCEF precludes design-based estimation. In
comparison, the proposed model draws on the broader TIU dataset, to deliver
CPC and BCEF high-resolution maps and areal estimates by stratum. A key
contribution here is demonstrating how a two-stage hierarchical Bayesian model
allows information to be shared between multiple model components and how
sparse data is used to inform parameter estimates within and across
components, and ultimately to prediction. This model uses information where
available to improve prediction accuracy and precision. For example, biomass
prediction near G-LiHT measurements is seen in Figures 5a and 4b as stripes of
higher fidelity and increased precision, respectively. Away from G-LiHT flight
lines, biomass prediction retreats to the stratum mean and its variance is
dominated by that of the non-spatial residual process.

Our model is trained by the observed plots where both biomass and G-LiHT
measurements are reported. We anticipate future data collection efforts in
interior Alaska (beyond the TIU) will also align G-LiHT collection with all
inventory plots. Since model-based inference is not limited to probability
sampling, future work could explore adaptive sampling designs given objective
functions that balance data acquisition cost at each model level (e.g., LiDAR
and plot measurements) with inferential objectives (e.g., maximizing
prediction accuracy and precision) (see, e.g., Xia et al.,, 2006; Mateu and
Müller,, 2012). Given a posited model, like the one presented here, such
adaptive sampling designs could facilitate cost-effective data collection
efforts in remote regions while meeting inferential objectives.

Finally, we anticipate substantial opportunities for expanding our modeling
framework in future research. A natural extension would be to build
multivariate process models (see, e.g., Ch.9 in Banerjee et al.,, 2014) that
will jointly model biomass and an ensemble of data on canopy structures. These
models, while offering richness in modeling through the use of cross-
covariance functions or linear combinations of latent processes (Gelfand and
Banerjee,, 2010), will need to reckon with a possibly large dimension of the
response vector comprising all measurements on biomass and all of the canopy
structures being considered. Here, the expanding literature on multivariate
spatial factor models on very large domains, including for spatially
misaligned data, offer exciting possibilities for further investigation (Ren
and Banerjee,, 2013; Taylor-Rodriguez et al.,, 2019; Zhang and Banerjee,,
2022).

## Acknowledgments

Funding was provided by: NASA Carbon Monitoring System (CMS) grants Hayes (CMS
2020) and Cook (CMS 2018); National Science Foundation (NSF) DMS-1916395;
joint venture agreement with the USDA Forest Service Forest Inventory and
Analysis 22-CA-11221638-201; USDA Forest Service, Region 9, Forest Health
Protection, Northern Research Station; Michigan State University
AgBioResearch.

## Supplemental Material

Parameter | ${}_{\text{(L. 95\%)}}\;\;\text{Mean}\;\;_{\text{(U. 95\%)}}$
---|---
$\sigma^{2}_{0}$ | ${}_{(0.332)}\;\;1.029\;\;_{(4.767)}$
$\sigma^{2}_{CH}$ | ${}_{(1.181)}\;\;3.124\;\;_{(10.130)}$
$\phi_{w}$ | ${}_{(0.026)}\;\;0.079\;\;_{(0.118)}$
$\sigma^{2}_{w}$ | ${}_{(6.767)}\;\;16.326\;\;_{(41.274)}$
Table S1: Process parameter estimates for the full model with $x_{CH}$. Biomass model (2) process parameter estimates. Associated parameter estimates are given in Table 4. | $\alpha_{CH,0}+\tilde{\alpha}_{CH,0,j}$ | $\alpha_{CH,TC}+\tilde{\alpha}_{CH,TC,j}$ | $\gamma^{2}_{j}$
---|---|---|---
Stratum | ${}_{\text{(L. 95\%)}}\;\;\text{Mean}\;\;_{\text{(U. 95\%)}}$ | ${}_{\text{(L. 95\%)}}\;\;\text{Mean}\;\;_{\text{(U. 95\%)}}$ | ${}_{\text{(L. 95\%)}}\;\;\text{Mean}\;\;_{\text{(U. 95\%)}}$
Conifer | ${}_{(2.431)}\;\;2.502\;\;_{(2.587)}$ | ${}_{(2.416)}\;\;2.484\;\;_{(2.559)}$ | ${}_{(2.078)}\;\;2.147\;\;_{(2.221)}$
Deciduous | ${}_{(5.898)}\;\;6.035\;\;_{(6.179)}$ | ${}_{(5.903)}\;\;6.025\;\;_{(6.146)}$ | ${}_{(11.500)}\;\;12.081\;\;_{(12.710)}$
Mixed | ${}_{(3.689)}\;\;3.852\;\;_{(4.013)}$ | ${}_{(3.687)}\;\;3.825\;\;_{(3.978)}$ | ${}_{(6.792)}\;\;7.307\;\;_{(7.864)}$
Other | ${}_{(1.555)}\;\;1.617\;\;_{(1.695)}$ | ${}_{(1.533)}\;\;1.597\;\;_{(1.668)}$ | ${}_{(0.315)}\;\;0.340\;\;_{(0.362)}$
Table S2: Mean canopy height $x_{CH}$ model (3) parameter estimates, with $j$ indexing stratum. Associated process parameter estimates are given in Table S3. Parameter | ${}_{\text{(L. 95\%)}}\;\;\text{Mean}\;\;_{\text{(U. 95\%)}}$
---|---
$\nu^{2}_{CH,0}$ | ${}_{(0.797)}\;\;1.927\;\;_{(6.345)}$
$\nu^{2}_{CH,TC}$ | ${}_{(0.826)}\;\;1.907\;\;_{(5.836)}$
$\phi_{CH,u}$ | ${}_{(0.560)}\;\;0.585\;\;_{(0.599)}$
$\nu^{2}_{CH,u}$ | ${}_{(4.666)}\;\;4.823\;\;_{(4.994)}$
Table S3: Mean canopy height $x_{CH}$ model (3) process parameter estimates.
Associated parameter estimates are given in Table S2. Figure S1: Stratum
specific and TIU-wide biomass posterior predictive distribution mean and
standard deviation (SD) for a range of prediction grid densities. The x-axis
ranges from one prediction location per 6.25 (ha) to one prediction location
per 50 000 (ha). Estimates in Table 6 are based on one prediction location per
6.25 (ha). These figures show stratum specific mean and SD estimates are
stable for densities of one prediction location per $\sim$50 (ha) and greater
and hence support the choice of prediction grid density of one location per
6.25 (ha) for results reported in Section 6.

## References

  * (1) Andersen, H.-E., Barrett, T., Winterberger, K., Strunk, J., and Temesgen, H. (2009a).  Estimating forest biomass on the western lowlands of the kenai peninsula of alaska using airborne lidar and field plot data in a model-assisted sampling design.  IUFRO Div. 4 Symposium, Extending Forest Inventory and Monitoring.
  * (2) Andersen, H.-E., Clarkin, T., Winterberger, K., and Strunk, J. (2009b).  An Accuracy Assessment of Positions Obtained Using Survey- and Recreational-Grade Global Positioning System Receivers across a Range of Forest Conditions within the Tanana Valley of Interior Alaska.  Western Journal of Applied Forestry, 24(3):128–136.
  * Andersen et al., (2024) Andersen, H.-E., Ståhlhl, G., Finley, A. O., and Cook, B. D. (2024).  Poststratified ratio estimators for model-assisted estimation in sample-based airborne surveys.  Should be on Arxiv by proof stage.
  * Andersen et al., (2022) Andersen, H.-E., Strunk, J., and McGaughey, R. J. (2022).  Using high-performance global navigation satellite system technology to improve forest inventory and analysis plot coordinates in the pacific region.  Technical report.
  * Babcock et al., (2018) Babcock, C., Finley, A. O., Andersen, H.-E., Pattison, R., Cook, B. D., Morton, D. C., Alonzo, M., Nelson, R., Gregoire, T., Ene, L., Gobakken, T., and Næsset, E. (2018).  Geostatistical estimation of forest biomass in interior alaska combining landsat-derived tree cover, sampled airborne lidar and field observations.  Remote Sensing of Environment, 212:212–230.
  * Banerjee, (2017) Banerjee, S. (2017).  High-Dimensional Bayesian Geostatistics.  Bayesian Analysis, 12(2):583 – 614.
  * Banerjee, (2023) Banerjee, S. (2023).  Finite population survey sampling: An unapologetic Bayesian perspective.
  * Banerjee et al., (2014) Banerjee, S., Carlin, B. P., and Gelfand, A. E. (2014).  Hierarchical modeling and analysis for spatial data.  CRC Press, Boca Raton, FL.
  * Bechtold and Patterson, (2005) Bechtold, W. A. and Patterson, P. L. (2005).  The enhanced forest inventory and analysis program-national sampling design and estimation procedures.  Gen. Tech. Rep. SRS-80. Asheville, NC: US Department of Agriculture, Forest Service, Southern Research Station. 85 p., 80.
  * Breidenbach and Astrup, (2012) Breidenbach, J. and Astrup, R. (2012).  Small area estimation of forest attributes in the Norwegian National Forest Inventory.  European Journal of Forest Research, 131(4):1255–1267.
  * Breidt and Opsomer, (2017) Breidt, F. J. and Opsomer, J. D. (2017).  Model-Assisted Survey Estimation with Modern Prediction Techniques.  Statistical Science, 32(2):190 – 205.
  * Burrill et al., (2021) Burrill, E. A., DiTommaso, A. M., Turner, J. A., Pugh, S. A., Christensen, G., Perry, C. J., and Conkling, B. L. (2021).  The forest inventory and analysis database: Database description and user guide for phase 2 (version 9.0).
  * Cahoon and Baer, (2022) Cahoon, S. M. and Baer, K. C. (2022).  Forest resources of the tanana unit, alaska: 2018.  Technical report.
  * Chan-Golston et al., (2022) Chan-Golston, A., Banerjee, S., Belin, T. R., Roth, S. E., and Prelip, M. L. (2022).  Bayesian finite-population inference with spatially correlated measurements.  Japanese Journal of Statistics and Data Science, 5:407–430.
  * Chan-Golston et al., (2020) Chan-Golston, A. M., Banerjee, S., and Handcock, M. S. (2020).  Bayesian inference for finite populations under spatial process settings.  Environmetrics, 31(3):e2606.
  * Cook et al., (2013) Cook, B. D., Corp, L. A., Nelson, R. F., Middleton, E. M., Morton, D. C., McCorkel, J. T., Masek, J. G., Ranson, K. J., Ly, V., and Montesano, P. M. (2013).  Nasa goddard’s lidar, hyperspectral and thermal (g-liht) airborne imager.  Remote Sensing, 5(8):4045–4066.
  * Coulston et al., (2021) Coulston, J. W., Green, P. C., Radtke, P. J., Prisley, S. P., Brooks, E. B., Thomas, V. A., Wynne, R. H., and Burkhart, H. E. (2021).  Enhancing the precision of broad-scale forestland removals estimates with small area estimation techniques.  Forestry: An International Journal of Forest Research, 94(3):427–441.
  * Datta et al., (2016) Datta, A., Banerjee, S., Finley, A. O., and Gelfand, A. E. (2016).  Hierarchical nearest-neighbor gaussian process models for large geostatistical datasets.  Journal of the American Statistical Association, 111(514):800–812.  PMID: 29720777.
  * Ekström and Nilsson, (2021) Ekström, M. and Nilsson, M. (2021).  A comparison of model-assisted estimators, with and without data-driven transformations of auxiliary variables, with application to forest inventory.  Frontiers in Forests and Global Change, 4.
  * Finley et al., (2011) Finley, A. O., Banerjee, S., and MacFarlane, D. W. (2011).  A hierarchical model for quantifying forest variables over large heterogeneous landscapes with uncertain forest areas.  Journal of the American Statistical Association, 106(493):31–48.  PMID: 26139950.
  * Finley et al., (2019) Finley, A. O., Datta, A., Cook, B. D., Morton, D. C., Andersen, H. E., and Banerjee, S. (2019).  Efficient algorithms for bayesian nearest neighbor gaussian processes.  Journal of Computational and Graphical Statistics, 28(2):401–414.
  * Gelfand and Banerjee, (2010) Gelfand, A. E. and Banerjee, S. (2010).  Multivariate spatial process models.  In Gelfand, A., Diggle, P., Fuentes, M., and Guttorp, P., editors, Handbook of Spatial Statistics, pages 495–516. CRC Press, Boca Raton, FL.
  * Gelman et al., (2014) Gelman, A., Hwang, J., and Vehtari, A. (2014).  Understanding predictive information criteria for bayesian models.  Statistics and Computing, 24(6):997–1016.
  * Ghosh, (2012) Ghosh, M. (2012).  Finite population sampling: A model-design synthesis.  Statistics in Transition new series, 13(2):235–242.
  * Green et al., (2020) Green, E., Finley, A., and Strawderman, W. (2020).  Introduction to Bayesian Methods in Ecology and Natural Resources.  Springer International Publishing.
  * Gregoire, (1998) Gregoire, T. G. (1998).  Design-based and model-based inference in survey sampling: appreciating the difference.  Canadian Journal of Forest Research, 28(10):1429–1447.
  * Gregoire et al., (2011) Gregoire, T. G., Ståhl, G., Næsset, E., Gobakken, T., Nelson, R., Nelson, R., and Holm, S. (2011).  Model-assisted estimation of biomass in a lidar sample survey in hedmark county, norway.  Canadian Journal of Forest Research, 41(1):83–95.
  * Hansen et al., (2013) Hansen, M. C., Potapov, P. V., Moore, R., Hancher, M., Turubanova, S. A., Tyukavina, A., Thau, D., Stehman, S. V., Goetz, S. J., Loveland, T. R., Kommareddy, A., Egorov, A., Chini, L., Justice, C. O., and Townshend, J. R. G. (2013).  High-resolution global maps of 21st-century forest cover change.  Science, 342(6160):850–853.
  * Homer et al., (2015) Homer, C., Dewitz, J., Yang, L., Jin, S., Danielson, P., Xian, G., Coulston, J., Herold, N., Wickham, J., and Megown, K. (2015).  Completion of the 2011 national land cover database for the conterminous united states - representing a decade of land cover change information.  Photogrammetric Engineering and Remote Sensing, 81:346–354.
  * Horvitz and Thompson, (1952) Horvitz, D. G. and Thompson, D. J. (1952).  A Generalization of Sampling Without Replacement From a Finite Universe.  Journal of the American Statistical Association, 47(260):663–685.
  * Hou et al., (2021) Hou, Z., Domke, G. M., Russell, M. B., Coulston, J. W., Nelson, M. D., Xu, Q., and McRoberts, R. E. (2021).  Updating annual state-and county-level forest inventory estimates with data assimilation and FIA data.  Forest Ecology and Management, 483:118777.
  * Köhl et al., (2006) Köhl, M., Magnussen, S., and Marchetti, M. (2006).  Sampling methods, remote sensing and GIS multiresource forest inventory.  Springer.
  * Lister et al., (2020) Lister, A. J., Andersen, H., Frescino, T., Gatziolis, D., Healey, S., Heath, L. S., Liknes, G. C., McRoberts, R. E., Moisen, G. G., Nelson, M., et al. (2020).  Use of Remote Sensing Data to Improve the Efficiency of National Forest Inventories: A Case Study from the United States National Forest Inventory.  Forests, 11(12):1364.
  * Little, (2004) Little, R. J. (2004).  To model or not to model? Competing modes of inference for finite population sampling.  Journal of the American Statistical Association, 99(466):546–556.
  * Magnussen et al., (2018) Magnussen, S., Nord-Larsen, T., and Riis-Nielsen, T. (2018).  Lidar supported estimators of wood volume and aboveground biomass from the danish national forest inventory (2012–2016).  Remote Sensing of Environment, 211:146–153.
  * Mateu and Müller, (2012) Mateu, J. and Müller, W. (2012).  Spatio-temporal Design: Advances in Efficient Data Acquisition.  Statistics in Practice. Wiley.
  * May et al., (2023) May, P., McConville, K. S., Moisen, G. G., Bruening, J., and Dubayah, R. (2023).  A spatially varying model for small area estimates of biomass density across the contiguous united states.  Remote Sensing of Environment, 286:113420.
  * May et al., (2024) May, P. B., Finley, A. O., and Dubayah, R. O. (2024).  A spatial mixture model for spaceborne lidar observations over mixed forest and non-forest land types.  Journal of Agricultural, Biological and Environmental Statistics.
  * McConville et al., (2017) McConville, K. S., Breidt, F. J., Lee, T. C. M., and Moisen, G. G. (2017).  Model-Assisted Survey Regression Estimation with the Lasso.  Journal of Survey Statistics and Methodology, 5(2):131–158.
  * McGaughey et al., (2017) McGaughey, R. J., Ahmed, K., Andersen, H.-E., and Reutebuch, S. E. (2017).  Effect of occupation time on the horizontal accuracy of a mapping-grade GNSS receiver under dense forest canopy.  Photogrammetric Engineering and Remote Sensing, 83(12):861–868.
  * McRoberts, (2010) McRoberts, R. E. (2010).  Probability- and model-based approaches to inference for proportion forest using satellite imagery as ancillary data.  Remote Sensing of Environment, 114(5):1017–1025.
  * Peruzzi et al., (2021) Peruzzi, M., Banerjee, S., Dunson, D. B., and Finley, A. O. (2021).  Grid-parametrize-split (grips) for improved scalable inference in spatial big data analysis.
  * Rao and Molina, (2015) Rao, J. N. and Molina, I. (2015).  Small area estimation.  John Wiley & Sons.
  * Ren and Banerjee, (2013) Ren, Q. and Banerjee, S. (2013).  Hierarchical factor models for large spatially misaligned data: A low-rank predictive process approach.  Biometrics, 69(1):19–30.
  * Saarela et al., (2016) Saarela, S., Holm, S., Grafström, A., Schnell, S., Næsset, E., Gregoire, T. G., Nelson, R. F., and Ståhl, G. (2016).  Hierarchical model-based inference for forest inventory utilizing three sources of information.  Annals of Forest Science, 73(4):895–910.
  * Saarela et al., (2018) Saarela, S., Holm, S., Healey, S. P., Andersen, H.-E., Petersson, H., Prentius, W., Patterson, P. L., Næsset, E., Gregoire, T. G., and Ståhl, G. (2018).  Generalized hierarchical model-based estimation for aboveground biomass assessment using gedi and landsat data.  Remote Sensing, 10(11).
  * Särndal et al., (2003) Särndal, C., Swensson, B., and Wretman, J. (2003).  Model Assisted Survey Sampling.  Springer Series in Statistics. Springer New York.
  * Särndal et al., (1978) Särndal, C.-E., Thomsen, I., Hoem, J. M., Lindley, D. V., Barndorff-Nielsen, O., and Dalenius, T. (1978).  Design-based and model-based inference in survey sampling [with discussion and reply].  Scandinavian Journal of Statistics, 5(1):27–52.
  * Schroeder et al., (2014) Schroeder, T. A., Healey, S. P., Moisen, G. G., Frescino, T. S., Cohen, W. B., Huang, C., Kennedy, R. E., and Yang, Z. (2014).  Improving estimates of forest disturbance by combining observations from Landsat time series with US Forest Service Forest Inventory and Analysis data.  Remote Sensing of Environment, 154:61–73.
  * Shirota et al., (2022) Shirota, S., Finley, A. O., Cook, B. D., and Banerjee, S. (2022).  Conjugate sparse plus low rank models for efficient bayesian interpolation of large spatial data.  Environmetrics, page e2748.
  * Smith, (2002) Smith, W. B. (2002).  Forest inventory and analysis: a national inventory and monitoring program.  Environmental pollution, 116:S233–S242.
  * Spiegelhalter et al., (2002) Spiegelhalter, D. J., Best, N. G., Carlin, B. P., and Van Der Linde, A. (2002).  Bayesian measures of model complexity and fit.  Journal of the Royal Statistical Society Series B, 64(4):583–639.
  * (53) Ståhl, G., Holm, S., Gregoire, T. G., Gobakken, T., Næsset, E., and Nelson, R. (2011a).  Model-based inference for biomass estimation in a lidar sample survey in hedmark county, norwaythis article is one of a selection of papers from extending forest inventory and monitoring over space and time.  Canadian Journal of Forest Research, 41(1):96–107.
  * (54) Ståhl, G., Holm, S., Gregoire, T. G., Gobakken, T., Næsset, E., and Nelson, R. (2011b).  Model-based inference for biomass estimation in a lidar sample survey in hedmark county, norwaythis article is one of a selection of papers from extending forest inventory and monitoring over space and time.  Canadian Journal of Forest Research, 41(1):96–107.
  * Ståhl et al., (2016) Ståhl, G., Saarela, S., Schnell, S., Holm, S., Breidenbach, J., Healey, S. P., Patterson, P. L., Magnussen, S., Næsset, E., McRoberts, R. E., and Gregoire, T. G. (2016).  Use of models in large-area forest surveys: comparing model-assisted, model-based and hybrid estimation.  Forest Ecosystems, 3(1):5.
  * Strunk et al., (2019) Strunk, J., Packalen, P., Gould, P., Gatziolis, D., Maki, C., Andersen, H.-E., and McGaughey, R. J. (2019).  Large area forest yield estimation with pushbroom digital aerial photogrammetry.  Forests, 10(5):397.
  * Taylor-Rodriguez et al., (2019) Taylor-Rodriguez, D., Finley, A. O., Datta, A., Babcock, C., Andersen, H.-E., Cook, B. D., Morton, D. C., and Banerjee, S. (2019).  Spatial factor models for high-dimensional and large spatial data: An application in forest variable mapping.  Statistica Sinica, 29(3):1155–1180.
  * Vehtari et al., (2017) Vehtari, A., Gelman, A., and Gabry, J. (2017).  Practical bayesian model evaluation using leave-one-out cross-validation and waic.  Statistics and Computing, 27(5):1413–1432.
  * Watanabe, (2010) Watanabe, S. (2010).  Asymptotic equivalence of bayes cross validation and widely applicable information criterion in singular learning theory.  J. Mach. Learn. Res., 11:3571–3594.
  * Westfall et al., (2011) Westfall, J. A., Patterson, P. L., and Coulston, J. W. (2011).  Post-stratified estimation: within-strata and total sample size recommendations.  Canadian Journal of Forest Research, 41(5):1130–1139.
  * Wurtzebach et al., (2020) Wurtzebach, Z., DeRose, R. J., Bush, R. R., Goeking, S. A., Healey, S., Menlove, J., Pelz, K. A., Schultz, C., Shaw, J. D., and Witt, C. (2020).  Supporting national forest system planning with forest inventory and analysis data.  Journal of Forestry, 118(3):289–306.
  * Xia et al., (2006) Xia, G., Miranda, M. L., and Gelfand, A. E. (2006).  Approximately optimal spatial design approaches for environmental health data.  Environmetrics, 17(4):363–385.
  * Zhang and Banerjee, (2022) Zhang, L. and Banerjee, S. (2022).  Spatial factor modeling: A Bayesian matrix-normal approach for misaligned data.  Biometrics, 78(2):560–573.
  * Zhang, (2016) Zhang, X. (2016).  An optimized blas library based on gotoblas2.  https://github.com/xianyi/OpenBLAS/.  Accessed 2015-06-01.