Bayesian spatial models for estimating means of sampled and non-sampled small areas
Section 4. A simulation study

• Previous
• Next

In this section, we compare prediction performances of the independent FH model and the four spatial models in the absence of informative covariates with multiple non-sampled areas. Excluding Hawaii and Alaska, we consider contiguous $m=49$  states of the U.S., including the District of Columbia. To evaluate the quality of prediction in the absence of direct estimates, we do not simulate direct estimates of randomly chosen $m_1=⎣0.15m⎦\mathrm{<mo>=</mo>}7$ states. These areas are Delaware, Massachusetts, Michigan, Nebraska, Rhode Island, South Dakota, and Texas. This results in $m_2=42$ areas which have direct estimates.

To make simulation settings realistic, we mimic the 1989 4-person family median income (median income) data described in Section 5. We generate replicated datasets so that Moran’s I values of each replicated small area means, ${\theta }_1\mathrm{,}\dots \mathrm{,}{\theta }_m,$ are approximately centered around 0.44, the Moran’s I value of the 1990 Census median income. Direct estimates are generated with the sampling variances $D_1\mathrm{,}\dots \mathrm{,}{D}_{m_2}$ of the 1990 Current Population Survey (CPS) estimates. These sampling variances of sampled small areas range from 1.95 to 25.03 with the mean of 9.08, where dollar amounts are scaled by $1,000. For each setting, we consider $S=100$ replicated datasets.

Data generation: Let $\overline{D}=m_2^{-1}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>}1}^{m_2}}{D}_i$ and $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>}1}^{m_2}.$ We set $\rho =$ 0.85 and ${\sigma }_v^2=\overline{D}/2$ and consider two independent covariates $x_1$ and $x_2$ with SAR spatial dependence, i.e., $x_1\mathrm{,}{x}_2~N_m\left(0_m\mathrm{,}{\left\{{\Omega }_3\mathrm{<mo\; stretchy="false">(</mo>}\rho \mathrm{<mo\; stretchy="false">)</mo>}\right\}}^{-1}\right).$ Then, letting ${\mathrm{<mo\; stretchy="false">(</mo>}{\beta }_1\mathrm{,}{\beta }_2\mathrm{<mo\; stretchy="false">)</mo>}}^T={\mathrm{<mo\; stretchy="false">(</mo>2,}\mathrm{1<mo\; stretchy="false">)</mo>}}^T$ and $\mu ={\beta }_1{x}_1+{\beta }_2{x}_2,$ we generate small area means and direct estimates from the following independent FH model:

$$\theta ~N_m\mathrm{<mo\; stretchy="false">(</mo>}\mu \mathrm{,}{\sigma }_v^2{I}_m\mathrm{<mo\; stretchy="false">)</mo>,}{Y}_{\mathrm{<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}}|{\theta }_{\mathrm{<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}}~N_{m_2}\mathrm{<mo\; stretchy="false">(</mo>}{\theta }_{\mathrm{<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}}\mathrm{,}{D}_{\mathrm{<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>,}$$

where the components of $Y_{\mathrm{<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}}$ and ${\theta }_{\mathrm{<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}}$ correspond to the $m_2$ sampled small areas as defined below equation (2.13). The covariate $x_1\mathrm{<mo\; stretchy="false">(</mo>}{x}_2\mathrm{<mo\; stretchy="false">)</mo>}$ introduces stronger (weaker) spatial pattern to the ${\theta }_i\text{’s},$ and accordingly, we call $x_1$ $\mathrm{<mo\; stretchy="false">(</mo>}{x}_2\mathrm{<mo\; stretchy="false">)</mo>}$ as the strong (weak) covariate. Moran’s I values of 100 replicated small area means range from 0.115 to 0.713 with mean 0.449.

We consider two different covariate settings to examine how the spatial models can capture extra variability introduced by the spatial dependence from a missing covariate, i.e., $X=\mathrm{<mo\; stretchy="false">[</mo>}{1}_m\mathrm{,}{x}_2\mathrm{<mo\; stretchy="false">]</mo>}$ and $X=\mathrm{<mo\; stretchy="false">[</mo>}{1}_m\mathrm{,}{x}_1\mathrm{<mo\; stretchy="false">]</mo>},$ where $1_m$ represents the $m$ -component vector of ones. Excluding any of the covariates from the fitted model will leave the spatial variation of that covariate to the residual. We do not consider the full model involving both the covariates since that model will fully capture $\mu$ and leave no spatial variability unexplained; consequently, the independent FH model will be sufficient to capture the variability of the i.i.d. random effects.

Posterior simulations: For all models proposed, independent posterior samples can be obtained by the rejection sampling scheme outlined in Section 3. The sampling procedure begins with the rejection sampling from the marginal posterior distribution of $\mathrm{<mo\; stretchy="false">(</mo>}{\sigma }_v^2\mathrm{,}\rho \mathrm{<mo\; stretchy="false">)</mo>}$ and continues with successive samplings of the rest of the parameters from the conditional posterior distributions. However, when the marginal posterior density of $\mathrm{<mo\; stretchy="false">(</mo>}{\sigma }_v^2\mathrm{,}\rho \mathrm{<mo\; stretchy="false">)</mo>}$ is concentrated at or around the boundaries of $\rho ,$ a proposal distribution must be carefully specified to have a sufficiently high acceptance rate, which may require adaptive specification a proposal distribution for each replicated dataset. To avoid such difficulties, we use Hamiltonian Monte Carlo algorithm with the R package rstan (Stan Development Team, 2018). We fit the HB model (2.11)-(2.13) with each combination of covariates for $k=\mathrm{1,}\dots \mathrm{,}5.$ For each model, we run four parallel Hamiltonian Monte Carlo chains (No-U-Turn Sampler) for 2,500 iterations after 5,000 burn-in iterations. We keep every 10^th iteration and concatenate the four chains to obtain a posterior sample of size 10,000. The R codes implementing the rejection sampling step described in Section 3 and the stan models are available at https://github.com/heech31/spatial_sae.

Measures of performance: With the posterior sample under each model, we predict the true small area mean vector, ${\theta }^{\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}}={\mathrm{<mo\; stretchy="false">(</mo>}{\theta }_1^{\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}\dots \mathrm{,}{\theta }_m^{\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}}^T,$ of the $s^{\text{th}}$ replicated dataset using the posterior mean, which we denote by ${\widehat{\theta }}^{\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}}={\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{\theta }}_1^{\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}},\dots \mathrm{,}{\widehat{\theta }}_m^{\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}}^T.$ Let $A$ be a subset of $\mathrm{<mo>\{</mo>1,}\dots \mathrm{,}m\mathrm{<mo>\}</mo>},$ which is determined by indices only of the sampled or non-sampled small areas. For a given subset $A,$ we calculate the mean squared prediction error, ${\text{MSPE}}^{\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}}={\displaystyle {\sum }_{i\in A}}{\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{\theta }}_i^{\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}}-{\theta }_i^{\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}}^2/m_A,$ where $m_A=\left|A\right|$ is the number of areas in $A.$ We then average ${\text{MSPE}}^{\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo><mtext>\hspace{0.17em}</mtext>}}$ over $S$ replications to compute the average empirical mean squared prediction error (AeMSPE), where

$$\text{AeMSPE}=\frac{1}{S}{\displaystyle \sum _{s\mathrm{<mo>=</mo>1}}^S}{\text{MSPE}}^{\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}}=\frac{1}{S}{\displaystyle \sum _{s\mathrm{<mo>=</mo>}1}^S}\frac{1}{m_A}{\displaystyle \sum _{i\in A}}{\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{\theta }}_i^{\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}}-{\theta }_i^{\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}}^2\mathrm{.}(4.1)$$

We also evaluate the uncertainty of the predictions using the average posterior standard deviation (APSD) defined as $S^{-1}{\displaystyle {\sum }_{s\mathrm{<mo>=</mo>1}}^S}{m}_A^{-1}{\displaystyle {\sum }_{i\in A}}\text{sd}\mathrm{<mo\; stretchy="false">(</mo>}{\theta }_i^{\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>},$ where $\text{sd}\mathrm{<mo\; stretchy="false">(</mo>}{\theta }_i^{\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}$ is the posterior standard deviation of ${\theta }_i.$ By setting the independent FH model as a reference model, we consider the following ratios:

$$\text{AeMSPE}-{\text{Ratio}}_k=\frac{{\text{AeMSPE}}_k}{{\text{AeMSPE}}_1}\mathrm{,}\text{APSD}-{\text{Ratio}}_k=\frac{{\text{APSD}}_k}{{\text{APSD}}_1}\mathrm{,}(4.2)$$

where the subscript $k$ indicates that the quantity is calculated from the posterior sample under the $k^{\text{th}}$ model. These ratios measure improvements in AeMSPE and APSD achieved by fitting a spatial model over the independent FH model. A ratio less than 1 indicates superiority of the spatial model, otherwise the independent FH model is better. The smaller the ratio is, the more superior a spatial model is.

Model comparison: Various plots of Figure 4.1 summarize the ratios when the strong covariate $x_1$ is excluded from the fitted models. We categorize AeMSPE-Ratios and APSD-Ratios under each replication into three groups based on the Moran’s I values of ${\theta }_i\text{’s},$ namely the lowest, middle, and biggest thirds, respectively. The first row summarizes the prediction results for the seven non-sampled areas. As expected, spatial models show remarkable improvement in AeMSPE and APSD, and the improvements are greater when Moran’s I values are larger. In terms of AeMSPE, the SAR and LCAR models produce at least 20%, 30%, and 50% more accurate predictions when Moran’s I values are in the first, second, and third groups, respectively. In terms of the uncertainty of prediction, predictions of the SAR and LCAR models have 10% smaller APSD for the first and second groups. In the third group, the SAR model shows more than 25% reduction in APSD.

For the sampled areas with direct estimates, the improvements in AeMSPE are less than 10% for the first group but more than 15% and 25% for the second and third groups (grouped by Moran’s I values), respectively. Additionally, predictions from spatial models have a lower level of uncertainty, and for the third group, the SAR model predictions have more than 10% smaller APSD.

Similarly, various plots of Figure 4.2 summarize the ratios when the weak covariate $x_2$ is excluded from the fitted models. In general, the spatial models continue to produce better predictions over the independent FH model. The LCAR model performs the best overall for the seven non-sampled areas, resulting in a reduction of over 10% in AeMSPE for the first and second groups, respectively, and around 25% in AeMSPE for the third group. In terms of uncertainty, the spatial models show more than 5% but less than 10% smaller APSD. For the sampled small areas, the SAR and LCAR models show an AeMSPE reduction of around 5%-13%, but APSD reductions are less than 5%. Unlike the previous results with the weak covariate, the improvements in AeMSPE and APSD are comparable across three groups that are categorized by Moran’s I values. This is because the Moran’s I values of small area means are mostly determined by the strong covariate, and once the strong covariate is present in the model to explain the spatial variability, the spatial variability in the residuals do not vary markedly across the three groups. We regroup the ratios in terms of the Moran’s I values of the residuals obtained by regressing the strong covariate on the small area means, and summarize the ratios in Figure 4.3. Under this categorization, the LCAR model shows the best performance illustrating 5%, 10%, and 15% AeMSPE reductions for the first, second, and third groups, respectively. It can be also seen that the bigger the Moran’s I value is, the more improvement the spatial models achieve.

• Previous
• Article main page
• Next