Global-Local Priors for Spatial Small Area Estimation

1. Introduction

With a modest beginning in the seventies, small area estimation has now emerged as a discipline of its own, catering to the needs of both theoretical and applied survey researchers. To our knowledge, the term ‘small area’ was first coined by Gonzalez and Hoza^[10] in the context of finding unemployment and housing estimates. A small area may refer to small geographical areas such as counties, census tracts, school districts and others. But it may also refer to small domains cross classified by age, sex, ethnicity and other demographic characteristics. In these days, small area estimates have become virally indispensable for policy formulation by both the public and private sectors. A notable example is Small Area Income and Poverty Estimation (SAIPE) by the United States Bureau of the Census.

Direct survey estimates of small areas are usually unreliable owing to inadequate sample size, resulting in large standard errors and coefficients of variation. The reason behind this is that the original survey was designed to attain prescribed accuracy at a higher level of aggregation than those of small or local areas. Owing to limited resources, there cannot be a new survey, and data from the original survey need to be used to obtain small area estimates, and such data will typically be unable to attain the desired level of precision. This necessitates ‘borrowing strength from the ensemble’, or in other words, link the different small areas through a model.

While classical small area estimates such as ‘synthetic’ or ‘composite’ estimates did not make explicit use of models, they could often be justified using a model-based approach. A major breakthrough came in the classic article of Fay and Herriot^[8] who considered small area estimation of per capita income. In essence, the Fay-Herriot model is a linear mixed effects model with random area effects. The said model has been widely used for more than four decades, and it is anticipated that this will continue in the years ahead.

Over the years, there have been several generalizations of the Fay-Herriot model to address issues related to robustness, departures from normality, techniques such as jackknife and bootstrap to estimate mean squared errors of small area estimates, handling binary and count data and many others. The wide variety of research has enriched the small area literature to a very large extent. A very comprehensive account appeared in Rao and Molina^[22] which captured most of the recent work up until that time. However, a phenomenal advancement of the subject has happened since then, making it quite a challenge to keep up with this growth.

The Fay-Herriot model uses a single area level variance component across all small areas. This works well for small or moderate number of small areas. But the assumption is likely to be violated when the number of small areas is very large, for example when the small areas are more than 3,000 counties of the United States. With areas varying considerably in size, it is impossible to visualize homogeneity across all small areas.

What we want to consider in this article is that the area level variances are products of two components. One is a global parameter which is the same across all areas as in the Fay-Herriot model, but there are also local parameters which vary across all areas. The global parameter is designed to shrink all the random effects towards zero so that the small random effects are captured, while the local parameters try to neutralize the shrinkage effect for areas that need large random effects. The degree of this neutralizing effects is closely related with the tails of the priors on the local parameters. With appropriately heavy-tailed priors, both small and large random effects can be well-captured. This was implemented in Tang et al.^[26].

The global-local shrinkage priors were introduced in a series of articles (Carvalho et al.^[4]; Polson and Scott).^[17] They have been extended into a richer class. Some recent inventions are the three parameter beta normal priors Armagan et al.,^[1] which includes the now famous horseshoe prior Carvelho et al.^[4] and the generalized double Pareto prior Armagan et al.,^[2].

The objective of this article is to extend the work of Tang et al.^[26] to spatial small area models. The spatial aspect is quite natural when, for example, one thinks of adjoining counties or sometimes even areas that are structurally very similar. This requires modifying the priors proposed in Tang et al.^[26] in order to incorporate the spatial structure.

Spatial models have started appearing in the small area literature. Among others, we may refer to Petrucci and Salvati^[16], Pratesi and Salvati^[21], You and Zhou,^[27] Schmid and Mnnich^[24], Mercer et al.^[12], Porter et al.^[20] and Porter et al.^[19]. However, the introduction of global-local priors in this context seems to be new.

The work of Tang et al.^[26] was motivated by some previous work of Datta et al.^[6], Datta and Mandal^[7] and Chakraborty et al.^[5] The objective of Datta et al.^[6] was to investigate whether one needed a random area effect for all small areas, or whether a fixed effect model should suffice for some of the areas. To this end, these authors suggested a test-based approach where the null hypothesis is that the common random effect variance is zero. The test was based on a discrepancy statistic measuring the lack of fit of the multiple regression model of small area means on certain covariates. This work was followed in subsequent articles of Molina et al.^[13] and Morales et al.^[14]

The above approach works well when the number of small areas is moderately large, but, as it often happens in practice, the number of small areas is very large, for example, when one considers all counties in the United States. In such situations, the null hypothesis of no random effects is very likely to be rejected, owing to significant departure of direct estimates from the regression estimates caused by a few large residuals.

In order to rectify the above problem, Datta and Mandal^[7] proposed to model random effects with spike-and-slab priors. These priors put a positive mass at zero for the area level variance component (the spike part), but also a normal prior (the slab part) to account for variance components which differ significantly from zero. In areas where random effects are not needed, the variance component can be very small, while, in areas where this is not the case, the variance component can be large. These priors use a single variance component which can be zero with a positive probability. A slight modification was provided later in Chakraborty et al.^[5] which allows two variance components, one large and one small. The approach of Tang et al.^[26] seems more informative, as it allows multiple variance components and provides posterior estimates of the same.

The rest of the article is organized as follows. Section 2 of this article introduces the spatial model, proves the propriety of the resulting posterior under some conditions, and discusses the implementation of the proposed method. Some simulation results are provided in Section 3, while Section 4 contains analyses of real datasets. Some final remarks are made in Section 5.

2. The Model

Consider m small areas. For i = 1, … , m , let y i denote the direct survey estimate for the small area mean θ i of area i . We assume

y i | θ i ∼ N ( θ i , V i ) and θ i = x i ⊤ β + u i ,

(2.1)

where V i is the known variance of the sampling error, x i is a p -dimensional covariate vector, β is the responding coefficient vector, and u i is a random effect describing the variation in θ i not explained by x i . We will use the notations θ = ( θ 1 , … , θ m ) ⊤ , u = ( u 1 , … , u m ) ⊤ , V = diag { V 1 , … , V m } , y = ( y 1 , … , y m ) ⊤ and X = ( x 1 , ... , x m ) ⊤ and assume rank( X ) = p ≤ m . Model (2.1) is consistent with the classic Fay-Herriot (FH) model where the random effects are further assumed to follow independent normal distributions with a common variance as

u | σ 2 ∼ N ( 0 , σ 2 I ).

(2.2)

In our model, we assume

u | ρ , λ , τ ∼ N ( 0 , τ 2 Λ ( D − ρ W ) − 1 Λ ),

(2.3)

where λ = ( λ 1 , … , λ m ) , Λ = d i a g { λ 1 , … , λ m } , W is an m × m adjacency matrix whose i j -th element, w i j , is 1 if areas i and j are neighbors and 0 otherwise, and D is a diagonal matrix with the diagonal elements being the row sums of W . This model is able to account for both the spatial dependence and the wide variation in the random effects. To better see this, we express (2.3) as

u i = τ λ i u ⌣ i , i = 1, … , m ;

(2.4)

u ⌣ ∼ N ( 0, ( D − ρ W ) − 1 ),

(2.5)

where u ⌣ = ( u ⌣ 1 , … , u ⌣ m ) ⊤ . In (2.4), parameters τ and λ i play similar roles as the global parameter and local parameter, respectively, in the global-local (GL) model proposed in Tang et al.^[26] The global parameter τ shared across areas controls the overall scale of the random effects while the local parameter λ i , which is unique for each area, controls the scale of individual random effects. In fact, (2.4) is equivalent to the GL model if u ⌣ ∼ N ( 0, I ) . If both small and large random effects (in absolute value) are expected, then we need a small τ to constrain the overall scale of random effects to be small. For those areas with large random effects, we need a large λ i to compensate the effect of τ .

The spatial dependence among the random effects is introduced through (2.5) which imposes a Conditional Auto-Regressive (CAR) model Besag^[3] for u ⌣ . The CAR model is often used to introduce dependence among neighbors in spatial statistics. It is can be alternatively formulated as

u ⌣ i | u ⌣ − i ∼ N ρ ∑ j = 1 m w i j u ⌣ j ∑ j = 1 m w i j , 1 ∑ j = 1 m w i j , i = 1, … , m ,

(2.6)

where the dependence of u ⌣ i on the corresponding components in the neighboring areas is more evident. The parameter ρ ∈ ( 0,1 ) controls the strength of spatial dependence. A larger ρ produces stronger dependence. If ρ = 0 , then u ⌣ 1 , … , u ⌣ m are independent. Since our model is a combination of the GL model and the spatial CAR model, we call our model Global-Local Spatial (GLSP) model.

Under our GLSP model described in (2.3),

E ( θ | y , β , λ , τ , ρ ) = ( I − Q ) y + Q X β = y + Q ( y − X β ),

(2.7)

where Q = ( V − 1 + U − 1 ) − 1 U − 1 and U = τ 2 Λ ( D − ρ W ) − 1 Λ . Similar to the FH and GL models, the conditional posterior mean of θ can be expressed as a weighted average of the direct estimator y and the systematic estimator X β . It shrinks y towards X β with a matrix shrinkage factor Q . However, Q is not a diagonal matrix unlike the one in the FH and GL models. The conditional posterior mean of θ i depends on the data in not only area i but also its neighboring areas due to the spatial dependence in the random effects.

Prior specification. To complete the model, we specify priors on the parameters β , λ , τ 2 , ρ . We consider independent priors on these parameters

π ( β , λ , τ 2 , ρ ) = ∏ i = 1 m π λ i π ( β ) π ( τ 2 ) π ( ρ ).

We choose π ( β ) ∝ 1 following the convention in Bayesian SAE models. A beta prior Beta ( c 0 , d 0 ) is imposed on ρ and an inverse gamma prior IG ( a 0 , b 0 ) is imposed on τ 2 . We use a 0 = b 0 = 10 − 10 and c 0 = d 0 = 1 in the simulation experiments and real data examples presented later. Recently, the half-Cauchy prior has been recommended as a default prior for scale parameters as opposed to the inverse gamma priors IG ( a 0 , b 0 ) (Gelman^[9]; Polson and Scott^[18]) as the performance with the latter prior may be sensitive to the choice of a 0 and b 0 . Based on our exploration and also the results in Tang et al.^[26], the half Cauchy prior and the chosen inverse gamma prior do not produce significantly different results. Therefore, we choose to use the inverse gamma prior for τ 2 for convenience. Since λ i ’s in the GLSP model and in the GL model play similar roles, we borrow the priors used in Tang et al.^[26] here. Table 1 of Tang et al.^[26] provides a list of possible choices of the prior on λ i 2 . These priors can be divided into two groups, one with exponential tails and other one with polynomial tails. We consider the Laplace (LA) prior π ( λ i ) = exp ( − λ i ) and the famous horseshoe (HS) prior (Carvalho et al.^[4]) π ( λ i ) = ( 2 / π ) ( 1 + λ i 2 ) − 1 as the representatives of the two groups in the GLSP model.

The joint posterior is then given by

π u ˜ , β , λ , τ 2 , ρ ∣ y ∝ exp − 1 2 ( y − X β − Λ u ˜ ) ⊤ V − 1 ( y − X β − Λ u ˜ ) × τ 2 − m 2 | D − ρ W | 1 2 exp − 1 2 τ 2 u ˜ ⊤ ( D − ρ W ) u ˜ × τ 2 − a 0 − 1 e − b 0 τ 2 ρ c 0 − 1 ( 1 − ρ ) d 0 − 1 ∏ i = 1 m π λ i ,

(2.8)

OPEN IN VIEWER

Table 1.

DIC for various models fitted to the two datasets.

Model	FH	GL-HS	GL-LA	SP	GLSP-HS	GLSP-LA
Data 1	273.29	273.09	275.92	276.77	268.70	274.96
Data 2	−2093.14	−2068.44	−2097.30	−2089.30	−2065.26	−2104.44

where u ˜ = τ u ⌣ . Although we use an improper prior for β , the posterior propriety is guaranteed under mild conditions as the theorem below shows. The proof of the theorem is deferred to the appendix.

Theorem 1. The posterior derived in (2.8) is proper if X is of full column rank and π ( λ i ) is proper.

Computation. The implementation of the proposed Bayesian method is achieved via Markov chain Monte Carlo (MCMC). More specifically, we use a Gibbs sampler to obtain posterior samples of the parameters and random effects. To this end, we write the full conditionals as

π ( β | y , u ˜ , λ , τ 2 , ρ ) = π ( β | y , u ˜ , λ ) ∝ exp − 1 2 ( β − μ β ) ⊤ Σ β − 1 ( β − μ β ) , where μ β = Σ β X V − 1 z , Σ β = ( X ⊤ V − 1 X ) − 1 , and z = y − Λ u ˜ .

The full conditionals of β , u ˜ , and τ 2 suggest that these components can be sampled directly from N ( μ β , Σ β ) , N ( μ u , Σ u ) , and I G ( a 0 + m 2 , b 0 + 1 2 ⊤ ( D − ρ W ) u ˜ ) , respectively. Other components, λ and ρ , cannot be sampled directly. We use a Metropolis-Hastings step to draw samples for ρ . The proposal distribution is the same as the prior distribution B e t a ( c 0 , d 0 ) . The components of λ can be sampled independently. The full conditional of λ i can be expressed as

π λ i | y i , β , u ˜ i ∝ exp − 1 2 σ λ 2 i λ i − μ λ i 2 π ( λ i ),

(2.9)

where μ λ i = ( y i − x i ⊤ β ) / u ˜ i and σ λ i 2 = V i / u ˜ i 2 . The first factor on the right-hand side of (2.9) is in the form of a normal density but with the constraint λ i > 0 . Following Ray et al.,^[23] we approximate the indicator function 1 ( 0, ∞ ) ( λ i ) by a sigmoid function 1 / ( 1 + exp ( − η λ i ) ) with a large constant η > 0 and use elliptical sliced sampling (Hahn et al.^[11] 2019; Murray et al.^[15]) to draw samples for λ i .

In this section, we demonstrate the performance of the proposed model in estimating the small area means using simulations. In our experiments, datasets of m small areas are generated according to model (2.1) with m = 100,200,500 , p = 2 , and β = ( 10,1 ) ⊤ . The covariate matrix X = ( x 1 , … , x m ) ⊤ consists of a column of ones and a randomly generated column whose elements are drawn independently from the standard normal distribution. The following models are used for generating the random effects u :

u 1 , … , u m ∼ i i d N ( 0,4 ) .

The random effects in Model 1 are independent while those in Model 2 are spatially dependent. In Models 3–5, various levels of sparsity is introduced on top of spatial dependence. The five models are labeled as Normal, Spatial, Sparse Spatial (SSP) 0.8, SSP 0.5, and SSP 0.2, respectively.

Among the two choices of W , the first one has a simpler structure. The elements w i j = 1 if and only if | i − j | ≤ 1 . The second choice considers a more complex structure. In this choice, we assume the m areas are represented by the nodes of a 10 × ( m / 10 ) regular square grid. Areas that are directly connected by an edge are treated as adjacent areas.

For each setting, 100 datasets are generated. The proposed model with two choices of the local priors, HS and LA, is fitted to each dataset. The two models are denoted by GLSP-HS and GLSP-LA, respectively. We also fit the proposed model with the local parameter λ i 2 set to 1, so that only spatial dependence is considered. Besides, FH and the GL model without spatial dependence are also fit for comparison. The choices of the local priors in the non-spatial GL models are also HS and LA. The two resulting GL models are denoted by GL-HS and GL-LA, respectively. For each model, the corresponding MCMC algorithm is run for 20,000 iterations. The first half is discarded as burn-in. The posterior sample means computed from the remaining 10,000 iterations are used as point estimates of the small area means. The estimation performance is evaluated using the average absolute deviation (AAD), average squared deviation (ASD), average absolute relative deviation (ARB), and average squared relative deviation (ASRB). They are defined as

AAD = 1 m ∑ i = 1 m | θ ^ i − θ i |, ASD = 1 m ∑ i = 1 m ( θ ^ i − θ i ) 2 ,

ARB = 1 m ∑ i = 1 m | ( θ ^ i − θ i ) / θ i |, ASRB = 1 m ∑ i = 1 m ( θ ^ i − θ i ) 2 / θ i 2 .

The averages of the deviation measures across 100 datasets for different settings with the two choices of W are presented in Figures 1 and 2. As the figures show, the proposed GLSP models usually outperform their GL counterparts when spatial dependence exists, and they have similar performances when spatial dependence does not exist.

OPEN IN VIEWER

Figure 1.

Average deviation measures of estimated small area means under different settings for the first choice of W . The horizontal axis is the number of small areas.

OPEN IN VIEWER Figure 2.

Average deviation measures of estimated small area means under different settings for the second choice of W . The horizontal axis is the number of small areas.

According to the results in the first two columns of the figures, FH and SP models usually have the best performance for data generated from the Normal and Spatial settings, respectively. This is not surprising as the fitted model matches the data generation model. Although the proposed GLSP models are usually not the best in these settings, the performance of GLSP-LA is often close to the best. It is even slightly better than the SP model in the Spatial setting presented in Figure 1.

The last three columns in Figures 1 and 2 show that if the local priors are chosen appropriately, the proposed GLSP models often outperform other models when both spatial dependence and moderate to high levels of sparsity (proportion of zero random effects) exist in the generated random effects. Overall, GLSP-LA tends to have a better performance than GLSP-HS regardless of the deviation measures. However, in the high sparsity level SSP setting (SSP 0.2), GLSP-HS performs better in terms of AAD and ARB. This seems to suggest GLSP-HS tends to have a better performance when the sparsity level is higher. One possible reason for GLSP-HS not showing advantages over GLSP-LA in terms of ASD and ASRB in this setting is that the HS prior shrinks small random effects towards zero more aggressively than the LA prior. This effect causes relatively large discrepancies between the estimated and true small area means in a few areas. These discrepancies stand out under the measures involving squared errors. Similar phenomena are also observed when comparing the GLSP-LA and the SP models in the SSP 0.5 setting.

In this section, we apply the proposed model to two real datasets that have been analyzed by Datta and Mandal^[7] and Tang et al.^[26]. The first dataset provides the direct estimates of child (age 5–17) poverty ratio obtained from 1999 current population survey for 50 states and the District of Columbia of the United States ( m = 51 ). Three covariates, the number of child exemptions, Internal Revenue Service no-filer rate, and the residuals from fitting a model for the 1989 census poverty data on the previous two covariates, are also included in the dataset. The second dataset includes 5-year (2007–2011) American Community Survey (ACS) estimates of overall poverty rates for 3,141 counties in the United States. For simplicity, we only consider 414 counties in the West Census Region excluding Alaska and Hawaii. Same as in Tang et al.^[26], foodstamp participation rate is used as the covariate.

We fit FH, GL-HS, GL-LA, SP, GLSP-HS, and GLSP-LA models to both datasets. In the GLSP models, the adjacency matrix W is constructed such that two areas are neighbors if they share boundaries. The MCMC algorithms are run for 100,000 iterations. The samples from the first 50,000 iterations are discarded as the burn-in period. For the remaining iterations, samples are kept for every ten iterations. Following Tang et al.^[26], we use the Deviance Information Criterion (DIC) proposed in Spiegelhalter et al.^[25] to select the most appropriate model for each of the two datasets. Among several models, the one with the smallest DIC is selected as the most appropriate. Let θ ( g ) = ( θ 1 ( g ) , … , θ m ( g ) ), g = 1, … , G , denote the posterior samples of the small area means based on an aforementioned model. The DIC of the model is estimated by

DIC ^ = 2 G ∑ g = 1 G ∑ i = 1 m ( y i − θ i ( g ) ) 2 V i − ∑ i = 1 m ( y i − θ ^ i ) 2 V i ,

where θ ^ i = 1 G ∑ g = 1 G θ i ( g ) is the posterior sample mean of θ i . According to the DIC results presented in Table 1, GLSP models are selected for both datasets but with different priors for the local parameters λ i 2 . The GLSP-HS model is selected for the state-level data, while the GLSP-LA model is selected for the county-level data. Previous analyses of the state-level dataset showed that only Massachusetts has a significantly nonzero random effect (Datta and Mandal^[7]; Tang et al.^[26]), which indicates the high sparsity level of random effects. As a result, GLSP-HS is preferred over GLSP-LA. Tang et al.^[26] considered the GL models without spatial dependence for the county-level data with all 3,141 counties and showed that the LA prior is more appropriate for the dataset than the HS prior. Our results are consistent with theirs.

In the following, we take a closer look at the results for the county-level dataset as the dataset contains more small areas. In Figure 3, we compare the results from the GL-LA and the GLSP-LA models. The left panel shows the small area mean estimates from the two models are very close. The right panel shows the potential benefits of GLSP models for reducing the variation of the estimation. Although the GLSP-LA model does not produce a smaller coefficient of variation (CV) for all small areas, it reduces CV in 50% of the counties and often by a large amount while producing a slightly higher CV for the remaining 50% of counties. The reduction can be as high as 0.17 and is more than 0.04 in 15 counties, while the increase is at most 0.042.

OPEN IN VIEWER

Figure 3.

Left: Estimated poverty rates from GL-LA and GLSP-LA. Right: Coefficients of Variation from GL-LA and GLSP-LA.

Finally, we present a map of small mean estimates from the GLSP-LA models in Figure 4. Among the counties we consider, Douglas County in Colorado has the lowest poverty ratio 0.036, while Boise County in Idaho has the highest poverty ratio 0.347.

OPEN IN VIEWER Figure 4.

Estimated poverty rates of 414 counties in West Census Region.

The article considers GL priors in the context of small area estimation where the area level random effects exhibit a spatial structure. The proposed GLSP model combines the popular CAR model and the GL model of Tang et al.^[26] to accommodate both the natural spatial dependence among geographically connected areas and the wide variation of random effects for a large number of small areas. The simulation experiments show that the GLSP model can outperform the GL model under various settings with spatial dependence. In the real data analysis, we show that the GLSP model can significantly reduce the CV for some areas without substantially increasing the CV for other areas.

In our MCMC algorithm for drawing posterior samples for the proposed model, the random effects u ˜ are sampled jointly from a multivariate normal distribution. This step requires computing the variance–covariate matrix Σ u in every iteration. Such computation is time-consuming when a large number of small areas are involved. For example, fitting the SP, GLSP-LA, or GLSP-HS model on the county level dataset considered in Section 4 takes about 6 hours on a personal laptop with 8 GB memory and 2.3 GHz processor while fitting the FH or the non-spatial GL model takes about 20 and 280 seconds, respectively. As a remedy to the heavy computation cost, one could consider using the alternative formulation (2.6) to perform elementwise sampling of u ˜ .

In terms of future work, it will be interesting to see how alternative spatial models perform in a similar context. Also of interest will be a generalization to spatial–temporal models.

Appendix: Proof of Theorem 1

Proof. Let X * = V − 1 2 X , P X * = ( X * T X * ) − 1 ( X * ) T β ^ = ( X * T X * ) − 1 X * T ( θ − Λ u ˜ . , We begin with the identity

∑ i = 1 m y i − x i ⊤ β − λ i u ˜ i 2 / V i = β ⊤ X ∗ X ∗ β − 2 β ⊤ X ∗ ⊤ ( y − Λ u ˜ ) + ( y − Λ u ˜ ) ⊤ ( y − Λ u ˜ ) = ( β − β ^ ) ⊤ ( β − β ^ ) + ( y − Λ u ˜ ) ⊤ V − 1 2 I − P X ∗ V − 1 2 ( y − Λ u ˜ )

(A.1)

Integrating π ( u ˜ , β , λ , τ 2 , ρ | y ) with respect to β , it follows from (2.8) and (A.1),

π ( u ˜ , λ , τ 2 , ρ | y ) = K exp − 1 2 ( y − Λ u ˜ ) ⊤ V − 1 / 2 ( I − P X * ) V − 1 / 2 ( y − Λ u ˜ ) × ( τ 2 ) − m 2 | D − ρ W | 1 2 exp − 1 2 τ 2 u ˜ ⊤ ( D − ρ W ) u ˜ π ( τ 2 ) π ( ρ ) ∏ i = 1 m π ( λ i ) ≤ K ( τ 2 ) − m 2 | D − ρ W | 1 2 exp − 1 2 τ 2 u ˜ ⊤ ( D − ρ W ) u ˜ π ( τ 2 ) π ( ρ ) ∏ i = 1 m π ( λ i )

(A.2)

where K > 0 is a generic constant. Since the priors on ρ , τ 2 and λ i are proper, integrating with respect to u ˜ and these parameters gives finiteness of the integral in the right hand side of (A.2).