Robust Bayesian small area estimation
Section 3. Application

• Previous
• Next

3.1  Real data analysis

The data set is selected by a 10% random sampling of households in each area from a test demographic census completed in one municipality in Brazil. The municipality consists of 38,740 households in 140 small areas in total, and the number of households per area in the population ranges from 57 to 588. Thus the area sample sizes in the data set range from 6 to 59. We are interested in estimating the 140 population means of the head of household’s income. The response variable $y_i$ denotes the average income of the heads of households in $i^{\text{th}}$ area.

This data set includes two centered auxiliary covariates which are the respective small area population means of the educational attainment of the head of households (ordinal scale of $0-5)$ and the average number of rooms in households $\mathrm{<mo\; stretchy="false">(</mo>\; 1}-11+\mathrm{<mo\; stretchy="false">)</mo>}$ . Lastly, the data set contains the respective sampling variances which are calculated in the usual way. Since only area level data were provided to us and the true area means are known, we can compare the 140 small area predictions with the true area means respectively. The analysis suggests that our model performs better than other models where random effects are based on the normal distribution. For comparison, we use three other models.

The first one is the Fay-Herriot model, referred to as FH, with known sampling variances.

$$\begin{array}{l}{y}_i|{\theta }_i\text{}\mathrm{,}{v}_i\stackrel{\text{ind}\text{.}}{{\displaystyle \sim }}N\left({\theta }_i\mathrm{,}{v}_i\right)\mathrm{,}{\theta }_i|\beta \mathrm{,}{\sigma }_{\delta }^2\stackrel{\text{ind}\text{.}}{{\displaystyle \sim }}N\left(x_i^T\beta \mathrm{,}{\sigma }_{\delta }^2\right)\hfill \\ \pi \left(\beta \right)\propto \mathrm{1,}\pi \left({\sigma }_{\delta }^2\right)\sim \text{IG}\left(a_0\mathrm{,}{b}_0\right)\mathrm{,}\hfill \end{array}$$

where $a_0$ and $b_0$ are chosen to be 0.0001 (a small constant) to reflect the vague knowledge of ${\sigma }_{\delta }^2\mathrm{.}$

The second model suggested by You and Chapman (2006), referred to as YC, is a hierarchical Bayesian model given by

$$\begin{array}{l}{y}_i|{\theta }_i\text{}\mathrm{,}{v}_i\stackrel{\text{ind}\text{.}}{{\displaystyle \sim }}N\left({\theta }_i\mathrm{,}{v}_i\right)\mathrm{,}{\theta }_i|\beta \mathrm{,}{\sigma }_{\delta }^2\stackrel{\text{ind}\text{.}}{{\displaystyle \sim }}N\left(x_i^T\beta \mathrm{,}{\sigma }_{\delta }^2\right)\hfill \\ s_i^2|v_i\stackrel{\text{ind}}{{\displaystyle \sim }}G\left(\frac{n_i-1}{2}\mathrm{,}\frac{1}{2v_i}\right)\hfill \\ \pi \left(v_i\right)\sim \text{IG}\left(a_i\mathrm{,}{b}_i\right)\mathrm{,}\pi \left(\beta \right)\propto \mathrm{1,}\pi \left({\sigma }_{\delta }^2\right)\sim \text{IG}\left(a_0\mathrm{,}{b}_0\right)\mathrm{,}\hfill \end{array}$$

where $a_i\mathrm{,}{b}_i\mathrm{,}{a}_0$ and $b_0$ are also chosen to be 0.0001.

The third model is a Bayesian multi-stage small area model proposed by Sugasawa et al. (2017), referred to as STK. The STK model produces shrinkage estimation of both means and variances.

$$\begin{array}{l}{y}_i|{\theta }_i\mathrm{,}{v}_i\stackrel{\text{ind}\text{.}}{{\displaystyle \sim }}N\left({\theta }_i\mathrm{,}{v}_i\right)\mathrm{,}{\theta }_i|\beta \mathrm{,}{\sigma }_{\delta }^2\stackrel{\text{ind}\text{.}}{{\displaystyle \sim }}N\left(x_i^T\beta \mathrm{,}{\sigma }_{\delta }^2\right)\hfill \\ s_i^2|v_i\stackrel{\text{ind}}{{\displaystyle \sim }}G\left(\frac{n_i-1}{2}\mathrm{,}\frac{1}{2v_i}\right)\mathrm{,}{v}_i|\gamma \sim \text{IG}\left(a_i\mathrm{,}{b}_i\gamma \right)\mathrm{,}\hfill \\ \pi \left(\beta \mathrm{,}{\sigma }_{\delta }^2\text{}\mathrm{,}\gamma \right)\mathrm{<mo>=</mo>\; 1,}\hfill \end{array}$$

where $a_i\mathrm{<mo>=</mo>\; 2}$ and $b_i\mathrm{<mo>=</mo>}1/n_i$ as suggested by authors for a reasonable choice.

We compare the small area means predicted by FH, YC, STK, and our model, hereafter referred to as RTS model. For the MCMC implementation, we generate a chain with a burn-in length of 50,000 and the sampling size of $G\mathrm{<mo>=</mo>}$50,000. The estimates of the ${\theta }_i$ are given by

$${\widehat{\theta }}_i\mathrm{<mo>=</mo>}\frac{1}{G}{\displaystyle \sum _{g\mathrm{<mo>=</mo>\; 1}}^G}\left({\gamma }_i^{\left(g\right)}{y}_i+\left(1-{\gamma }_i^{\left(g\right)}\right)x_i^T{\beta }^{\left(g\right)}\right)$$

where

$${\gamma }_i^{\left(g\right)}\mathrm{<mo>=</mo>}\frac{v_i^{-1\left(g\right)}}{v_i^{-1\left(g\right)}+{\sigma }_{\delta }^{-2\left(g\right)}{\eta }_i^{\left(g\right)}}\mathrm{.}$$

The comparison criteria are the average squared deviation (ASD), average absolute bias (AAB), average squared relative bias (ASRB), and average relative bias (ARB). They are defined as follows;

$$\text{ASD}=\frac{1}{m}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^m}{\left({\widehat{\theta }}_i-{\theta }_i\right)}^2\text{}\mathrm{,}\text{AAB}=\frac{1}{m}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^m}\left|{\widehat{\theta }}_i-{\theta }_i\right|\mathrm{,}\text{ASRB}=\frac{1}{m}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^m}{\left(\frac{{\widehat{\theta }}_i-{\theta }_i}{{\theta }_i}\right)}^2\mathrm{,}$$

and

$$\text{ARB}\mathrm{<mo>=</mo>}\frac{1}{m}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^m}\left|\frac{{\widehat{\theta }}_i-{\theta }_i}{{\theta }_i}\right|\mathrm{,}$$

where ${\widehat{\theta }}_i$ and ${\theta }_i$ are the estimated and true values respectively in the $i^{\text{th}}$ area. Table 3.1 compares the four models. Recall the prior distribution of ${\sigma }_{\delta }^2\mathrm{,}$ which is $\pi \left({\sigma }_{\delta }^2\right)\sim \text{IG}\left(\frac{p}{2}\mathrm{,}\frac{a}{2}\right)\mathrm{.}$ With the shape parameter, $\frac{p}{2}\mathrm{<mo>=</mo>\; 1,}$ we consider several values of $a\mathrm{.}$ If we choose $a$ to be close to $p\mathrm{,}$ RTS model fits better than the rest under all four criteria. When we choose $a\mathrm{<mo>=</mo>\; 1,}$ RTS model performs best. YC model performs worse than the other three models. If we choose very small $a,$ such as 0.01 or 0.001, then RTS model performs the worst.

Additionally, as suggested by a referee, we compute the residuals fitting a regression with the true area means and covariates to see the distribution of the random effects for this real data. Figure 3.1 shows that the distribution departs from the normal distribution.

3.2  Simulation study

In this section, we set up a simulation close to Maiti et al. (2014) (or Sugasawa et al. (2017)) to compare the accuracy of our estimators to other estimators, specifically those from You and Chapman (2006) and Sugasawa et al. (2017). We generate observations for each small area from the model

$$y_{ij}\mathrm{<mo>=</mo>}{\beta }_0+{\beta }_1{x}_i+u_i+e_{ij}\mathrm{,}j\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}{n}_i\text{}\mathrm{,}i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}m\mathrm{,}$$

where $u_i\sim t_{\nu }\left(\mathrm{0,}{\sigma }_{\delta }\right)$ and $e_{ij}\sim N\left(\mathrm{0,}{n}_i{v}_i\right)\mathrm{.}$ Then the random effects model for the small area mean is

$$y_i\mathrm{<mo>=</mo>}{\beta }_0+{\beta }_1{x}_i+u_i+e_i\mathrm{,}i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}m\mathrm{,}$$

where $y_i\mathrm{<mo>=</mo>}{\overline{y}}_i\mathrm{<mo>=</mo>}{n}_i^{-1}{\displaystyle {\sum }_{j\mathrm{<mo>=</mo>\; 1}}^{n_i}}{y}_{ij}$ and $e_i\mathrm{<mo>=</mo>}{n}_i^{-1}{\displaystyle {\sum }_{j\mathrm{<mo>=</mo>\; 1}}^{n_i}}{e}_{ij}\text{}\mathrm{.}$ Hence, $y_i|{\theta }_i\sim N\left({\theta }_i\mathrm{,}{v}_i\right)$ where ${\theta }_i\mathrm{<mo>=</mo>}{\beta }_0+{\beta }_1{x}_i+u_i;$ ${\theta }_i\sim t_{\nu }\left({\beta }_0+{\beta }_1{x}_i\text{}\mathrm{,}{\sigma }_{\delta }\right)\mathrm{,}$ and $e_i\sim N\left(\mathrm{0,}{v}_i\right)\mathrm{.}$ The interest parameter is the mean ${\theta }_i\mathrm{,}$ for the $i^{\text{th}}$ small area. Also, the direct estimator of $v_i$ is

$$\frac{1}{n_i\left(n_i-1\right)}{\displaystyle \sum _{j\mathrm{<mo>=</mo>\; 1}}^{n_i}}{\left(y_{ij}-{\overline{y}}_i\right)}^2\text{}\mathrm{.}$$

We set $m\mathrm{<mo>=</mo>\; 30}$ and $n_i\mathrm{<mo>=</mo>\; 7}$ for all areas, and generate covariates $x_i$ from the uniform distribution on (2, 8). The true parameter values are set as ${\beta }_0\mathrm{<mo>=</mo>}$0.5, ${\beta }_1\mathrm{<mo>=</mo>}$0.8, ${\sigma }_{\delta }\mathrm{<mo>=</mo>\; 1,}$ $\nu \mathrm{<mo>=</mo>\; 3}$ and $v_i\sim \text{IG}\left(\mathrm{10,5}\text{exp}\left(\text{0}\text{.3}{x}_i\right)\right)\mathrm{.}$ Also, we chose $a\mathrm{<mo>=</mo>\; 3}$ for all simulations.

For the MCMC implementation, we generated 5,000 posterior samples after discarding the first 1,000 for $R\mathrm{<mo>=</mo>}$2,000 simulation runs. Table 3.2 provides comparison among the four models. The comparison criteria are ASD, AAB, and BIAS, the latter being defined as

$$\text{BIAS}=\frac{1}{mR}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^m}\left|{\displaystyle \sum _{r\mathrm{<mo>=</mo>\; 1}}^R}\left({\widehat{\theta }}_i^{\left(r\right)}-{\theta }_i^{\left(r\right)}\right)\right|\mathrm{.}$$

RTS model performs better than others under ASD, AAB, and BIAS criteria for the mean. While RTS shows small improvements over other models for AAB and BIAS criteria, it shows approximate 23.5%, 10%, and 35.7% improvements over STK, FH and YC models respectively for ASD criteria. For the variance, RTS and STK models perform better than YC model.

The following two tables provide the simulation results when one sets the degrees of freedom as $\nu \mathrm{<mo>=</mo>\; 2}$ and $\nu \mathrm{<mo>=</mo>\; 4.}$

With $\nu \mathrm{<mo>=</mo>\; 2,}$ RTS model performs better than others under ASD, AAB, and BIAS criteria for the mean. In this simulation, the ASD values for STK and YC models are very large compared with RTS and FH models. RTS model shows improvements of about 78.6% over STK model, 82.2% over YC model, and 15.8% over FH model. For AAB and BIAS, the values of RTS model are smaller than those of other models. When considering the variance, ASD for RTS model gives smallest value.

With $\nu \mathrm{<mo>=</mo>\; 4,}$ RTS model also shows better performance over others. Especially, ASD and BIAS values indicate that RTS model improves results when compared with STK and YC model.

The next two tables consider the situation where one assumes normality of the random effects. Here RTS model performs slightly worse than the other models.

• Previous
• Article main page
• Next