Bayesian Inference for Finite Populations Under Spatial Process Settings

TL;DR

This paper introduces a Bayesian spatial process model for finite population estimation that improves inference accuracy and spatial mapping by accounting for geographic dependence in sampling designs.

Contribution

It presents a novel framework integrating spatial dependence into finite population Bayesian inference, enhancing model fit and spatial interpolation capabilities.

Findings

Process-based models outperform non-spatial models in fit and inference.

Spatial models provide richer, interpolated maps of the population.

Application to California groundwater demonstrates scalable, improved estimation.

Abstract

We develop a Bayesian model-based approach to finite population estimation accounting for spatial dependence. Our innovation here is a framework that achieves inference for finite population quantities in spatial process settings. A key distinction from the small area estimation setting is that we analyze finite populations referenced by their geographic coordinates. Specifically, we consider a two-stage sampling design in which the primary units are geographic regions, the secondary units are point-referenced locations, and the measured values are assumed to be a partial realization of a spatial process. Estimation of finite population quantities from geostatistical models do not account for sampling designs, which can impair inferential performance, while design-based estimates ignore the spatial dependence in the finite population. We demonstrate using simulation experiments that process-based finite population sampling models improve model fit and inference over models that fail to account for spatial correlation. Furthermore, the process based models offer richer inference with spatially interpolated maps over the entire region. We reinforce these improvements and demonstrate scaleable inference for ground-water nitrate levels in the population of California Central Valley wells by offering estimates of mean nitrate levels and their spatially interpolated maps.

Figures

Tables

Table 1. Table 1: Comparison of Parameter Estimation and Model Fit in Two Exact Models

	Two-Stage	Spatial
$\nu$ (2)	2.60 (1.55, 3.58)	2.79 (1.36, 4.23)
${\delta }^2$ (4)	6.36 (5.47, 7.45)	3.84 (3.35, 4.43)
${\tau }^2$ (9)	—	8.65 (7.53, 9.96)
$\overline{y}$ (2.60)	2.66 (7.53, 9.96)	2.92 (7.53, 9.96)
WAIC	1803.83 (25.52)	1686.7 (27.6)
D = G + P	66100.83=34203.26+31897.56	45335.54=22875.45+22460.08

Table 2. Table 2: Results of Data Analysis

Model	FP Mean (95% CI)	WAIC (SE)
1. Two-Stage	26.6 (19.2, 35.1)	4701.6 (78.6)
2. Spatial	32.6 (27.2, 38.7)	4858.4 (78.0)
3. Two-Stage + Spatial	31.2 (24.2, 37.3)	2697.0 (150.5)
4. Regional Spatial	26.2 (18.7, 34.4)	4536.6 (93.1)

Table 3. Table 3: Comparison of Estimates and Model Fit from MCMC

Model	FP Mean (95% CI)	$\nu$ (95% CI)	WAIC (SE)
1. Two-Stage	2.84 (1.93, 3.67)	2.84 (1.82, 3.82)	1870.26 (35.00)
2. Spatial	3.04 (1.33, 4.56)	3.24 (1.20, 5.01)	1912.70 (26.36)
3. Two-Stage + Spatial	3.11 (-0.86, 6.70)	3.22 (-1.23, 7.33)	455.67 (17.03)
4. Regional Spatial	2.56 (1.74, 3.37)	2.44 (1.49, 3.47)	1202.08 (38.99)

Equations

$$[\text{mbox}data|\text{mbox}process]\times [\text{mbox}process|\text{mbox}parameters]\times [\text{mbox}parameters],$$

$$\begin{bmatrix}y_{s}\\ y_{ns}\end{bmatrix}=\begin{bmatrix}X_{s}\\ X_{ns}\end{bmatrix}\beta+\begin{bmatrix}\epsilon_{s}\\ \epsilon_{ns}\end{bmatrix}\;;\\ \begin{bmatrix}\epsilon_{s}\\ \epsilon_{ns}\end{bmatrix}\sim N\Bigg{(}\begin{bmatrix}0\\ 0\end{bmatrix},\begin{bmatrix}V_{s}(\theta)&V_{s,ns}(\theta)\\ V_{ns,s}(\theta)&V_{ns}(\theta)\end{bmatrix}\Bigg{)}\;.$$

$$\text{E}[{\alpha }^{\top }y|y_s]=\{{\alpha }_s^{\top }+{\alpha }_{ns}^{\top }[V_{ns,s}+{\alpha }_{ns}^{\top }Q{V}_{\beta |y_s}{X}_s^{\top }]V_s^{-1}\}{y}_s.$$

$$\text{Var}[\text{E}[{\alpha }^{\top }y|y_s]]={\alpha }_s^{\top }{V}_s{\alpha }_s+2{\alpha }_s^{\top }{B}_V^{\top }{\alpha }_{ns}+{\alpha }_{ns}^{\top }{B}_V{V}_s^{-1}{B}_V^{\top }{\alpha }_{ns}.$$

$$\text{Var}[{\alpha }^{\top }y|y_s]={\alpha }_{ns}^{\top }(Q{V}_{\beta |y_s}{Q}^{\top }+V_{ns}-V_{ns,s}{V}_s^{-1}{V}_{s,ns}){\alpha }_{ns}.$$

$$\text{E}[\alpha^{\top}y\,|\,y_{s}]=\sum_{i=1}^{n}\Bigg{(}\alpha_{i}+\frac{\sum_{i=n+1}^{N}\alpha_{i}}{\frac{\sigma^{2}}{\xi^{2}}+n}\Bigg{)}y_{i}$$

$$\text{E}[\alpha^{\top}y\,|\,y_{s}]=\sum_{i=1}^{n}\sum_{j=1}^{m_{i}}\Bigg{(}\alpha_{ij}+\Bigg{[}\alpha_{i}+\frac{\sum_{i=1}^{N}\alpha_{i}(1-\lambda_{i})}{1/\tilde{\gamma}^{2}+\sum_{i=1}^{n}\lambda_{i}}\Bigg{]}\frac{\lambda_{i}}{m_{i}}\Bigg{)}y_{ij}\\$$

$$y({\ell }_{ij})=\mu ({\ell }_{ij})+\omega ({\ell }_{ij})+ϵ({\ell }_{ij});\omega \sim N(0,\Omega );ϵ\sim N(0,V^{(\sigma )}),$$

$$\begin{split}\text{E}[\alpha^{\top}y\,|\,y_{s}]&=\sum_{i=1}^{n}\sum_{j=1}^{m_{i}}\Bigg{(}\alpha_{ij}+\alpha_{ns}^{\top}\Bigg{\{}\Omega_{ns,s}Q_{s}+\;[X_{ns}-\Omega_{ns,s}Q_{s}^{-1}X_{s}]\times\begin{pmatrix}\frac{1}{\delta^{2}}I_{N}+X_{s}^{\top}Q_{s}^{-1}X_{s}\end{pmatrix}^{-1}\\ &\quad\;\times\Bigg{[}X_{s}^{\top}Q_{s}^{-1}+\frac{\frac{1}{\delta^{2}}1_{N}1_{N}^{\top}X_{s}^{\top}(\delta^{2}X_{s}X_{s}^{\top}+Q_{s})^{-1}}{\frac{1}{\gamma^{2}}+\sum_{i=1}^{n}\lambda_{i}^{*}}\Bigg{]}\Bigg{\}}q_{ij}\Bigg{)}y(\ell_{ij})\;,\end{split}$$

$$p(\theta ,\nu ,\beta |y_s)\propto p(\theta )\times N(\nu |0,V_{\nu })\times N(\beta |A\nu ,V_{\beta })\times N(y_s|X_s\beta ,V_s(\theta )).$$

$$p(y_{ns}|y_s)=\int p(y_{ns}|y_s,\theta ,\nu ,\beta )\times p(\theta ,\nu ,\beta |y_s)d\theta d\nu d\beta .$$

$$IG({\delta }^2|a,b)\times N(\nu |0,{\delta }^2{\tilde{V}}_{\nu })\times N(\beta |A\nu ,{\delta }^2{\tilde{V}}_{\beta })\times N(y_s|X\beta ,{\delta }^2{\tilde{V}}_s).$$

$$\nu |y_s,{\delta }^2\sim N(\nu |c,{\delta }^{2}d),$$

$$\delta^{2}\,|\,y_{s}\sim IG\Bigg{(}a+\frac{k}{2},b+\frac{1}{2}\Big{[}y_{s}^{\top}(X_{s1}X_{s1}^{\top}+\tilde{V}_{s}^{(\sigma)})^{-1}y_{s}+\frac{c^{2}}{d}\Big{]}\Bigg{)}\;.$$

$$\underset{{\tilde{\gamma }}^2\to \infty }{\lim }c=\frac{{\sum }_{i=1}^n{\lambda }_i{\bar{y}}_i}{{\sum }_{i=1}^n{\lambda }_i}\text{and}\underset{{\tilde{\gamma }}^2\to \infty }{\lim }{\delta }^{2}d=\frac{{\delta }^2}{{\sum }_{i=1}^n{\lambda }_i}.$$

$${\mu }_s|y_s,\nu ,{\delta }^2\sim N({\mu }_s|c_{\ast },{\delta }^2{d}_{\ast }),$$

$$\delta^{2}\,|\,y_{s}\sim IG\left[a+\frac{k}{2},b+\frac{1}{2}y_{s}^{\top}\Big{(}\tilde{\Omega}_{s}^{-1}-\tilde{\Omega}_{s}^{-1}1_{k}V_{\tilde{\Omega}_{s}}1_{k}^{\top}\tilde{\Omega}_{s}^{-1}\Big{)}y_{s}\right]\text{and}$$

$$\mu_{s}\,|\,y_{s},\delta^{2}\sim N\Bigg{[}V_{\tilde{\Omega}_{s}}1_{k}^{\top}\tilde{\Omega}_{s}^{-1}y_{s},\delta^{2}V_{\tilde{\Omega}_{s}}\Bigg{]}\;.$$

$$a^{\ast }$$

$$b^{\ast }$$

$$M_{\nu }$$

$$m_{\nu }$$

$$M_{\beta }$$

$$m_{\beta }$$

$$b$$

$$B^{-1}$$

$$b_{\ast }$$

$$B_{\ast }^{-1}$$

$$B_{\ast }{b}_{\ast }$$

$$={\left[\begin{array}{c}(1-{\lambda }_1)\nu +{\lambda }_1{\bar{y}}_1,\cdots ,(1-{\lambda }_n)\nu +{\lambda }_n{\bar{y}}_n\end{array}\right]}^{\top }.$$

$$\text{E}[{\alpha }^{\top }y|y_s]$$

$$={\alpha }_s^{\top }{y}_s+{\alpha }_{ns}^{\top }\text{E}[X_{ns}\beta +V_{ns,s}{V}_s^{-1}(y_s-X_s\beta )|y_s]$$

$$={\alpha }_s^{\top }{y}_s+{\alpha }_{ns}^{\top }{V}_{ns,s}{V}_s^{-1}{y}_s+Q{V}_{\beta |y_s}{X}_s^{\top }{V}_s^{-1}{y}_s$$

$$=\{{\alpha }_s^{\top }+{\alpha }_{ns}^{\top }[V_{ns,s}+{\alpha }_{ns}^{\top }Q{V}_{\beta |y_s}{X}_s^{\top }]V_s^{-1}\}{y}_s,$$

$$\text{Var}[\text{E}[{\alpha }^{\top }y|y_s]]$$

$$={\alpha }_s^{\top }{V}_s{\alpha }_s+2{\alpha }_{ns}^{\top }{B}_V{\alpha }_s+{\alpha }_{ns}^{\top }{B}_V{V}_s^{-1}{B}_V^{\top }{\alpha }_{ns}\text{, and}$$

$$\text{Var}[\text{E}[{\alpha }^{\top }y|y_s]]$$

$$=\text{Var}[\text{E}[{\alpha }_{ns}^{\top }{y}_{ns}|\beta ,y_s]|y_s]+\text{E}[\text{Var}[{\alpha }_{ns}^{\top }{y}_{ns}|\beta ,y_s]|y_s]$$

Full text

Key words: Finite population inference; Bayesian modeling; Spatial process; Two-stage sampling; Hierarchical models.

I.   Introduction

Finite population survey sampling concerns statistical modeling and inference on finite populations from sampling designs; see, for example, Cochran (1977), Hartley and Sielken Jr. (1975), Royall (1970), and Horvitz and Thompson (1952). In this article, we will concern ourselves with Bayesian inference for finite populations when the sampling units are spatially oriented. For instance, one may consider estimating the total biomass in a forest given a sample of trees, the average income of a city given a sample of individuals and their addresses, or the total amount of air pollution attributable to cars on a freeway given a sample of pollution measurements. Additionally, these techniques could be used to conduct Bayesian inference for finite population survey sampling is discussed in great detail in Gelman (2007), Little (2004), Ghosh and Meeden (1997), and Ericson (1969). In this domain, there is a substantial literature on small area estimation for regionally aggregated data (see, e.g., Rao, 2003; Ghosh et al., 1998; Ghosh and Rao, 1994; Clayton and Kaldor, 1987), where interest lies in modeling dependencies across regions.

Unlike the aforementioned literature on small area estimation, where the sampling units are regions such as counties, states or census-tracts, spatial process models consider quantities that, at least conceptually, exist in continuum over the entire domain. The process assigns a probability law to an uncountable subset within a $d$-dimensional Euclidean domain. In general, spatial process modeling (Banerjee et al. 2014; Cressie and Wikle 2011; and Ripley 2004) follows the generic paradigm

[TABLE]

which accommodates complex dependencies and multiple sources of variation.

With regard to finite population sampling in spatial process settings, the literature appears to be considerably more scant than small area estimation. Here, Ver Hoef (2002) discuss connections between geostatistical models and classical design-based sampling and develop methods for executing model-based block kriging. Cicchitelli and Montanari (2012) present a spline regression model-assisted, design-based estimator of the mean for use on a random sample from both finite and infinite spatial populations. A linear spatial interpolator is used by Bruno et al. (2013) to create a design-based predictor of values at unobserved locations which outperforms non-spatial predictors. While related to these developments, none of these techniques have presented a Bayesian approach. In this manuscript, we pursue a fully Bayesian, model-based approach as in (1) and carry out inference on the finite population quantities and the spatial process.

Bayesian finite population survey sampling is essentially model-based (see, e.g., Little 2004). The population units are themselves assumed to be endowed with a probability distribution. In a Gaussian setting, Scott and Smith (1969) devised Bayesian hierarchical models for inferring with two-stage designs, while Malec and Sedransk (1985) extended this framework to general multi-stage (more than two-stages) models and also discussed handling unknown variances. Our current contribution focuses on incorporating survey sampling designs within (1). We extend this framework to spatial process settings under the context of ignorable sampling designs (Rubin 1976; Sugden and Smith 1984), where the probability of element selection is assumed independent of the measured outcome given the design variables. We specifically develop the distribution theory and algorithms for implementing (1) in the context of two-stage designs that encompass simple random, cluster and stratified sampling (as defined in Cochran, 1977) as special cases. Extension of this work to multi-stage present no new methodological difficulties, building upon Malec and Sedransk (1985).

The remainder of this paper evolves as follows. In Section II, we review a general framework for Bayesian modeling for multi-stage sampling and how simple, two-stage, and stratified random sampling designs arise as special cases. Section III presents modeling strategies for spatially correlated data sampled using a two-stage design, the implementation of which, along with the model proposed by Scott and Smith (1969), is discussed in Section IV using Bayesian exact and Markov chain Monte Carlo (MCMC) sampling. These models are then applied to simulated data in Section V and then used in an analysis of nitrate levels in California groundwater in Section VI. The paper concludes with a brief discussion of the results in Section VII.

II.   Bayesian modeling of multi-stage sampling

Suppose that $n$ samples are randomly drawn from a finite population of size $N$, $n\leq N$ and for the $i$-th sampled unit, the outcome $y_{i}$ is measured. Without loss of generality, let the set of outcomes from the finite population be stacked into a vector $y=[y_{s}^{\top}\colon y_{ns}^{\top}]^{\top}$, where $y_{s}=[y_{1},\dots,y_{n}]^{\top}$ and $y_{ns}=[y_{n+1},\dots,y_{N}]^{\top}$ are vectors of outcome values from the sampled and nonsampled units, respectively. This vector has corresponding design matrix $X=[X_{s}^{\top}\colon X_{ns}^{\top}]^{\top}$, which observed for the entire finite population and denotes group membership.

From a superpopulation perspective, we consider the finite population to be a random sample of size $N$ from an infinitely large population. This superpopulation is assumed to follow a Gaussian distribution with mean $\nu$ and a covariance function defined by parameters $\theta$. In general, we can construct the following linear regression model:

[TABLE]

Bayesian specifications further model $\beta\sim N(A\nu,V_{\beta})$, where $A$ and $\nu$ are a vector with length equal to the number of groups and a scalar, respectively, and $V_{\beta}$ is the variance of $\beta$. The hierarchy continues with probabilistic specifications on $\nu$ and $\theta$.

Our goal is to estimate linear finite population quantities of the form $\alpha^{\top}y$, where $\alpha$ is a given, fixed vector of weights defined for the entire population. Suppose $\nu\sim\text{N}(0,\gamma^{2})$. Define $V_{\beta\,|\,y_{s}}=[(\gamma^{2}AA^{\top}+V_{\beta})^{-1}+X_{s}^{\top}V_{s}^{-1}X_{s}]^{-1}$ and $Q=X_{ns}-V_{ns,s}V_{s}^{-1}X_{s}$. Fixing the variance parameters, the posterior expectation of the finite population quantity is:

[TABLE]

Defining $B_{V}=V_{ns,s}+QV_{\beta\,|\,y_{s}}X_{s}^{\top}$, the variance of this expectation is:

[TABLE]

Additionally, the posterior variance of $\alpha^{\top}y$ is:

[TABLE]

For the special case of a census, in which all members of the population are sampled, e.g. $y_{s}=y$, the conditional expectation of the finite population quantity is finite population consistent, in the sense that $\text{E}[\alpha^{\top}y\,|\,y_{s}]=\text{E}[\alpha^{\top}y\,|\,y]=\alpha^{\top}y$. Different sampling designs can incorporated by appropriately structuring the sampled and nonsampled elements. We provide a few examples below. All derivations can be located in the supplementary materials.

Example 1**.**

Simple Random Sampling

In simple random sampling, $n$ units are randomly drawn from a population of size $N$, where each unit in the population is independent and identically distributed with mean $\mu$ and variance $\sigma^{2}$. To express this as in (2), define $y_{s}=[y_{1},\dots,y_{n}]^{\top}$ and $y_{ns}=[y_{n+1},\dots,y_{N}]^{\top}$, with corresponding design matrices $X_{s}=1_{n}$ and $X_{ns}=1_{N-n}$, respectively, where $1_{n}$ represents the $n\times 1$ vector of ones. Let $A=1$ and take $\beta$ to be a scalar $\mu$ with mean $\nu=0$ and variance $V_{\beta}=\xi^{2}$. Additionally, let $V_{s}=\sigma^{2}I_{n}$, $V_{ns}=\sigma^{2}I_{N-n}$, and $V_{s,ns}=V_{ns,s}^{\top}=O$, where $O$ is a matrix of zeroes of appropriate order. Define finite population weights $\alpha=[\alpha_{1},\dots,\alpha_{N}]^{\top}$. Fixing the variance parameters, the posterior expectation of $\alpha^{\top}y$ is:

[TABLE]

$\square$

Example 2**.**

Two-Stage Sampling

In a more complex case, suppose that the population is divided into $N$ distinct groups defined by geography or other characteristics, with the $i$-th group of size $M_{i}$. Assume that within the $i$-th group, each unit is independent and identically distributed with mean $\mu_{i}$ and variance $\sigma_{i}^{2}$. The group means $\mu_{1},\dots,\mu_{N}$ are independent and follow a normal distribution centered at $\nu$ with a variance of $\delta^{2}$, hence $\beta=\mu=[\mu_{1},\dots,\mu_{N}]^{\top}$, $A=1_{N}$ and $V_{\beta}=\delta^{2}I_{N}$. Suppose that only $n$ of the $N$ groups are randomly sampled, where $n\leq N$. Without loss of generality, take the first $n$ groups to be sampled, and then within the chosen $i$-th group, $m_{i}$ units are randomly selected; $m_{i}\leq M_{i}$, $i=1,\dots,n$. As $N-n$ groups are not sampled, the number of observed units in these groups are zero, e.g. $m_{i}=0$, $i=n+1,\dots,N$. Hence, we can define the number of sampled units as $k=\sum_{i=1}^{N}m_{i}=\sum_{i=1}^{n}m_{i}$, the number of unsampled units as $K=\sum_{i=1}^{N}(M_{i}-m_{i})$, and the population total to be $T=K+k=\sum_{i=1}^{N}M_{i}$.

To examine (2) in the context of a two-stage design, define outcome vectors $y_{s}=[y_{1}^{\top},\dots,y_{n}^{\top}]^{\top}$ and $y_{ns}=[y_{1}^{*\top},\dots,y_{N}^{*\top}]^{\top}$, with $i^{\text{th}}$ components $y_{i}=[y_{i1},\dots,y_{im_{i}}]^{\top}$ and $y_{i}^{*}=[y_{im_{i}+1},\dots,y_{iM_{i}}]^{\top}$, respectively. The design matrix for the sampled units can be modified by fixing the $k\times N$ matrix $X_{s}=[\oplus_{i=1}^{n}1_{m_{i}}\colon O]$, reflecting that $N-n$ of the $N$ sites are unobserved. Similarly, for the unobserved units, define $X_{ns}$ as a block diagonal, $K\times N$ matrix with upper block $[\oplus_{i=1}^{n}1_{M_{i}-m_{i}}]$ and lower block $[\oplus_{i=n+1}^{N}1_{M_{i}-m_{i}}]$. For notational convenience, we also define $X_{s1}=[\oplus_{i=1}^{n}1_{m_{i}}]$ and divide the group mean vector $\mu$ into sampled, $\mu_{s}=[\mu_{1},\dots,\mu_{n}]^{\top}$, and nonsampled, $\mu_{ns}=[\mu_{n+1},\dots,\mu_{N}]^{\top}$, components such that $\mu=[\mu_{s}^{\top},\mu_{ns}^{\top}]^{\top}$. Note that distributional mean of $y_{s}$, $X_{s}1_{N}\mu$, can be simplified to $X_{s1}1_{n}\mu_{s}$, as the mean of the sampled units does not depend on $\mu_{ns}$. Define the sampled and nonsampled covariance matrices to be $V_{s}=V_{s}^{(\sigma)}=[\oplus_{i=1}^{n}\sigma_{i}^{2}I_{m_{i}}]$ and $V_{ns}=V_{ns}^{(\sigma)}=[\oplus_{i=1}^{N}\sigma_{i}^{2}I_{M_{i}-m_{i}}]$, respectively, and set $V_{s,ns}=V_{ns,s}^{\top}=O$. Additionally, define $V^{(\sigma)}=\begin{bmatrix}V_{s}^{(\sigma)}&O\\ O&V_{ns}^{(\sigma)}\end{bmatrix}$.

To make this model fully Bayesian, let $\nu\sim N(0,\gamma^{2})$ and $\delta^{2}\sim IG(a,b)$. As our interest lies in estimating $\alpha^{T}y$, we can derive the posterior distributions of $p(\delta^{2}\,|\,y_{s})$ and $p(\nu\,|\,y_{s})$ for exact sampling of the superpopulation parameters, the details of which are provided in Section II. This approach yields results similar to those derived by Scott and Smith (1969), but has the added strength of including a prior distribution on $\nu$, and Ghosh and Meeden (1997), who replaced distributional assumptions with the assumption of posterior linearity and fixed the variance parameters.

In the two-stage case, for a set of weights $\alpha=[\alpha_{11},\dots,\alpha_{NM_{N}}]^{\top}$, define the group mean of sampled units as $\bar{y}_{i}=\frac{1}{m_{i}}\sum_{j=1}^{m_{i}}y_{ij}$, $i=1,\dots,n$, and the group weight of nonsampled units as $\alpha_{i}=\sum_{j=m_{i}+1}^{M_{i}}\alpha_{ij}$, $i=1,\dots,N$. Also, let $\tilde{\gamma}^{2}=\gamma^{2}/\delta^{2}$ and define $\lambda_{i}=\delta^{2}/(\delta^{2}+\sigma_{i}^{2}/m_{i})$ if $i\in\{1,\dots,n\}$ and $\lambda_{i}=0$ if $i\in\{n+1,\dots,N\}$. Fixing all variance parameters, the expected value of the finite population estimate is

[TABLE]

Additionally, the two-stage case can be extended to a three-stage case by assuming that the $j^{\text{th}}$ element of the $i^{\text{th}}$ group has $m_{ij}^{*}$ subelements. Malec and Sedransk (1985) derive posterior distributions for the means for a three-stage sampling scheme and provide a framework to extend this to data with $t$ stages of sampling. $\square$

Example 3**.**

Stratified Random Sampling

Stratified sampling is a special case of two-stage sampling where all groups are sampled (e.g. $n=N$ and $m_{i}>0$, $i=1,\dots,N$), and, therefore, considering the same population described in Example 2, the number of sampled units is $k=\sum_{i=1}^{N}m_{i}$, the number of nonsampled units is $K=\sum_{i=1}^{N}(M_{i}-m_{i})$, and the population total is again $T=K+k=\sum_{i=1}^{N}M_{i}$. Thus, to express this design as (2), let $y_{s}=[y_{1}^{\top},\ldots,y_{N}^{\top}]^{\top}$ and $y_{ns}=[y_{1}^{*\top},\ldots,y_{N}^{*\top}]^{\top}$, with $i^{\text{th}}$ components, $y_{i}=[y_{i1},\dots,y_{im_{i}}]^{\top}$ and $y_{i}^{*}=[y_{im_{i}+1},\dots,y_{iM_{i}}]^{\top}$, respectively. To reflect a membership to one of $N$ groups, we take $X_{s}=[\oplus_{i=1}^{N}1_{m_{i}}]$, $X_{ns}=[\oplus_{i=1}^{N}1_{M_{i}-m_{i}}]$, $\beta=\mu=[\mu_{1},\dots,\mu_{N}]^{\top}$, and $A=1_{N}$. The variance components also reflect this and are defined as $V_{\beta}=\delta^{2}I_{N}$, $V_{s}=[\oplus_{i=1}^{N}\sigma_{i}^{2}I_{m_{i}}]$, $V_{ns}=[\oplus_{i=1}^{N}\sigma_{i}^{2}I_{M_{i}-m_{i}}]$, and $V_{s,ns}=V_{ns,s}^{\top}=O$.

The posterior expectation of $\alpha^{\top}y$ is given by (4), noting that $n=N$ and $\lambda_{i}=\frac{\delta^{2}}{\delta^{2}+\sigma_{i}^{2}/m_{i}}$, $i=1,...,N$, is well-defined as $m_{i}>0$ for all $i$. In fact, if non-informative priors are taken for the means, e.g. $\gamma^{2}\to\infty$ and $\delta^{2}\to\infty$, then $\lambda_{i}\to 1$, $i=1,\dots,n$ and the stratified finite population mean is $\displaystyle\text{E}\left[\frac{1}{T}1_{T}^{\top}y\,|\,y_{s}\right]=\sum_{i=1}^{N}\frac{M_{i}}{T}\bar{y}_{i}$, (see, e.g. Little, 2004). $\square$

III.   Bayesian spatial process modeling for multi-stage sampling

Data believed to be correlated as a function of geographic distance is typically described using a spatial process model. The data is assumed to be a partial realization of a Gaussian process with dependencies between elements defined by an isotropic covariance function, $C(d)$, where $d$ is the distance between any two points. Several choices for $C(d)$ are available (see, e.g. Banerjee et al., 2014), but a versatile family is the Matérn, defined as $C(d_{ab})=\sigma^{2}+\tau^{2}$ if $d_{ab}=0$ and C(d_{ab})=\tau^{2}\frac{2^{1-\eta}}{\Gamma(\eta)}\big{(}\sqrt{2\eta}d_{ab}\phi\big{)}^{\eta}K_{\eta}\big{(}\sqrt{2\eta}d_{ab}\phi\big{)} if $d_{ab}>0$, where $K_{\eta}(\cdot)$ is the modified Bessel function, $d_{ab}$ is the distance between two locations $\ell_{a}$ and $\ell_{b}$. Here $\sigma^{2}$ captures variation due to measurement error or micro-resolution variation, $\tau^{2}$ is the spatial variance, $\phi$ is a decay parameter which determines the rate of decline in spatial association, and $\eta$ is a smoothness parameter. The exponential covariance function is a special case of Matérn when $\eta=1/2$. In this specific instance, the decay parameter is used to calculate the effective spatial range, which is the distance where spatial correlation between two points drops below 0.05.

Extending Example 2 to a geographic context, our spatial domain comprises $N$ regions. Let $\ell_{ij}$ denote the $j$-th location in region $i$. The finite population is described by values $y(\ell_{ij})$, $i=1,\ldots,N$ and $j=1,\ldots,M_{i}$. Let $y_{s}$ be the $k\times 1$ vector corresponding to measurements from the sampled locations and $y_{ns}$ be the $K\times 1$ vector of unsampled measurements. Consider the following spatial regression model for the two-stage finite population,

[TABLE]

where $\mu(\ell_{ij})$ is the mean of the outcome at $\ell_{ij}$, $\omega=[\omega_{s}^{\top}:\omega_{ns}^{\top}]^{\top}$ and $\epsilon=[\epsilon_{s}^{\top}:\epsilon_{ns}^{\top}]^{\top}$ are $T\times 1$ vectors formed by stacking up $\omega(\ell_{ij})$’s and $\epsilon(\ell_{ij})$’s, respectively (analogous to $y$ in Example 2). Here, $\omega$ accounts for spatial effects and $\Omega$ is the $T\times T$ spatial covariance matrix constructed with $C(d_{ab})$ and is partitioned as $\Omega=\begin{bmatrix}\Omega_{s}&\Omega_{s,ns}\\ \Omega_{ns,s}&\Omega_{ns}\end{bmatrix}$. Introducing spatial effects in Example 2 yields $\mu(\ell_{ij})=\mu_{i}$, $V_{s}=\Omega_{s}+V_{s}^{(\sigma)}$, $V_{ns}=\Omega_{ns}+V_{ns}^{(\sigma)}$, and $V_{s,ns}=V_{ns,s}^{\top}=\Omega_{s,ns}$ in (2). This also accommodates spatial versions of Examples1 and 3 by setting $N=1$ and $N=n$, respectively.

Analogous to (4), the posterior estimate of a linear function of the population values is

[TABLE]

where $\lambda^{*\top}=[\lambda_{1}^{*}\dots,\lambda_{N}^{*}]=1_{N}^{\top}X_{s}^{\top}(\delta^{2}X_{s}X_{s}^{\top}+\Omega_{s}+V_{s}^{(\sigma)})^{-1}X_{s}$, $Q_{s}=\Omega_{s}+V_{s}^{(\sigma)}$, and $q_{ij}$ is a set of $k$ indicator vectors of length $k$, $i=1,\dots,n$, $j=1,\dots,m_{i}$. For $i=1$, $q_{1j}=j$ and 0 elsewhere, and if $i>1$, $q_{ij}$ is 1 at element $\sum_{i=1}^{i-1}m_{i}+j$ and 0 elsewhere. This two-stage spatial model, (5), can be written as an intercept-only spatial model by setting $\mu(\ell_{ij})=\mu$ and $\sigma_{i}^{2}=\sigma^{2}$, $i=1,\dots,N$, i.e., simplifying $V^{(\sigma)}$ to $\sigma^{2}I_{T}$. As region is not accounted for, the design matrices $X_{s}$ and $X_{ns}$ are replaced with $1_{k}$ and $1_{K}$, respectively.

However, as the size of the finite population, $T$, grows, the scaleability of (5) diminishes due to an increased computational burden stemming from the inversion of the $T\times T$ matrix $\Omega$. To address this, we also consider a more computationally efficient model which also allows for region specific means, but specifies that each region is defined by its own process parameters and is independent from all other regions. To reflect this regional independence, we specify the covariance function specifying the spatial process $\omega(\ell)$ in (5) to be [math] for any two points in different regions, and equal to the value of the Matérn covariance function for any two points within the same region.

Comparing the finite populations estimates given in (4) and (6), it is evident that accounting for spatial variation results in a more complex equation, as all observed and unobserved outcome values in a population can no longer be assumed to be independent conditional on the group means. This can also be seen in the calculation of the $\lambda$ parameters, which in the two-stage model, are a simple ratio of variances. In the spatial case, however, the complexity of the parameters is increased by the addition of the spatial covariance matrix.

IV.   Model Implementation and Assessment

I.   General framework

A Bayesian linear model corresponding to the likelihood of the sampled data in (2) is

[TABLE]

We use Markov chain Monte Carlo algorithms (see, e.g. Robert and Casella, 2004) for sampling from (7). Subsequeny Bayesian inference for $y_{ns}$ is available in posterior predictive fashion by drawing samples from

[TABLE]

Using the conditional independence of parameters in (2), we obtain $p(y_{ns}\,|\,y_{s},\theta,\nu,\beta)=N(y_{ns}\,|\,\mu_{ns\,|\,s},V_{ns\,|\,s})$, where $\mu_{ns\,|\,s}=X_{ns}\beta+V_{ns,s}(\theta)V_{s}(\theta)^{-1}(y_{s}-X_{s}\beta)$ and $V_{ns\,|\,s}=V_{ns}(\theta)-V_{ns,s}(\theta)V_{s}(\theta)^{-1}V_{s,ns}(\theta)$. Therefore, sampling from (8) is achieved by drawing one $\{\beta,\theta\}$ from (7) followed by one $y_{ns}\sim N(X_{ns}\beta,V_{ns}(\theta))$, for each posterior sample of $\{\beta,\theta\}$. The resulting samples provide inference on the nonsampled group means $\mu_{ns}$ and Bayesian imputation for the nonsampled population units, $y_{ns}$.

These samples from the posterior predictive distribution can be used to obtain posterior finite population estimates of the form $\alpha^{\top}y$. We consider four models using (7).

Model 1. Two-Stage

For the model provided in Example 2, we take $\theta=[\gamma^{2},\delta^{2},\sigma_{1}^{2},\dots,\sigma_{N}^{2}]^{\top}$, $p(\theta)=IG(\gamma^{2}\,|\,a_{\gamma},b_{\gamma})\times IG(\delta^{2}\,|\,a_{\delta},b_{\delta})\times\prod_{i=1}^{N}IG(\sigma_{i}^{2}\,|\,a_{\sigma_{i}},b_{\sigma_{i}})$ in (7), and follow the other specifications as in the two-stage setting in Example 2.

As the priors have been chosen to be fully conjugate, one can derive the full posterior conditional distributions for each of the parameters. Specifically, the variance parameters will have posterior distributions of the form $IG(a^{*},b^{*})$, while the rest of the parameters will have posterior distributions of the form $N(Mm,M)$. However, as only $n$ of the $N$ groups are observed, the variance terms of the unsampled groups, $\sigma_{n+1}^{2},\dots,\sigma_{N}^{2}$, must either be fixed or given informative priors. If not, draws from the posterior predictive distribution corresponding to units in the nonsampled groups will have arbitrary variability and could spuriously dominate the finite population estimates.

Model 2. Spatial

Under (7), the intercept-only spatial model defines $\theta=[\phi,\delta^{2},\sigma^{2},\tau^{2}]^{\top}$ with corresponding prior distribution $p(\theta)=p(\phi)\times IG(\delta^{2}\,|\,a_{\delta},b_{\delta})\times IG(\sigma^{2}\,|\,a_{\sigma},b_{\sigma})\times IG(\tau^{2}\,|\,a_{\tau},b_{\tau})$ and $V_{\beta}=\delta^{2}$, $V_{s}=\Omega_{s}+\sigma^{2}I_{k}$. As there are no group terms, replace $X_{s}$ with $1_{k}$ and take $\nu=0$ with probability 1, e.g. $V_{\nu}^{-1}=0$.

Unlike model 1, regardless of the prior distribution placed on $p(\phi)$, a closed-form posterior distribution cannot be found for $\phi$. In practice, $\phi$ is often fixed using an estimate found from a variogram and then full posterior conditional distributions can be found using the same techniques described for the non-spatial case. However, MCMC can still be implemented by specifying a prior distribution for $\phi$ (Banerjee et al. 2014), which is often taken to be a uniform distribution.

To recover the spatial effects $\omega$ absorbed into the variance parameter of $y$, note that $\beta=\mu_{s}$ and $p(\omega\,|\,y,\theta,\mu_{s})\propto N(\omega\,|\,0,\Omega)\times N(y\,|\,1_{T}\mu_{s}+\omega,\sigma^{2}I_{T})\propto N(M_{\omega}m_{\omega},M_{\omega})$, where $m_{\omega}=\frac{1}{\sigma^{2}}(y-1_{T}\mu_{s})$ and $M_{\omega}=(\Omega^{-1}+\frac{1}{\sigma^{2}}I_{T})^{-1}$. Thus, drawing one $\omega\sim N(M_{\omega}m_{\omega},M_{\omega})$ for each posterior sample of $\{\theta,\mu_{s},y_{ns}\}$ will result in a set of posterior samples of $\omega$.

Model 3. Two-Stage + Spatial

The spatial model in (5) can be rewritten using (7) by letting $\theta=[\phi,\gamma^{2},\delta^{2},\tau^{2},\sigma_{1}^{2},\dots,\sigma_{N}^{2}]^{\top}$, with $p(\theta)=p(\phi)\times IG(\gamma^{2}\,|\,a_{\gamma},b_{\gamma})\times IG(\delta^{2}\,|\,a_{\delta},b_{\delta})\times IG(\tau^{2}\,|\,a_{\tau},b_{\tau})\times\prod_{i=1}^{N}IG(\sigma_{i}^{2}\,|\,a_{\sigma_{i}},b_{\sigma_{i}})$, $V_{\nu}=\gamma^{2}$, $V_{\beta}=\delta^{2}I_{N}$, $A=1_{N}$, $V_{s}=\Omega_{s}+V_{s}^{(\sigma)}$, and $\beta=\mu$. After posterior samples of $\{\theta,\mu,y_{ns}\}$ are drawn as described in (8), posterior samples of the spatial effects can be recovered by sampling one $\omega\sim N(M_{\omega 2}m_{\omega 2},M_{\omega 2})$ for each posterior sample of $\{\theta,\mu,y_{ns}\}$, where $m_{\omega 2}=V^{(\sigma)-1}(y-X\mu)$ and $M_{\omega 2}=(\Omega^{-1}+V^{(\sigma)-1})^{-1}$.

Model 4. Regional Spatial

To rewrite the region-specific spatial model given using (7), let $V_{\nu}=\gamma^{2}$, $V_{\beta}=\delta^{2}I_{N}$, $A=1_{N}$, $V_{s}=\Omega_{s*}+V_{s}^{(\sigma)}$, and $\beta=\mu$. Also take $\theta=[\phi_{1},\dots,\phi_{N},\gamma^{2},\delta^{2},\tau_{1}^{2},\dots,\tau_{N}^{2},\sigma_{1}^{2},\dots,\sigma_{N}^{2}]^{\top}$ with $p(\theta)=\prod_{i=1}^{N}p(\phi_{i})\times IG(\gamma^{2}\,|\,a_{\gamma},b_{\gamma})\times IG(\delta^{2}\,|\,a_{\delta},b_{\delta})\times\prod_{i=1}^{N}IG(\tau_{i}^{2}\,|\,a_{\tau_{i}},b_{\tau_{i}})\times\prod_{i=1}^{N}IG(\sigma_{i}^{2}\,|\,a_{\sigma_{i}},b_{\sigma_{i}})$. Similar to model 1, as not all locations are sampled, informative priors must be placed on the $\phi_{n+1},\dots,\phi_{N}$ spatial decay parameters. Additionally, to recover posterior samples of $\omega$, sample one $\omega\sim N(M_{\omega 3}m_{\omega 3},M_{\omega 3})$ for each posterior sample of $\{\theta,\mu,y_{ns}\}$ drawn using (8), where $m_{\omega 3}=V^{(\sigma)-1}(y-X\mu)$ and $M_{\omega 3}=(\Omega_{*}^{-1}+V^{(\sigma)-1})^{-1}$.

To achieve computation efficiency, redefine $y=[y_{1}^{\top},y_{1}^{*\top},\dots,y_{n}^{\top},y_{n}^{*\top},y_{n+1}^{*\top},\dots y_{N}^{*\top}]^{\top}$ so that the outcome is organized by region and then $\Omega_{*}$ becomes a $T\times T$ block diagonal matrix composed of $N$ blocks. This allows us to instead invert $N$ covariance matrices of size $M_{1}\times M_{1},\dots,M_{N}\times M_{N}$, rather than one $T\times T$ matrix, in the estimation of $\omega$.

II.   Exact Monte Carlo Estimation

If we are able to provide reasonable fixed values of the parameters, (2) can be simplified into a conjugate Bayesian linear model resembling:

[TABLE]

For a model such as this, the components $a$, $b$, $\tilde{V}_{\nu}$, $\tilde{V}_{\beta}$, and $\tilde{V}_{s}$ are fixed, reducing the model to three unknown parameters, $\delta^{2}$, $\nu$, and $\beta$. Thus, we can avoid MCMC and sample from the joint posterior $p(\delta^{2},\nu,\beta\,|\,y_{s})$ using the following steps. First sample $\delta^{2}$ from $p(\delta^{2}\,|\,y-S)=IG(a^{*},b^{*})$ and then for each $\delta^{2}$ drawn, draw a corresponding $\nu$ from $N(M_{\nu}m_{\nu},\delta^{2}M_{\nu})$. Next, for each pair of $\{\delta^{2},\nu\}$, draw $\beta$ from $N(M_{\beta}m_{\beta},\delta^{2}M_{\beta})$ (see the Supporting Information for details). As an example, we recast each model presented in Section I in the form of (9) and derive the posterior conditional distributions for model 1 and model 2, details of which are provided in the Supporting Information.

Model 1. Two-Stage

To create a conjugate Bayesian model such as (9) from the non-spatial model, define $A=1_{N}$, $\tilde{V}_{\nu}=\tilde{\gamma}^{2}=\frac{\gamma^{2}}{\delta^{2}}$, $\tilde{V}_{\beta}=I_{N}$, and $\tilde{V}_{s}=\tilde{V}_{s}^{(\sigma)}=[\oplus_{i=1}^{n}\frac{\sigma_{i}^{2}}{\delta^{2}}I_{m_{i}}]$. Noting that $p(\nu\,|\,y_{s})\propto N(\nu\,|\,0,\delta^{2}\tilde{\gamma}^{2})\times N(y_{s}\,|\,X_{s1}1_{n}\nu,\delta^{2}[X_{s1}X_{s1}^{\top}+\tilde{V}_{s}^{(\sigma)}])$ a little algebra reveals

[TABLE]

where $c=\frac{\sum_{i=1}^{n}\lambda_{i}\bar{y}_{i}}{\frac{1}{\tilde{\gamma}^{2}}+\sum_{i=1}^{n}\lambda_{i}}$ and d=\Big{[}\frac{1}{\tilde{\gamma}^{2}}+\sum_{i=1}^{n}\lambda_{i}\Big{]}^{-1}. The mean of the posterior distribution, $c$, is the weighted average of the sampled group means, where each mean is weighted by a function of each group’s element-wise variance. Integrating out $\nu$ and $\mu$ from $p(\delta^{2},\nu,\mu\,|\,y_{s})$ yields $p(\delta^{2}\,|\,y_{s})$, which is:

[TABLE]

Taking the limits of $c$ and $d$ as $\tilde{\gamma}^{2}\to\infty$ (e.g. $\gamma^{2}\to\infty$) we recover the findings of Scott and Smith (1969), who assigned $p(\nu)\propto 1$:

[TABLE]

As $p(\mu_{s}\,|\,y_{s},\delta^{2},\nu)\propto N(\mu_{s}\,|\,\nu 1_{n},\delta^{2}I_{n})\times N(y_{s}\,|\,X_{s1}\mu_{s},\delta^{2}\tilde{V}_{s}^{(\sigma)})$ we have that:

[TABLE]

where $c_{*}=\begin{bmatrix}(1-\lambda_{1})\nu+\lambda_{1}\bar{y}_{1}\\ \vdots\\ (1-\lambda_{n})\nu+\lambda_{n}\bar{y}_{n}\end{bmatrix}$ and $d_{*}=\begin{bmatrix}\oplus_{i=1}^{n}(1-\lambda_{i})\end{bmatrix}$. The posterior mean is appealing for interpretation, as its $i$-th element is the weighted average of the $i$-th group’s sample mean and the superpopulation mean estimate. Finally, note that $\mu_{ns}\,|\,y_{s},\nu,\delta^{2}\sim N(\mu_{ns}\,|\,\nu,\delta^{2}I_{N-n})$, as $y_{s}$ provides no information pertaining to the nonsampled groups.

Model 2. Spatial

The spatial model can be recast as (9) by defining $\tilde{V}_{\nu}^{-1}=0$ and $\tilde{V}_{s}=\tilde{\Omega}_{s}=\frac{1}{\delta^{2}}\Omega_{s}+I_{k}$, where $\tilde{V}_{\beta}$ is fixed to $1$. Defining $V_{\tilde{\Omega}_{s}}=(1+1_{k}^{\top}\tilde{\Omega}_{s}^{-1}1_{k})^{-1}$, the posterior conditionals are:

[TABLE]

[TABLE]

Model 3. Two-Stage + Spatial

The form of (9) is achieved by defining $\tilde{V}_{\nu}=\tilde{\gamma}^{2}$, $\tilde{V}_{\beta}=I_{N}$, $A=1_{N}$, and $\tilde{V}_{s}=\frac{1}{\delta^{2}}\Omega_{s}+\tilde{V}_{s}^{(\sigma)}$.

Model 4. Regional Spatial

The form of (9) is achieved by defining $\tilde{V}_{\nu}=\tilde{\gamma}^{2}$, $\tilde{V}_{\beta}=I_{N}$, $A=1_{N}$, and $\tilde{V}_{s}=\frac{1}{\delta^{2}}\Omega_{s*}+\tilde{V}_{s}^{(\sigma)}$.

III.   Model Comparison and Assessment

Model fit was evaluated in two ways. In general, consider a sample of size $k$ drawn from a population of size $T$ with outcome $y=[y_{s}^{\top}\colon y_{ns}^{\top}]^{\top}$. Without loss of generality, say $y_{h}\in y_{s}$ if $h=1,\dots,k$ and $y_{h}\in y_{ns}$ if $h=k+1,\dots,T$. First we evaluate the predictive accuracy of the models using the Watanabe-Akaike Information Criteria (WAIC), which is expressed as $WAIC=-2\widehat{elpd}=-2(\widehat{lpd}+\hat{p}_{WAIC})$ in Vehtari et al. (2017), where $\widehat{elpd}$ is the estimated expected log pointwise predictive density and is multiplied by $-2$ to be on the deviance scale. To calculate this, at each iteration, $l=1,\dots,L$, $p(y_{h}\,|\,\Theta^{(l)})$ is computed; the likelihood of each observed value conditional on that iteration’s parameters. The estimated log pointwise predictive density is the sum of the log average likelihood for each observation, \widehat{lpd}=\sum_{h=1}^{k}\log\Big{[}\frac{1}{L}\sum_{l=1}^{L}p(y_{h}\,|\,\Theta^{(l)})\Big{]}. The sample variance of the log-likelihood for each observation is s_{lp(y_{h})}^{2}=\frac{1}{L-1}\sum_{l=1}^{L}\Big{[}\log(p(y_{h}\,|\,\Theta^{(l)}))-\frac{1}{L}\sum_{l=1}^{L}\log(p(y_{h}\,|\,\Theta^{(l)}))\Big{]}^{2} and the estimated effective number of parameters is the sum of these variances: $\hat{p}_{WAIC}=\sum_{h=1}^{k}s_{lp(y_{h})}^{2}$. To calculate the standard error of the WAIC, rewrite -2(\widehat{lpd}+\hat{p}_{WAIC})=-2\sum_{h=1}^{k}\widehat{elpd}_{h}=\sum_{h=1}^{k}\Big{\{}\log\Big{[}\frac{1}{L}\sum_{l=1}^{L}p(y_{h}\,|\,\Theta^{(l)})\Big{]}+s_{lp(y_{h})}^{2}\Big{\}}. Under the assumption that each $\widehat{elpd}_{h}$ is independent, the sample variance of each individual $\widehat{elpd}_{h}$ is s_{elpd,ind}^{2}=\frac{1}{N-1}\sum_{h=1}^{k}\Big{[}\widehat{elpd}_{h}-\frac{1}{k}\sum_{h=1}^{k}\widehat{elpd}_{h}\Big{]}. Then $SE(WAIC)=\sqrt{Var(-2\sum_{h=1}^{k}\widehat{elpd}_{h})}=2\sqrt{nVar(\widehat{elpd}_{h})}=2s_{elpd,ind}\sqrt{n}$.

Second, for simulated data the true values $y=[y_{s}:y_{ns}]^{\top}$ are known and so we compare these with replicated datasets, $y_{rep}^{(l)}=[y_{rep,1}^{(l)}\dots y_{rep,k}^{(l)}]^{\top}$, generated from the pointwise posterior predictive distribution at each iteration $l$. These are used to formulate the goodness of fit measurement $D=G+P$ described in Gelfand and Ghosh (1998), composed of an error sum of squares term and a penalty term for large predictive variances. For $L$ iterations, $G=\sum_{h=1}^{k}(y_{h}-\text{E}[y_{rep,h}\,|\,y_{s}])^{2}$ and $P=\sum_{h=1}^{k}var(y_{rep,h}\,|\,y_{s})$. We approximate $\text{E}[y_{rep,h}\,|\,y_{s}]\approx\frac{1}{L}\sum_{l=1}^{L}y_{rep,h}^{(l)}$ and $var(y_{rep,h}\,|\,y_{s})\approx\frac{1}{L-1}\sum_{l=1}^{L}(y_{rep,h}^{(l)}-\frac{1}{L}\sum_{l=1}^{L}y_{rep,h}^{(l)})^{2}$. For non-simulated datasets, where $y_{ns}$ is unknown, $D$ can still be calculated by restricting the replicate datasets to the observed units, $y_{s}$, e.g. $y_{rep}^{(l)}=[y_{rep,1}^{(l)}\dots y_{rep,k}^{(l)}]^{\top}$, at the $l$-th iteration.

V.   Simulation

I.   Data Generation

To simulate spatial correlation and allow for two-stage random sampling, a unit square was divided into 100 equally sized square regions and 2,500 locations were randomly drawn from the unit square. Data was simulated from the intercept-only spatial model described in Model 2 with $\mu=2$. A distance matrix for all points was constructed and used to create a covariance matrix that reflects an exponential covariance function described in Section III, where $\phi$ was assigned a value of 10, reflecting an effective spatial range of 3/10. The spatial variance, $\tau^{2}$, was fixed at 9, while the non-spatial variance, $\sigma^{2}$, was set to 4. After a dataset was generated, a cluster random sampling scheme was implemented. 25 regions were randomly selected and then in each cluster, a random number of individuals were selected (the minimum and maximum percent of those selected from a region was set to be 20% and 90%, respectively). 20 datasets containing information of both the sampled and nonsampled units were generated in this way. To examine Models 1 and 4 for larger dataset, this process was then repeated with the same parameters to generate 20 datasets with 8,100 locations from 324 regions, where 81 regions were randomly sampled. All data generation and analyses were performed using R version 3.5.1 R Core Team 2018.

II.   Exact Monte Carlo Simulation

To perform the two-stage procedure using the conditional distributions and methods described in Section II, sample means, $\hat{\mu}_{i}$, and sample variances, $\hat{\sigma}_{i}^{2}$, were calculated from each observed cluster, $i=1,\dots,n$. The variance matrix of the sampled units was fixed to be $V_{s}^{(\sigma)}=\begin{bmatrix}\oplus_{i=1}^{n}\frac{\hat{\sigma}_{i}^{2}}{\text{Var}(\hat{\mu})}I_{m_{i}}\end{bmatrix}$, where $\text{Var}(\hat{\mu})$ represents the sample variance of the observed sample means. Similarly, the variance matrix of the nonsampled units was fixed at $V_{ns}=\begin{bmatrix}\oplus_{i=1}^{N}\frac{\tilde{\sigma}_{i}^{2}}{\text{Var}(\hat{\mu})}I_{M_{i}-m_{i}}\end{bmatrix}$, where $\tilde{\sigma}_{i}^{2}=\hat{\sigma}_{i}^{2}$ if $i\in\{1,\dots,n\}$ and $\tilde{\sigma}_{i}^{2}=\frac{1}{n}\sum_{i=1}^{n}\hat{\sigma}_{i}^{2}$ if $i\in\{n+1,\dots,N\}$. The value of $\gamma^{2}$ was fixed to be half of the value of $\delta^{2}$, reflecting the belief that there was less variability in the population mean than between group means. The prior distribution for $\delta^{2}$ was assigned to be $IG(3,5)$. Sampling from the posterior was performed using the conditional distributions and methods described in Section II. At each iteration $g$, the population mean estimate for that iteration was then calculated as \bar{y}^{(g)}=\frac{1}{T}\Big{(}\sum_{i=1}^{n}\sum_{j=1}^{m_{i}}y_{ij}+\sum_{i=1}^{N}\sum_{j=m_{i}+1}^{M_{i}}y_{ij}^{(g)}\Big{)}. Details of this iterative procedure can be found in the supplementary materials.

To perform the spatial random effect procedure, $\phi$ was set to its true value of 10 and the ratio of $\delta^{2}/\tau^{2}$ to its true value of $4/9$. The posterior conditionals of $\delta^{2}\,|\,y_{s}$ and $\mu_{s}\,|\,y_{s},\delta^{2}$, (13) and (14) respectively, were sampled as outlined in Section III. This sampling and the prediction of $y_{ns}$ was performed using commands from the spBayes R package (Finley et al., 2015, 2007). The population mean estimate was calculated using the technique described in the non-spatial sampling case above.

Figure 1 plots population average-centered mean estimates and 95% credible intervals from both methods applied to the twenty simulated datasets. While the spatial cases consistently have a smaller credible interval, their point estimates are similar to the two-stage case. However, as the ratio of spatial and non-spatial variance is fixed for this method, this may result in a reduction in the overall variance of the population mean. Posterior mean estimates and their associated 95% credible intervals of the superpopulation parameters and finite population mean, $\bar{y}$, from the first generated dataset are given in Table 1, along with the WAIC, its standard error, and $D$ values. While both models have similar estimates for $\nu$, the two-stage model overestimates the non-spatial variance. This is expected, as we know that there is additional variance due to spatial correlation that is not being accounted for otherwise in the model. Similarly, both measures of goodness of fit prefer the spatial model.

III.   Markov Chain Monte Carlo Simulation

To explore these findings further, we implemented the four models described in Section IV using the JAGS software in R on the same generated datasets. Models were run for 650 iterations with 50 burn-in, as examination of individual trace plots suggested sufficient mixing and convergence of the non-spatial parameters. At each iteration $g$, estimates of the nonsampled units were drawn and estimates for the population mean, \bar{y}^{(g)}=\frac{1}{T}\Big{(}\sum_{i=1}^{n}\sum_{j=1}^{m_{i}}y_{ij}+\sum_{i=1}^{N}\sum_{j=m_{i}+1}^{M_{i}}y_{ij}^{(g)}\Big{)} were calculated. All variance parameters ($\sigma^{2}$, $\tau^{2}$, and $\delta^{2}$, as well as site-specific variances such as $\sigma_{i}^{2}$ and $\tau_{i}^{2}$) were given an inverse-gamma prior with shape 2 and scale 10, reflecting a weakly-informative prior distribution with mean $10$. There is a substantial literature in theoretical spatial statistics regarding the identifiability, or lack there of, of the spatial process parameters, hence, non-informative or completely flat priors are excluded from consideration. The prior families and specifics we use are fairly customary in spatial modeling. They exploit some information about the spatial domain and the extent of spatial association that can be detected from finite samples using variograms. For example, in practice, given a real data set, we would pass the data through an exploratory analysis tool, such as a variogram, glean some information about the spatial variance component and the measurement error component, and use the weakly informative centered inverse-gamma priors to reflect these values. In addition, $\nu$ was given a flat prior to not inform the estimation of the mean and all $\phi$ parameters were given Uniform(5,15) priors to allow the spatial range to vary from 0.2 (3/15) to 0.6 (3/5). While our priors were chosen to be weakly-informative to be conservative in estimation, more informative priors could easily be added in a data analysis if additional information regarding the parameters was known. MCMC sampling was performed using the computer program JAGS (Plummer, 2017) in R.

When assessing model fit in the first realization of the data with WAIC, the spatial model performed slightly worse than the rest of the models with a value of 1,912.70 (SE = 26.36). This may be due to the additional variation which comes from varying the spatial range parameter. This was followed closely by the two-stage model with 1,870.26 (35.00) , which was outperformed by both the regional spatial model with 1,202.08 (38.99) and the two-stage + spatial model with 455.67 (17.03). It is interesting that while the data was generated by the spatial model and sampled by a two-stage framework, neither of these models perform better than the two models which take both the spatial correlation and study design into account. Additional comparisons of estimates from this first realization can be found in the supplementary materials.

Figure 2 shows the models’ posterior mean estimates of the finite population mean, which are centered at the true population mean and presented with 95% credible intervals for the 20 simulated datasets. While point estimates remain similar across models, the best fitting model, Model 3, has the widest credible intervals for the population mean. Accounting for only regional effects results in tight credible intervals in Models 1 and 4, which are narrow compared to Model 2, which fails to take into account region specific variability. Similar results were found when applying Models 1 and 4 to the larger simulated datasets and are provided in the supplementary materials.

VI.   Data Analysis: Nitrate in Central California Groundwater

In this section, we provide an analysis of groundwater nitrate content of the Tulare Lake Basin (TLB) in Central California from the California Ambient Spatio-Temporal Information on Nitrate in Groundwater (CASTING) Database, which is described in Harter et al. (2017) and Boyle et al. (2012) and is available as the University of California, Davis, nitrate data in the data repository of the Groundwater Ambient Monitoring and Assessment Program (2019). Interest lies in identifying regions in which ground water nitrate levels exceed 45 mg/L, which is the maximum contaminant level established by the EPA Boyle et al. (2012). At high levels, infants and pregnant women are more susceptible to nitrate poisoning, which makes it more difficult for oxygen to be distributed to body and can be fatal to infants less than six months old. Besides human sources such as sewage disposal, many sources of nitrate are agricultural, such as fertilizer for crops and animal waste (Harter and Lund, 2012). Because of this, regions with high agricultural activity, such as the TLB, have experienced rising levels of nitrate over the past few decades. As groundwater, and therefore nitrate levels in groundwater, can be assumed to be present at all areas of the Central Valley, we can assume that water samples taken from wells come from a spatial field. Therefore, given a sample of readings from various wells, our primary goal is to estimate of the finite population average of all known wells, which represents an overall measure of water-health. Additionally, plots of posterior predictive distribution may be useful in identifying high-risk regions which exceed the maximum contaminant level.

The CASTING Database is an extensive collection of nitrate readings from the TLB and Salinas Valley collected by multiple agencies, over 70% of which were collected between 2000 and 2011. Of these, most wells had repeated measurements taken over the time. As the Salinas Valley and the TLB are geographically separate regions of California, only the TLB was included. In order to avoid associations over time, the data was restricted to a single year. The year 2009 was selected as variogram plots suggested nitrate levels followed a roughly exponential distribution. While directional variograms suggested that the measurements may be anisotropic, we continued with the methods presented above, recognizing that a model accounting for directional spatial dependence may provide a better fit to the data.

While the data roughly covers the TLB, we construct a likely sampling scenario in which the TLB is separated into distinct geographic regions and due constraints (perhaps time or financial), a random subset of these regions are sampled. In our scenario, zip codes were used as it is common to collect such geographic information in large scale health surveys, but many alternatives such as cities or grid-based approach could have also been used. A map of California zip codes tabulation areas obtained from the tigris R package (Walker, 2018) was overlaid on the approximate geographic locations of each of the sampled wells, effectively assigning each well to one specific region, defined by a zip code. 1) Only the most recent observation was taken from each well so that each well was only represented once. 2) If unique wells had the same geographic coordinates, one was chosen at random to be removed. 3) Sparse zip codes with less than 10 wells were excluded to ensure that each selected zip code would have a large sample size and to avoid overfitting when modeling. These restrictions resulted in a dataset with 6,117 unique wells among 63 zip codes. Nitrate level had a mean of 37.9 mg/L, standard deviation of 52.3 mg/L, and ranged from 0.0 to 903.1 mg/L.

In order to recreate a cluster sampling scenario, 21 of the zip codes were randomly chosen and 50-90% of the population in that zip code was randomly sampled. This resulted in an observed sample size of 489 with a mean nitrate level of 34.2 mg/L and a standard deviation of 40.0 mg/L. The nitrate level ranged from 0.0 to 269.6. A plot of these sampled and non-sampled zip codes is available in the supplementary materials. Using this sampled data, all four models in Section III were implemented and the results are shown in Table  2. In all models, $\nu$ was given a flat prior. For model 1, the regional variance parameters were given an inverse-gamma prior with shape 2 and scale 1600 to reflect the sample variance. For the spatial models, a variogram was fit and spatial variances were given an inverse-gamma prior with shape 2 and scale 1800. Similarly, non-spatial variance terms were assigned an inverse-gamma prior with shape 2 and scale 1000. The variance of regional means was assigned an inverse-gamma prior with shape 2 and scale 10 to allow for small, localized deviations. All $\phi$ parameters were given Uniform(0.01,5) priors to reflect a spatial range varying between 0.6 km (3/5) and 300 km (3/0.01). MCMC sampling was performed using JAGS (Plummer, 2017) in R (R Core Team, 2018).

With respect to the estimate of the true mean nitrate level, only the intercept-only spatial model contained the true mean value within its 95% credible intervals. However, as evidenced by the larger mean, standard deviation, and range in the complete dataset, it appears that the sampled units did not capture some of the larger outliers, so it is unsurprising that the estimates of the population mean are lower than the truth. Comparing WAIC, we see results similar to those found in Section III. The spatial models which accounted for regional means had lower WAIC values than the two-stage model, which is evidence that this data is spatially correlated. However, the intercept-only spatial model did not fit the data as well as the two-stage model, which may be due to ignoring the study design. Additionally, the two-stage + spatial model again fits the model the best.

Figure 3 shows the interpolated population surface from the complete sample and the interpolated surface from posterior predictive samples. While there are common regions at high risk (nitrate level greater than 45 mg/L) in all the posterior predictive maps, the spatial and two-stage + spatial maps predict larger regions. Also seen in Table 2, it is clear that Model 3 estimates a population mean that is larger than Models 1 and 4, but smaller than Model 3. Spatial residual plots for Models 2, 3, and 4 are provided in the supplementary materials.

VII.   Discussion

This paper examines the implications of performing two-stage random sampling on point-referenced data which exists in a spatial field. While Scott and Smith (1969) and Malec and Sedransk (1985) provided a Bayesian model-based framework to account for such a study design, we have demonstrated that an analysis ignoring the underlying spatial correlation between locations or sampling design may lead to spurious inference and poorer model fit.

This work is a first step in developing an overarching framework for Bayesian finite population sampling from spatial process based populations. In our two-stage case, additional work may be done to further improve this model. For instance, CAR priors could be placed on regional parameters such as the $\mu_{i}$’s, the regional means, to induce additional spatial correlation in the model. Also, many of the models presented account for regional differences in variance but if other sources of heteroscedasticity are suspected, new approaches (such as Zangeneh and Little 2015) may be needed to account for this. While an exponential covariance function was employed in the analyses in this paper, other spatial covariance functions could be used to create similar simulations and data analyses, as well as account for anisotropy.

Future work is needed to establish a more general framework that can account for more sophisticated sampling designs in a spatial context. The sampling designs presented in this paper are said to be ignorable (Rubin 1976; Sugden and Smith 1984), which allows us to perform inference on the superpopulation parameters while ignoring the inclusion probability distribution. However, designs in which the data cannot be assumed to be missing at random or where parameters define both the outcome and inclusion distributions are referred to as nonignorable and must account for the inclusion probability distribution. One example of this in the spatial context is preferential sampling (Diggle et al. 2010, Gelfand et al. 2012), in which the measurement values and sampling strategy are assumed to stem from the same spatial process. While Pati et al. (2011) have analyzed such data using Bayesian hierarchical models, an overall framework is needed to account for this and other non-ignorable design types.

Additionally, the implications of study design on finite population estimates when sampling from a spatially correlation population over time are unknown. In order to better understand this, these Bayesian models must first be extended to account for both study design and spatio-temporal associations.

Finally, while this paper provided a scaleable model which can account for study design and spatial correlation in massive survey data by assuming regional independence, further work should be done to incorporate recent strategies in modeling large spatial data (Heaton et al., 2018) when analyzing survey data with spatial correlations, such as nearest neighbor processes (Datta et al., 2016), covariance tapering (Furrer et al., 2006), and metakriging (Guhaniyogi and Banerjee, 2018). Finite population models would particularly benefit from such techniques, as computation increases as a function of the population total, $T$, rather than the sample size, $k$.

Supporting Information

The work of the first and second authors were supported, in part, by the Division of Information and Intelligent Systems of the National Science Foundation under Grant 1562303, the Division of Mathematical Sciences of the National Science Foundation under Grant 1916349, and the National Institute of Environmental Health Sciences of the National Institutes of Health under Grants 1R01ES027027 and R01ES030210-01.

Appendix A

This section first provides the derivation of the empirical Bayesian estimators presented in Section II and then the finite population estimates presented in Sections II and III.

The values $a^{*}$, $b^{*}$, $M_{\nu}$, $m_{\nu}$, $M_{\beta}$, and $m_{\beta}$ for the general case (9) in Section II are presented below.

[TABLE]

The derivation of these values arise the conjugacy of the Normal and Inverse-Gamma distributions. We first derive (11) and (12). Take $\epsilon_{s}\sim N(0,\delta^{2}\tilde{V}_{s}^{(\sigma)})$, $\epsilon_{ns}\sim N(0,\delta^{2}\tilde{V}_{ns}^{(\sigma)})$, and $\nu\sim N(0,\delta^{2}\tilde{\gamma}^{2})$, where $\tilde{\gamma}^{2}=\frac{1}{\delta^{2}}\gamma^{2}$, $\tilde{V}_{s}^{(\sigma)}=\frac{1}{\delta^{2}}V_{s}^{(\sigma)}=\begin{bmatrix}\oplus_{i=1}^{n}\frac{\sigma_{i}^{2}}{\delta^{2}}I_{m_{i}}\end{bmatrix}$, and $\tilde{V}_{ns}^{(\sigma)}=\frac{1}{\delta^{2}}V_{ns}^{(\sigma)}=\begin{bmatrix}\oplus_{i=1}^{n}\frac{\sigma_{i}^{2}}{\delta^{2}}I_{(M_{i}-m_{i})}\end{bmatrix}$. Since the elements of $y$ are independent conditional on $\mu$, $V_{s,ns}=0$ and $V_{ns,s}=0$. Define observed group means as $\bar{y}_{i}=\frac{1}{m_{i}}\sum_{j=1}^{n}y_{ij}$ and the ratio of variances as $\lambda_{i}=\delta^{2}/(\delta^{2}+\sigma_{i}^{2}/m_{i})$ if $i\in(1,\dots,n)$ and $\lambda_{i}=0$ if $i\in(n+1,\dots,N)$. Also define the vector of observed group variances to be $\tilde{\sigma}^{2}=[\sigma_{1}^{2},\dots,\sigma_{n}^{2}]^{\top}$. Recall $\mu_{s}=\nu 1_{n}+\epsilon_{\mu_{s}};\epsilon_{\mu_{s}}\sim N(0,\delta^{2}I_{n})$, then $y_{s}=X_{s1}\mu_{s}+\epsilon_{s}=X_{s1}1_{n}\nu+e_{s}^{*}$, where $e_{s}^{*}\sim N(0,\delta^{2}[X_{s1}X_{s1}^{\top}+\tilde{V}_{s}^{(\sigma)}])$. Then we have that $p(\nu\,|\,y_{s})\propto N(\nu\,|\,0,\delta^{2}\tilde{\gamma}^{2})\times N(y_{s}\,|\,X_{s1}1_{n}\nu,\delta^{2}[X_{s1}X_{s1}^{\top}+\tilde{V}_{s}^{(\sigma)}])\propto N(\nu\,|\,Bb,\delta^{2}B)$, where

[TABLE]

Therefore, B=\Big{[}\frac{1}{\tilde{\gamma}^{2}}+\sum_{i=1}^{n}\lambda_{i}\Big{]}^{-1} and $Bb=\frac{\sum_{i=1}^{n}\lambda_{i}\bar{y}_{i}}{\frac{1}{\tilde{\gamma}^{2}}+\sum_{i=1}^{n}\lambda_{i}}$. To solve $p(\delta^{2}\,|\,y_{s})$, split the posterior conditional distribution of the superpopulation parameters; $p(\delta^{2},\nu\,|\,y_{s})=IG(\delta^{2}\,|\,a_{\delta}^{*},b_{\delta}^{*})\times N(\nu\,|\,Bb,\delta^{2}B)$, where $a_{\delta}^{*}=a_{\delta}+\frac{k}{2}$ and $b_{\delta}^{*}=b_{\delta}+\frac{1}{2}[y_{s}^{\top}(X_{s1}X_{s1}^{\top}+\tilde{V}_{s}^{(\sigma)})^{-1}y_{s}+b^{\top}Bb]$. As $p(\mu_{ns}\,|\,y_{s},\nu,\delta^{2})=p(\mu_{ns}\,|\,\nu,\delta^{2})$, $\mu_{ns}\,|\,\nu,\delta^{2}\sim N(1_{N-n}\nu,\delta^{2}I_{N-n})$. To solve $p(\mu_{s}\,|\,y_{s},\nu,\delta^{2})$, note that $p(\mu_{s}\,|\,y_{s},\nu,\delta^{2})\propto N(\mu_{s}\,|\,1_{n}\nu,\delta^{2}I_{n})\times N(y_{s}\,|\,X_{s1}\mu_{s},\delta^{2}\tilde{V}_{s}^{(\sigma)})\propto N(\mu_{s}\,|\,B_{*}b_{*},\delta^{2}B_{*})$, where

[TABLE]

To derive (13) and (14) define $\tilde{V}_{\nu}^{-1}=0$ and $\tilde{V}_{s}=\tilde{\Omega}_{s}=\frac{1}{\delta^{2}}\Omega_{s}+I_{k}$, and $\tilde{V}_{\mu_{s}}=1$. Note that $p(\mu_{s}\,|\,y_{s},\delta^{2})\propto N(\mu_{s}\,|\,0,\delta^{2})\times N(y_{s}\,|\,1_{k}\mu_{s},\delta^{2}\tilde{\Omega}_{s})\propto N(\mu_{s}\,|\,B_{\mu_{s}}b_{\mu_{s}},\delta^{2}B_{\mu_{s}})$, where $B_{\mu_{s}}=(1+1_{k}^{\top}\tilde{\Omega}_{s}^{-1}1_{k})^{-1}$ and $b_{\mu_{s}}=1_{k}^{\top}\tilde{\Omega}_{s}^{-1}y_{s}$. Splitting the posterior conditional distribution of the superpopulation parameters, $p(\delta^{2},\mu_{s}\,|\,y_{s})=IG(\delta^{2}\,|\,a_{\delta}^{**,}b_{\delta}^{**})\times N(\mu_{s}\,|\,B_{\mu_{s}}b_{\mu_{s}},\delta^{2}B_{\mu_{s}})$, where $a_{\delta}^{**}=a+\frac{k}{2}$ and b_{\delta}^{**}=b+\frac{1}{2}y_{s}^{\top}\Big{\{}\tilde{\Omega}_{s}^{-1}-\tilde{\Omega}_{s}^{-1}1_{k}\Big{(}1+1_{k}^{\top}\tilde{\Omega}_{s}^{-1}1_{k}\Big{)}^{-1}1_{k}^{\top}\tilde{\Omega}_{s}^{-1}\Big{\}}y_{s}.

We now continue by deriving the general cases presented in Section II. As $\beta\,|\,\nu\sim\text{N}(A\nu,V_{\beta})$ and $\nu\sim\text{N}(0,\gamma^{2})$, we have that $\beta\sim\text{N}(0,\gamma^{2}AA^{\top}+V_{\beta})$. Note that $p(\beta\,|\,y_{s})\propto N(0,\gamma^{2}AA^{\top}+V_{\beta})\times N(X_{s}\beta,V_{s})\propto N(V_{\beta\,|\,y_{s}}X_{s}^{\top}V_{s}^{-1}y_{s},V_{\beta\,|\,y_{s}})$, where $V_{\beta\,|\,y_{s}}=[(\gamma^{2}AA^{\top}+V_{\beta})^{-1}+X_{s}^{\top}V_{s}^{-1}X_{s}]^{-1}$. Defining $B_{V}=V_{ns,s}+(X_{ns}-V_{ns,s}V_{s}^{-1}X_{s})V_{\beta\,|\,y_{s}}X_{s}^{\top}$ and $Q=X_{ns}-V_{ns,s}V_{s}^{-1}X_{s}$, we have that:

[TABLE]

To derive the estimate given in Example 1, note that $p(\mu\,|\,y_{s})\propto N(\mu\,|\,0,\xi^{2})\times N(y_{s}\,|\,1_{n}\mu,\sigma^{2}I_{n})\propto N(\mu\,|\,B_{srs}b_{srs},B_{srs})$, where $B_{srs}=(\frac{1}{\xi^{2}}+\frac{n}{\delta^{2}})^{-1}$, $b_{srs}=\frac{1}{\sigma^{2}}1_{n}^{\top}y_{s}$, and $B_{srs}b_{srs}=\frac{\frac{1}{\sigma^{2}}1_{n}^{\top}y_{s}}{\frac{1}{\xi^{2}}+\frac{n}{\sigma^{2}}}=\frac{\sum_{i=1}^{n}y_{i}}{\frac{\sigma^{2}}{\xi^{2}}+n}$. Fixing the variance components, the finite population estimate is

[TABLE]

To derive (4), it is helpful to first make a note regarding $X_{s1}$ vs $X_{s}$ in the calculation of $p(\nu\,|\,y_{s})$ and $p(\mu\,|\,\nu,y_{s})$. Specifically, $p(\nu\,|\,y_{s})$ does not change, since $X_{s}=[X_{s1}:0]$, $b=1_{n}^{\top}X_{s1}^{\top}(X_{s1}X_{s1}^{\top}+\tilde{V}_{s}^{(\sigma)})^{-1}y_{s}=1_{N}^{\top}X_{s}^{\top}(X_{s}X_{s}^{\top}+\tilde{V}_{s}^{(\sigma)})^{-1}y_{s}$. Similarly, $B^{-1}=\frac{1}{\tilde{\gamma}^{2}}+1_{n}^{\top}X_{s1}^{\top}(X_{s1}X_{s1}^{\top}+\tilde{V}_{s}^{(\sigma)})^{-1}X_{s1}1_{n}=\frac{1}{\tilde{\gamma}^{2}}+1_{N}^{\top}X_{s}^{\top}(X_{s}X_{s}^{\top}+\tilde{V}_{s}^{(\sigma)})^{-1}X_{s}1_{N}$. However, while computing $p(\mu_{s}\,|\,\nu,y_{s})$ using $X_{s1}$ is computationally convenient for interpretation, employing $X_{s}$ provides us the posterior distribution $p(\mu\,|\,\nu,y_{s})$. We have that $p(\mu\,|\,y_{s},\nu)\propto N(\mu\,|\,\nu 1_{N},\delta^{2}I_{N})\times N(y_{s}\,|\,X_{s}\mu,\delta^{2}\tilde{V}_{s}^{(\sigma)})\propto N(\mu\,|\,B_{**}b_{**},\delta^{2}B_{**})$. Some algebra simplifies the expressions for $b_{**}$ and $B_{**}b_{**}$ and matches the conclusions found by deriving $p(\mu_{s}\,|\,\nu,y_{s})$ and $p(\mu_{ns}\,|\,\nu,y_{s})$ separately:

[TABLE]

Using these derivations, define $\lambda=[\lambda_{1},\dots,\lambda_{N}]^{\top}$ and $\bar{y}=[\bar{y}_{1},\dots,\bar{y}_{n},0_{(N-n)}]^{\top}$.Then fixing the variance components, we have:

[TABLE]

Now consider the stratified case for estimating the population mean. Taking non-informative priors for the group means, $\mu$, is equivalent to letting $\delta^{2}\to\infty$ and $\gamma^{2}\to\infty$. Therefore $\lambda_{i}=\frac{\delta^{2}}{\delta^{2}+\sigma_{i}^{2}/m_{i}}\to 1$, for all $i=1,\dots,N$. Note $\alpha_{i}=\sum_{j=m_{i}+1}^{M_{i}}\frac{1}{T}=\frac{M_{i}-m_{i}}{T}$, $i=1,\dots,n$. We have that:

[TABLE]

To derive (6), note $p(\nu\,|\,y_{s},\tau^{2},\Omega_{s},V_{s}^{(\sigma)})\propto N(y_{s}\,|\,X_{s}1_{N}\nu,\delta^{2}X_{s}X_{s}^{\top}+\Omega_{s}+V_{s}^{(\sigma)})\times N(\nu\,|\,0,\gamma^{2})\propto N(\nu\,|\,B_{sp2}b_{sp2},B_{sp2})$, where $B_{sp2}=(\frac{1}{\gamma^{2}}+1_{N}^{\top}X_{s}^{\top}(\delta^{2}X_{s}X_{s}^{\top}+\Omega_{s}+V_{s}^{(\sigma)})^{-1}X_{s}1_{N})^{-1}$ and $b_{sp2}=1_{N}^{\top}X_{s}^{\top}(\delta^{2}X_{s}X_{s}^{\top}+\Omega_{s}+V_{s}^{(\sigma)})^{-1}y_{s}$. Consider the non-spatial case and define $\lambda^{\top}=[\lambda_{1},\dots,\lambda_{N}]=1_{N}^{\top}X_{s}^{\top}(X_{s}X_{s}^{\top}+\tilde{V}_{s}^{(\sigma)})^{-1}X_{s}$, then B=\Big{[}\frac{1}{\tilde{\gamma}^{2}}+1_{N}^{\top}X_{s}^{\top}(X_{s}X_{s}^{\top}+\tilde{V}_{s}^{(\sigma)})^{-1}X_{s}1_{N}\Big{]}^{-1}=\Big{[}\frac{1}{\tilde{\gamma}^{2}}+\lambda^{\top}1_{N}\Big{]}^{-1}=\Big{[}\frac{1}{\tilde{\gamma}^{2}}+\sum_{i=1}^{n}\lambda_{i}\Big{]}^{-1}, which agrees with our previous findings.

Similarly, define $\lambda^{*\top}=[\lambda_{1}^{*},\dots,\lambda_{N}^{*}]=1_{N}^{\top}X_{s}^{\top}(\delta^{2}X_{s}X_{s}^{\top}+\Omega_{s}+V_{s}^{(\sigma)})^{-1}X_{s}$.

Then B_{sp2}=\Big{[}\frac{1}{\tilde{\gamma}^{2}}+\lambda^{*\top}1_{N}\Big{]}^{-1}=\Big{[}\frac{1}{\tilde{\gamma}^{2}}+\sum_{i=1}^{n}\lambda_{i}^{*}\Big{]}^{-1}.

Additionally, $p(\mu\,|\,y_{s},\nu)\propto N(\mu\,|\,\nu 1_{N},\delta^{2}I_{N})\times N(y_{s}\,|\,X_{s}\mu,\Omega_{s}+V_{s}^{(\sigma)})\propto N(\mu\,|\,B_{sp2*}b_{sp2*},B_{sp2*})$. Here $B_{sp2*}=(\frac{1}{\delta^{2}}I_{N}+X_{s}^{\top}(\Omega_{s}+V_{s}^{(\sigma)})^{-1}X_{s})^{-1}$ and $b_{sp2*}=\frac{1}{\delta^{2}}1_{N}\nu+X_{s}^{\top}(\Omega_{s}+V_{s}^{(\sigma)})^{-1}y_{s}$.

Fixing the variance parameters, we have that:

[TABLE]

Appendix B

To perform the two-stage procedure implemented in Section II, sample means, $\hat{\mu}_{i}$, and sample variances, $\hat{\sigma}_{i}^{2}$, were calculated from each observed region, $i=1,\dots,25$. The variance matrix of the sampled units was fixed to be $\tilde{V}_{s}^{(\sigma)}=\begin{bmatrix}\oplus_{i=1}^{n}\frac{\hat{\sigma}_{i}^{2}}{\text{Var}(\hat{\mu})}I_{m_{i}}\end{bmatrix}$where $\text{Var}(\hat{\mu})$ represents the sample variance of the observed sample means. Similar to fixing $V_{s}^{(\sigma)}$, the variance matrix of the nonsampled units was fixed at $V_{ns}=\begin{bmatrix}\oplus_{i=1}^{N}\frac{\tilde{\sigma}_{i}^{2}}{\text{Var}(\hat{\mu})}I_{M_{i}-m_{i}}\end{bmatrix}$, where $\tilde{\sigma}_{i}^{2}=\hat{\sigma}_{i}^{2}$ if $i\in\{1,\dots,n\}$ and $\tilde{\sigma}_{i}^{2}=\frac{1}{n}\sum_{i=1}^{n}\hat{\sigma}_{i}^{2}$ if $i\in\{n+1,\dots,N\}$. The value of $\gamma^{2}$ was fixed to be half of the value of $\delta^{2}$, reflecting the belief that there was less variability in the population mean than between group means. The prior distribution for $\delta^{2}$ was assigned to be $IG(3,5)$. Sampling from the posterior was performed using the conditional distributions and methods described in Section II. As we have fixed the ratios of all variance components, we have also fixed $\tilde{\lambda}_{i}=\text{Var}(\hat{\mu})/(\text{Var}(\hat{\mu})+\hat{\sigma}_{i}^{2}/m_{i})$ if $i\in(1,\dots,n)$ and $\tilde{\lambda}_{i}=0$ if $i\in(n+1,\dots,N).$. Define $c=\frac{\sum_{i=1}^{n}\tilde{\lambda}_{i}\bar{y}_{i}}{\frac{1}{2}+\sum_{i=1}^{n}\tilde{\lambda}_{i}}$, d=\Big{[}\frac{1}{2}+\sum_{i=1}^{n}\tilde{\lambda}_{i}\Big{]}^{-1}, $c^{*(g)}=\begin{bmatrix}(1-\tilde{\lambda}_{1})\nu^{(g)}+\tilde{\lambda}_{1}\bar{y}_{1}\\ \vdots\\ (1-\tilde{\lambda}_{n})\nu^{(g)}+\tilde{\lambda}_{n}\bar{y}_{n}\end{bmatrix}$, and $d^{*}=\begin{bmatrix}\oplus_{i=1}^{n}(1-\tilde{\lambda}_{i})\end{bmatrix}$. The following procedure was implemented to produce posterior estimates of the population mean, $\bar{y}^{(g)}$, for $G$ iterations.

[TABLE]

Figure 4 recreates the centered mean plots presented in Figure 2 for the larger data case, in which the number of regions is 324. As in the $N=100$ case, the point estimates and 95% credible intervals are similar for the two models. Additionally, the regional spatial model still outperforms the two-stage model with a WAIC of 1,052 (SE = 38.89) compared to 1,869 (SE = 34.77). MCMC sampling was performed using JAGS (Plummer, 2017) in R (R Core Team, 2018).

Table 3 compares the finite population mean estimate, $\nu$ estimate, and model fit for the first simulated dataset. Recall that $\nu=2$ and the finite population mean was 2.60. Notice that while the true values are included in the credible intervals for all models, the credible intervals for $\nu$ are wider than those of the FP mean.

Appendix C

Figure 5 presents the 2010 zip code tabulation areas in the California Central Valley obtained from the tigris R package (Walker, 2018). Sixty-three zip codes were included in the analysis of a subset of observations from the UC Davis Nitrate Data in the data repository of the Groundwater Ambient Monitoring and Assessment Program (2019), presented in Section VI. These zip codes are denoted as either sampled or non-sampled, while all other zip codes are denoted as excluded.

The variograms provided in Figure 6 were fitted to the entire population of wells (left) and sampled wells (right). For the population variogram, the estimated values of the nugget, partial sill, and range were 2527.4, 319.9, and 2.29, respectively. For the sample variogram, the estimated values of the nugget, partial sill, and range were 1808.6, 1014.4, and 335.9, respectively.

Figure 7 provides spatial residual plots arising from the three spatial models. The spatial model which does not account for regional effects sees the most dispersed spatial effects, while the two-stage + spatial and regional spatial models show more localized spatial variability. MCMC sampling was performed using JAGS (Plummer, 2017) in R (R Core Team, 2018).