Bayesian spatial models for estimating means of sampled and non-sampled small areas
Section 2. Some spatial alternatives to the independent FH model

• Previous
• Next

2.1   Incorporating spatial random effects

Let $Y={\mathrm{<mo\; stretchy="false">(</mo>}{Y}_1\mathrm{,}\dots \mathrm{,}{Y}_m\mathrm{<mo\; stretchy="false">)</mo>}}^T$ be the $m$ -component vector with the direct estimates of $m$ small areas, and $D=\text{diag}{\mathrm{<mo>\{</mo>}{D}_i\mathrm{<mo>\}</mo>}}_{i=1}^m$ be the $m\times m$ diagonal matrix with the sampling variances of the direct estimates. We denote by $\theta ={\mathrm{<mo\; stretchy="false">(</mo>}{\theta }_1\mathrm{,}\dots \mathrm{,}{\theta }_m\mathrm{<mo\; stretchy="false">)</mo>}}^T$ the $m$ -component vector of small area means. Also, let $x_i\in {ℝ}^p$  be the $p$ -component vector of auxiliary variables (including the intercept term) for the $i^{\text{th}}$ small area, and $X={\left[x_1\mathrm{,}\dots ,x_m\right]}^T.$ A special case of the HB model given in (1.2)-(1.4) can be expressed as

$$Y|\theta \mathrm{,}\beta \mathrm{,}{\sigma }_v^2~N_m\mathrm{<mo\; stretchy="false">(</mo>}\theta \mathrm{,}D\mathrm{<mo\; stretchy="false">)</mo>,}(2.1)$$

$$\theta |\beta \mathrm{,}{\sigma }_v^2~N_m\mathrm{<mo\; stretchy="false">(</mo>}X\beta \mathrm{,}{\sigma }_v^2{I}_m\mathrm{<mo\; stretchy="false">)</mo>,}(2.2)$$

$$\pi \mathrm{<mo\; stretchy="false">(</mo>}\beta \mathrm{,}{\sigma }_v^2\mathrm{<mo\; stretchy="false">)</mo>}\propto \mathrm{1,}(2.3)$$

where $\beta$ is the $p$ -component regression coefficient vector, ${\sigma }_v^2$ is the model error variance, and $I_m$ is the identity matrix of order $m.$ The uniform prior (2.3) on the model parameters is a popularly used noninformative prior, and the resulting posterior pdf is proper provided that $m>p+2.$ See Berger (1985) and Datta and Smith (2003) for detailed discussion.

The model (2.2) assumes that ${\theta }_i,$ $i=\mathrm{1,}\dots \mathrm{,}m,$ are independently distributed over the small areas with common random effects variance ${\sigma }_v^2.$ In many small area estimation problems, however, the area characteristic of interest is closely related to geographical factors such as population size, ethnicity, age-group and education level. When available covariates do not fully explain such spatial association, the independence and equal variance assumptions of random effects fail, and inference based on the hierarchical model (2.1)-(2.3) may generate unreliable estimates, consequently resulting in erroneous decisions. Figure 2.1 illustrates spatial associations of the 1990 Census 4-person family median incomes of (scaled by $1,000) the $m=49$ states of the U.S., including the District of Columbia. Simulated data with the same Moran’s I value are also displayed for comparison, where Moran’s I is a measure for spatial autocorrelation. The simulated data are generated under the SAR model (defined below) with $\rho =$ 0.8 matching the location and scale with the Census data. The two panels demonstrate the existence of spatial dependence in the 1990 Census 4-person family median incomes. In practice, covariates capable of fully capturing existing spatial variation are not always available, and the problem can be exacerbated if lurking variables exist, as they introduce additional variability that cannot be explained by the independently and identically distributed (i.i.d.) random effects.

To address this issue, we propose to use spatially correlated random effects. Let $W={\left\{w_{ij}\right\}}_{ij},$  $1\le i\mathrm{,}j\le m,$ be the adjacency matrix which plays an important role in capturing spatial dependency. In particular, $w_{ij}=1$ if the $i^{\text{th}}$ and $j^{\text{th}}$ small areas are connected via geographical boundary or through other mechanisms (for example, air traffic), and $w_{ij}=0,$ otherwise. Also, $w_{ii}=0$ for $i=\mathrm{1,}\dots \mathrm{,}m.$ The off-diagonal entries, $w_{ij}\text{’s,}$ need not be binary; they can take other positive values, such as the “length” of the geographical border or volumes of air traffic between the two areas. Since the adjacency matrix $W$ is symmetric, its eigenvalues are real. We denote the $i^{\text{th}}$ largest eigenvalue of $W$ by ${\lambda }_i\mathrm{<mo\; stretchy="false">(</mo>}W\mathrm{<mo\; stretchy="false">)</mo>},$ such that ${\lambda }_m\mathrm{<mo\; stretchy="false">(</mo>}W\mathrm{<mo\; stretchy="false">)</mo>}\le \dots \le {\lambda }_1\mathrm{<mo\; stretchy="false">(</mo>}W\mathrm{<mo\; stretchy="false">)</mo>.}$ Since $W$ is non-null and ${\displaystyle {\sum }_{i\mathrm{<mo>=</mo>1}}^m}{w}_{ii}=0,$ we get as a result that ${\lambda }_m\mathrm{<mo\; stretchy="false">(</mo>}W\mathrm{<mo\; stretchy="false">)</mo>}<0<{\lambda }_1\mathrm{<mo\; stretchy="false">(</mo>}W\mathrm{<mo\; stretchy="false">)</mo>.}$ Let $w_{i\mathrm{.}}={\displaystyle {\sum }_{j\mathrm{<mo>=</mo>}1}^m}{w}_{ij}$ be the sum of the $i^{\text{th}}$ row of $W$ and $L=\text{diag}{\left\{w_{i\mathrm{.}}\right\}}_{i\mathrm{<mo>=</mo>}1}^m.$  Assuming that diagonal elements of $L$ are positive, i.e., all small areas have at least 1 neighboring area, we define $\tilde{W}=L^{-1}W.$ Since $\tilde{W}$ is a row stochastic matrix, all of its eigenvalues are between -1 and 1, with at least one of them is 1. Consequently, ${\lambda }_1\mathrm{<mo\; stretchy="false">(</mo>}\tilde{W}\mathrm{<mo\; stretchy="false">)</mo>}=1.$ Moreover, $\tilde{W}$ and $\text{diag}{\mathrm{<mo>\{</mo>}{w}_{i\mathrm{.}}^{-1/2}\mathrm{<mo>\}</mo>}}_{i\mathrm{<mo>=</mo>}1}^{m}W\text{diag}{\mathrm{<mo>\{</mo>}{w}_{i\mathrm{.}}^{-1/2}\mathrm{<mo>\}</mo>}}_{i\mathrm{<mo>=</mo>}1}^m$ have the same set of eigenvalues, and the latter matrix is symmetric. So all the eigenvalues of $\tilde{W}$ are real, and ${\lambda }_m\mathrm{<mo\; stretchy="false">(</mo>}\tilde{W}\mathrm{<mo\; stretchy="false">)</mo>}$ will be negative. We consider four alternative spatial dependencies associated with random effects, which are represented by the following positive definite precision matrices (excluding the scale parameter ${\sigma }_v^2):$

$$\begin{array}{llll}\text{SAR}\mathrm{<mo>:</mo>}\hfill & {\Omega }_2\mathrm{<mo\; stretchy="false">(</mo>}\rho \mathrm{<mo\; stretchy="false">)</mo>}={\mathrm{<mo\; stretchy="false">(</mo>}{I}_m-\rho \tilde{W}\mathrm{<mo\; stretchy="false">)</mo>}}^T\mathrm{<mo\; stretchy="false">(</mo>}{I}_m-\rho \tilde{W}\mathrm{<mo\; stretchy="false">)</mo>,}\hfill & \rho \in \mathrm{<mo\; stretchy="false">(</mo>}-\mathrm{1,}\mathrm{1<mo\; stretchy="false">)</mo>,}\hfill & (2.4)\hfill \\ \text{SCAR}\mathrm{<mo>:</mo>}\hfill & {\Omega }_3\mathrm{<mo\; stretchy="false">(</mo>}\rho \mathrm{<mo\; stretchy="false">)</mo>}=I_m-\rho W\mathrm{,}\hfill & \rho \in \mathrm{<mo\; stretchy="false">(</mo>}{\lambda }_m{\mathrm{<mo\; stretchy="false">(</mo>}W\mathrm{<mo\; stretchy="false">)</mo>}}^{-1}\mathrm{,}{\lambda }_1{\mathrm{<mo\; stretchy="false">(</mo>}W\mathrm{<mo\; stretchy="false">)</mo>}}^{-1}\mathrm{<mo\; stretchy="false">)</mo>,}\hfill & (2.5)\hfill \\ \text{CAR}\mathrm{<mo>:</mo>}\hfill & {\Omega }_4\mathrm{<mo\; stretchy="false">(</mo>}\rho \mathrm{<mo\; stretchy="false">)</mo>}=L-\rho W\mathrm{,}\hfill & \rho \in \mathrm{<mo\; stretchy="false">(</mo>}{\lambda }_m{\mathrm{<mo\; stretchy="false">(</mo>}\tilde{W}\mathrm{<mo\; stretchy="false">)</mo>}}^{-1}\mathrm{,}{\lambda }_1{\mathrm{<mo\; stretchy="false">(</mo>}\tilde{W}\mathrm{<mo\; stretchy="false">)</mo>}}^{-1}\mathrm{<mo\; stretchy="false">)</mo>,}\hfill & (2.6)\hfill \\ \text{LCAR}\mathrm{<mo>:</mo>}\hfill & {\Omega }_5\mathrm{<mo\; stretchy="false">(</mo>}\rho \mathrm{<mo\; stretchy="false">)</mo>}=\rho R+\mathrm{<mo\; stretchy="false">(</mo>1}-\rho \mathrm{<mo\; stretchy="false">)</mo>}{I}_m\mathrm{,}\hfill & \rho \in \mathrm{<mo\; stretchy="false">(</mo>0,}\mathrm{1<mo\; stretchy="false">)</mo>,}\hfill & (2.7)\hfill \end{array}$$

where $\rho$ is the spatial dependence parameter that represents the strength of spatial dependence (Hodges, 2019, Chapter 5.2) and $R$ is defined as $R={\Omega }_4\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}=L-W.$ Since the eigenvalues of $I_m-\tilde{W}$ are between 0 (the smallest eigenvalue) and $1-{\lambda }_m\mathrm{<mo\; stretchy="false">(</mo>}\tilde{W}\mathrm{<mo\; stretchy="false">)</mo>}$ (the largest eigenvalue, $>1),$ the matrix $R$ is nonnegative definite. Each precision matrix is guaranteed to be positive definite as long as $\rho$ is in the range specified in the respective definition.

The adjacency matrix $\tilde{W}$ of the simultaneous autoregressive (SAR) model (Whittle, 1954) is row-normalized so that $\rho$ can vary from -1 to 1 while preserving the positive definiteness (Banerjee, Carlin and Gelfand, 2003, Chapter 4.4). The model (2.5) is a simple version of conditional autoregressive (CAR) model (Rao and Molina, 2015, Chapter 9.6.2), where diagonal entries of the precision matrix are all equal to one. Even though the diagonal elements of a precision matrix are all equal, the diagonal elements of the inverse may not be all equal, leading to heteroscedasticity of random effects. We call this model the simple conditional autoregressive (SCAR) model. The model (2.6) is widely used conditional autoregressive model (CAR; Banerjee et al., 2003; Besag and Kooperberg, 1995; You and Zhou, 2011), where diagonal entries of the precision matrix are the number of neighborhoods of the corresponding area. The upper limit of $\rho$ is ${\lambda }_1{\mathrm{<mo\; stretchy="false">(</mo>}\tilde{W}\mathrm{<mo\; stretchy="false">)</mo>}}^{-1}=1,$ and in the case of $\rho =1,$ the model with ${\Omega }_4\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}$ is referred to as the intrinsic autoregressive (IAR) model (Banerjee et al., 2003, Chapter 4.3). The model (2.7) is a conditional autoregressive model, which we call Leroux’s conditional autoregressive (LCAR), whose precision matrix is given by the convex combination of $R={\Omega }_4\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}$ and $I_m.$ This model has been considered by Leroux, Lei and Breslow (2000); MacNab (2003); You and Zhou (2011), where the $i^{\text{th}}$ diagonal element of $R$ is the number of neighborhoods of the $i^{\text{th}}$ small area, and the ${\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{,}j\mathrm{<mo\; stretchy="false">)</mo>}}^{\text{th}}$ off-diagonal element is -1 if the $i^{\text{th}}$ and the $j^{\text{th}}$ small areas are connected and 0 otherwise.

The conditional autoregressive models, SCAR, CAR, and LCAR, assume that ${\theta }_i$ depends only on neighboring small area means. In other words, ${\theta }_i$ is correlated with ${\theta }_j\text{’s},$ $j\ne i,$ only through the means of surrounding areas. On the contrary, the SAR model assumes that ${\theta }_i$ is dependent on all other ${\theta }_j$ concurrently, $j\ne i,$ but has stronger (weaker) correlations for neighboring (remote) areas. The independent FH model can be viewed as a special case of the SAR, SCAR, or LCAR model with $\rho =0.$ For notational convenience, we include the independent FH model as part of our model by taking its precision matrix ${\Omega }_1\mathrm{<mo\; stretchy="false">(</mo>}\rho \mathrm{<mo\; stretchy="false">)</mo>}=I_m,$ although it is free from $\rho .$

We consider the following HB spatial models incorporating the five spatial dependencies defined in (2.4)-(2.7):

$$Y|\theta \mathrm{,}\beta \mathrm{,}{\sigma }_v^2\mathrm{,}\rho ~N_m\mathrm{<mo\; stretchy="false">(</mo>}\theta \mathrm{,}D\mathrm{<mo\; stretchy="false">)</mo>,}(2.8)$$

$$\theta |\beta \mathrm{,}{\sigma }_v^2\mathrm{,}\rho ~N_m\mathrm{<mo\; stretchy="false">(</mo>}X\beta \mathrm{,}{\sigma }_v^2{\left\{{\Omega }_k\mathrm{<mo\; stretchy="false">(</mo>}\rho \mathrm{<mo\; stretchy="false">)</mo>}\right\}}^{-1}\mathrm{<mo\; stretchy="false">)</mo>,}k=\mathrm{1,}\dots \mathrm{,}\mathrm{5,}(2.9)$$

where ${\sigma }_v^2$ is the model scale parameter, $g\mathrm{<mo\; stretchy="false">(</mo>}{\sigma }_v^2\mathrm{<mo\; stretchy="false">)</mo>}$ and $h\mathrm{<mo\; stretchy="false">(</mo>}\rho \mathrm{<mo\; stretchy="false">)</mo>}$ are suitable functions of ${\sigma }_v^2$ and $\rho ,$ $l_k$ and $u_k$ are the lower and upper bounds of $\rho$ under the $k^{\text{th}}$ model. We avoid the term “model error variance” for ${\sigma }_v^2$ as diagonal entries of ${\Omega }_k\mathrm{<mo\; stretchy="false">(</mo>}\rho \mathrm{<mo\; stretchy="false">)</mo>}$ vary across small areas and are not necessarily all one.

2.2   Estimation of population means of non-sampled small areas

In this section, we consider the case when, in the survey, there are several non-sampled small areas that have no direct estimates. In many applications, limited resources frequently preclude the inclusion of many subpopulations in the sample, resulting in non-sampled small areas. In this section, we consider the case when, in the survey, there are several non-sampled small areas that have no direct estimates. In many applications, limited resources frequently preclude the inclusion of many subpopulations in the sample, resulting in non-sampled small areas. Non-sampled small areas are sometimes referred to as misaligned areas (Trevisani and Gelfand, 2013) when they arise from domain misalignment between the direct estimate and auxiliary variables. For any of these non-sampled areas, the prediction of its mean from any non-spatial model is only based on its synthetic estimator. We propose to exploit spatial dependencies in predicting area means of non-sampled small areas. The predictions of the proposed models are obtained by modifying its synthetic estimator, using the vector of regression residuals, with more emphasis on the regression residuals of the neighboring areas.

Without loss of generality, let there be $m_1$ non-sampled small areas and $Y_{m_1+1}\mathrm{,}\dots \mathrm{,}{Y}_m$ be the direct estimates of the $m_2=m-m_1$ sampled small areas. Based on the direct estimates of $m_2$ sampled areas, we consider the following HB models:

$$Y_{\mathrm{<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}}|\theta \mathrm{,}\beta \mathrm{,}{\sigma }_v^2\mathrm{,}\rho ~N_{m_2}({\theta }_{\mathrm{<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}}\mathrm{,}{D}_{\mathrm{<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}})\mathrm{,}(2.11)$$

$$\theta |\beta \mathrm{,}{\sigma }_v^2\mathrm{,}\rho ~N_m\left(X\beta \mathrm{,}{\sigma }_v^2{\left\{{\Omega }_k\mathrm{<mo\; stretchy="false">(</mo>}\rho \mathrm{<mo\; stretchy="false">)</mo>}\right\}}^{-1}\right)\mathrm{,}k=\mathrm{1,}\dots \mathrm{,}\mathrm{5,}(2.12)$$

where $Y_{\mathrm{<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}}={\mathrm{<mo\; stretchy="false">(</mo>}{Y}_{m_1+1}\mathrm{,}\dots \mathrm{,}{Y}_m\mathrm{<mo\; stretchy="false">)</mo>}}^T,$ $D_{\mathrm{<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}}=\text{diag}{\mathrm{<mo>\{</mo>}{D}_i\mathrm{<mo>\}</mo>}}_{i\mathrm{<mo>=</mo>}{m}_1+1}^m,$ and ${\theta }_{\mathrm{<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}}={\mathrm{<mo\; stretchy="false">(</mo>}{\theta }_{m_1+1}\mathrm{,}\dots \mathrm{,}{\theta }_m\mathrm{<mo\; stretchy="false">)</mo>}}^T,$ which is the subvector of $\theta$ corresponding to the sampled areas.

2.3   Propriety of the posterior distributions

In this section, we establish propriety of the posterior distributions of spatial small area models given in (2.8)-(2.10) and (2.11)-(2.13). Let $I(\cdot )$  be the indicator function taking the value 1 when its argument is true and 0 otherwise. We first provide general conditions for the posterior propriety of the proposed models.

Theorem 1. For all the HB spatial models given in (2.8)-(2.10) and (2.11)-(2.13), the posterior probability density functions are proper if the following conditions hold for some positive constant $c>0:$

• (a)   ${\displaystyle {\int }_0^{\infty }}g\mathrm{<mo\; stretchy="false">(</mo>}{\sigma }_v^2\mathrm{<mo\; stretchy="false">)</mo>}I\mathrm{<mo\; stretchy="false">(</mo>}{\sigma }_v^2\le c\mathrm{<mo\; stretchy="false">)</mo>}d{\sigma }_v^2<\infty .$
• (b)  ${\displaystyle {\int }_0^{\infty }}{\mathrm{<mo\; stretchy="false">(</mo>}{\sigma }_v^2\mathrm{<mo\; stretchy="false">)</mo>}}^{-\mathrm{<mo\; stretchy="false">(</mo>}{m}^{*}-p\mathrm{<mo\; stretchy="false">)</mo>}/2}g\mathrm{<mo\; stretchy="false">(</mo>}{\sigma }_v^2\mathrm{<mo\; stretchy="false">)</mo>}I\mathrm{<mo\; stretchy="false">(</mo>}{\sigma }_v^2>c\mathrm{<mo\; stretchy="false">)</mo>}d{\sigma }_v^2<\infty .$
• (c)   ${\displaystyle {\int }_{l_k}^{u_k}}h\mathrm{<mo\; stretchy="false">(</mo>}\rho \mathrm{<mo\; stretchy="false">)</mo>}d\rho <\infty ,$

where $m^{*}=m$ for (2.8)-(2.10), and $m^{*}=m-m_1$ for (2.11)-(2.13).

If $g(\cdot )$  is a proper pdf, then (a) holds true automatically, and (b) is satisfied if $m^{*}\ge p.$ The condition $m^{*}\ge p$ is obvious since at least $p$ observations are needed to estimate $p$ components of $\beta$ when no substantive information about it is available. Also, any bounded function of $\rho$ satisfies $\mathrm{<mo\; stretchy="false">(</mo>}c\mathrm{<mo\; stretchy="false">)</mo>}$ in Theorem 1 as their supports are all bounded. In particular, under the popular family of noninformative priors

the posterior pdfs are proper under the following conditions.

Corollary 1. For any of the HB spatial models given in (2.8)-(2.9) and (2.11)-(2.12) with the prior in (2.14), the posterior pdf is proper as long as $\alpha <1$ and $m^{*}>p+2-2\alpha .$

For the uniform prior with $\alpha =0$ (which will be used in this paper), the propriety of the posterior distributions for models (2.8)-(2.9) are guaranteed as long as the number of small areas is greater than $p+2.$ For the models incorporating non-sampled areas given in (2.11)-(2.12), the second condition of Corollary 1.1 becomes $m-m_1>p+2,$ and thus, the posterior pdfs are proper as long as the number of non-sampled areas is fewer than $m-p-2,$ or at least $p+3$ areas have sample.

• Previous
• Article main page
• Next