Bayesian benchmarking of the Fay-Herriot model using random deletion
Section 4. Empirical studies

• Previous
• Next

The purposes of these empirical studies are twofold. First, it is demonstrated that the BFH model can be fit as stated in Section 2 and the deleting the last one benchmarking and random benchmarking methods are performed. Second, the benchmarking methods are compared in a simulation study that uses a well-used dataset in the small area literature.

In the data generation process, we use the data on corn and soybean acres in Battese, Harter and Fuller (1988), available for 12 counties (areas) in Iowa. The resulting county-level corn and soybean acreages are constructed using a number of segments sampled from the population (known number of segments). Landsat satellite data on the number of pixels of corn and soybean in the sampled segments (i.e., two covariates) are also available. The finite population means of the number of pixels classified as corn and soybean for each county are also reported. Starting with this dataset, we construct new datasets with any number of areas.

The data generation process has two steps. In the first step, the unit-level model $y_{ij}=x_{ij}^{\prime }\beta +e_{ij},$ $i=1,\dots ,l,j=1,\dots ,n_i,$ where $e_{ij}\stackrel{\text{iid}}{{\displaystyle \sim }}\left(0,{\sigma }^2\right),$ is fit to the data available for the $l=12$ counties in Iowa. The area sample sizes are $n_1=n_2=n_3=1,$ $n_4=2,n_5=n_6=n_7=n_8=3,$ $n_9=4,n_{10}=n_{11}=5,$ and $n_{12}=6.$ Using least squares, we estimate $\beta$ and ${\sigma }^2$ by $\widehat{\beta }$ and ${\widehat{\sigma }}^2,$ respectively. For the areas with sample size greater than one, we set $s_i^2$ equal to the estimated variance of the sample mean ${\overline{y}}_i\left({\overline{y}}_i={\displaystyle {\sum }_{j=1}^{n_i}{y}_{ij}/n_i}\right)$ and we let $S^2$ be their geometric mean. For the areas with sample size equal to one, we set $s_i^2$ equal to $S^2.$ The vector of covariates ${\overline{X}}_i$ has three elements, the integer one (for the intercept), followed by the population means of pixels classified as corn and soybean.

In the second step, the data generation process for any desired number $l$ of small areas is illustrated. The covariates $x_i,i=1,\dots ,l,$ are sampled with replacement from ${\overline{X}}_i,i=1,\dots ,12.$ Then, the area-level means are drawn using

$${\theta }_i\stackrel{\text{ind}}{{\displaystyle \sim }}\text{Normal}\left(x_i^{\prime }\widehat{\beta },{\widehat{\sigma }}^2\right),i=1,\dots ,l,$$

where $\widehat{\beta }$ and ${\widehat{\sigma }}^2$ are the least squares estimates defined above. The sample variances $s_i^2$ are generated in two steps. First, the sample sizes are drawn from a uniform distribution, $n_i\stackrel{\text{iid}}{{\displaystyle \sim }}\text{Uniform}\left(\mathrm{5,}25\right)\mathrm{,}$ $i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l.$ Second, let $s_i^2\mathrm{<mo>=</mo>}{S}^2{V}_i/\left(n_i-1\right),$ where $V_i\stackrel{\text{ind}}{{\displaystyle \sim }}{\chi }_{n_i-1}^2$ and $S^2$ defined above. Finally, the small area survey estimates are drawn using ${\widehat{\theta }}_i\stackrel{\text{ind}}{{\displaystyle \sim }}\text{Normal}\left({\theta }_i\mathrm{,}{s}_i^2\right),i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l.$ The benchmarking target is set equal to the sum of the ${\widehat{\theta }}_i$ and variants of this value, ${\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^l{\widehat{\theta }}_i}$ scaled up or down by 50%. In NASS’s practice, for crop county estimates, this target is an already set state value. To evaluate the benchmarking methods in extreme cases, we consider additional simulation scenarios, where an area sample size is set to 2 or 50, or where the factor $S^2$ is multiplied by ten.

In what follows, we report empirical results mostly for a simulation scenario using 12 areas. Examples using larger number of areas are briefly discussed. For example, Iowa has 99 counties, and one of NASS’s interests is in benchmarking county estimates for planted acres, harvested acres and production (bushels) to the predefined state-level total. For such small numbers of areas, no adjustment is needed to the benchmarking procedures, deleting the last one or random deletion, introduced in the previous sections. However, the computation may be intolerable for an extremely large number of areas (say, one million), and some adjustments would be needed to the current procedures.

It is pertinent to discuss the computations for the simulation scenario with 12 areas. For posterior inference under the BFH model, we have used 1,000 random draws, and this runs in just a few seconds. On the other hand, it is more difficult to run a Gibbs sampler for deleting one at a time or random deletion benchmarking. However, we have provided an efficient Gibbs sampler as follows. We used a long run of 20,000 iterations, with a “burn in” of the first 10,000 iterations, choosing every tenth iterate thereafter. This was obtained by trial and error that is gauged by the autocorrelations, the Geweke test for stationarity and the effective sample sizes. For the 1,000 selected iterations, the autocorrelations are all negligible. For random deletion benchmarking, the p-values of the Geweke test for the three regression coefficients and ${\delta }^2$ are, respectively, 0.651, 0.087, 0.828 and 0.699 (i.e., stationarity is not rejected), and the effective sample sizes are all 1,000. Also, the trace plots show no evidence of nonstationarity. Therefore, the Gibbs sampler is efficient, taking a few seconds despite the large number of runs.

The performance of benchmarking methods is assessed using a set of metrics that include posterior means (PM) and posterior standard deviations (PSD), and when it is convenient, posterior coefficients of variation (PCV), numerical standard errors (NSE) of the estimates and 95% highest posterior density intervals (95% HPD). Numerical results are presented in Tables 4.1-4.8.

A summarized version of the basic results is presented in Table 4.1, and serves for comparison of the average, standard error and coefficient of variation of the observed data with the PMs, PSDs, PCVs from the BFH model, benchmarking (deleting the last one, LO) model and random benchmarking (RD) model. The results in Table 4.1 apply to two simulation scenarios, where $S^2\mathrm{<mo>=</mo>\; 163},$ small variation in the observed data, and where $S^2\mathrm{<mo>=</mo>}\text{1,630},$ relatively larger variation in the observed data. When $S^2\mathrm{<mo>=</mo>\; 163},$ there are very little differences between the observed data and the posterior quantities from the BFH, LO and RD models. Given the small coefficients of variation for the survey estimates, it is difficult for any model to further reduce variability. Hence, the PCVs are comparable to the CVs of the survey estimates. On the other hand, three interesting points can be made for the scenario where $S^2\mathrm{<mo>=</mo>}\text{1,630}.$ First, the PMs under the BFH model can be very different from those of LO and RD models and these latter two PMs are very close. Second, the PSDs are much smaller than the standard errors of the observed data; there are substantial gains in precision under the BFH model. However, the PSDs are about four to five times smaller than those for the observed data and the PSDs under the LO and RD model are about twice those of the BFH model. The PCVs follow the same pattern. Third, LO and RD are very close in all three measures (PMs, PSDs, PCVs) with RD model having just slightly smaller PSDs. As expected, there is small difference between the LO model and the RD model. But one must also observe that benchmarking the BFH model is important because we can get answers that are different from the BFH model at least in terms of posterior standard deviations and coefficients of variation. Benchmarking is a jittering procedure, which helps to protect the model from misspecification, and therefore it must lead to increased variability in the small area estimates.

Under the basic simulation scenario, we compare the deletion benchmarking methods to one of the methods in DGSM that provides benchmarked posterior estimates without deletion. To match the notation in DGSM, the benchmarking equation must be rewritten as

$${\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^l{\omega }_i{\theta }_i}\mathrm{<mo>=</mo>}\frac{a}{l}\mathrm{<mo>=</mo>}t\mathrm{,}$$

where ${\omega }_i\mathrm{<mo>=</mo>}1/l\mathrm{,}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^l{\omega }_i}\mathrm{<mo>=</mo>\; 1.}$ Let ${\widehat{\theta }}_i{}^{\left(B\right)}$ denote the posterior means from the BFH model. Now, define ${\overline{\widehat{\theta }}}_B\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^l{\omega }_i}{\widehat{\theta }}_i^{\left(B\right)},$

$${\varphi }_i\mathrm{<mo>=</mo>}\frac{{\omega }_i}{{\widehat{\theta }}_i^{\left(B\right)}}\mathrm{,}{r}_i\mathrm{<mo>=</mo>}\frac{{\omega }_i}{{\varphi }_i}\mathrm{,}i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l\mathrm{,}$$

and $S^{*}\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^l}{\omega }_i^2/{\varphi }_i\text{}.$ Note that among the several specifications in DGSM, we have selected ${\varphi }_i$ at random (no preference). Then, the benchmarked Bayes estimators of DGSM are

$${\widehat{\theta }}_i^{\left(BM\right)}\mathrm{<mo>=</mo>}{\widehat{\theta }}_i^{\left(B\right)}+\left(t-{\overline{\widehat{\theta }}}_B\right)r_i/S^{*}\text{},i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l\mathrm{.}$$

Empirical results using the estimator ${\widehat{\theta }}_i^{\left(BM\right)}$ are presented in the note to Table 4.1. The largest difference between the benchmarked estimates under different benchmarking methods is for area 10 ( OB: 172.9; BFH: 129.2; LO: 145.5; RD: 145.3; DGSM: 126.4). In general, the PMs from LO and RD are closer to OB (observed data). Otherwise, these estimates compare reasonably well with the LO benchmarking and RD deletion although there are some small differences; DGSM does not provide posterior standard deviations and credible intervals.

More detailed results for $S^2\mathrm{<mo>=</mo>\; 163},$ are presented in Tables 4.2-4.8 and in Figures 4.1-4.4. Our interest is mainly to compare deletion of a single area (e.g., LO) and RD.

Using the results in Table 4.2, we conclude that the PMs from the BFH model (without benchmarking) are slightly different from the direct estimates, and as expected, larger than the smaller direct estimates and smaller than the larger ones. Except for two areas, as expected, the PSDs are smaller than the direct standard deviations. For example, the smallest direct estimate (76.997) has the largest shrinkage with a larger standard deviation (5.881 vs. 5.995); the results are consistent with the standard shrinkage that occurs in small area estimation. We note that the PCVs are all small and the NSEs are reasonably small, too.

The estimates from the BFH model with deleting the last area and with random deletion under a uniform prior (equal weights) are presented in Tables 4.3 and 4.4. The posterior weights barely differ from 0.083 with the largest one (0.097) of the last area and smallest one (0.056) of the $8^{\text{th}}$ area. Both random deletion and deleting the last one provide improved precision, as the PSDs of the benchmarked estimates are all smaller than the observed standard errors, for both benchmarking methods. The NSEs are larger than for no benchmarking, but this barely matters as these are errors of the PMs (the characteristic of the PM has three digits).

The three methods (BFH, RD, LO) are compared using the results in Table 4.5. The PMs are comparable, so that benchmarking (RD, LO) does not distort (shrink) the estimates much beyond the shrinkage under the BFH model. Also, the PSDs under LO and RD are almost always smaller than those under the BFH model. For eight of the twelve areas, RD has smaller PSDs than LO; in these areas, RD shows roughly 1% decrease in PSD over LO and roughly 4% over the PSDs from BFH.

To investigate how sensitive the PSDs are to different benchmarking targets, we present results using three choices of targets in Table 4.6. The PSDs change only slightly over different targets and are still better than the standard errors of the direct estimates.

As part of designing a complex set of simulations, we consider using unequal probabilities (weights) in the random deletion benchmarking, and present results in Table 4.7. Uniform weights (EW) are compared to weights inversely proportional (IW) to the sample sizes and to weights directly proportional (DW) to the samples sizes. Again, small differences are present among the three PMs and among the three PSDs. The PSDs are still smaller than those of the direct estimates.

Using the results in Table 4.8, we study how extreme sample sizes in the last county (to be deleted) affect posterior inference. For this, we set the sample size of the last county to be outside the simulation range (5-25), at 2 and 50. First, consider the case in which the sample size of the last county is 2. Consistent with previous findings, there are minor differences of the PMs over no benchmarking, deleting the last one and random deletion for all counties. The PSDs for LO and RD are smaller than those of BFH with nine of these PSDs for RD smaller than LO. However, for the last county, we observe relatively large posterior standard deviations (10.00, 8.771, 8.525), roughly 15% decrease in PSD of RD over no benchmarking. Next, consider the case in which the sample size of the last county is 50. The patterns are similar, except the PSDs for the last county are comparable to the others under BFH, LO and RD and again there is an approximately 10% decrease (6.282, 5.958, 5.702) in PSD of RD over no benchmarking. It appears that deliberately putting the county with the most extreme sample size (small or large) as the last county can affect the benchmarking procedure. In contrast, minor changes are observed when the areas with extreme sample size are not systematically deleted. When the sample size is 2, the new PMs and PSDs are the following, BFH: 124.307, 9.993; LO: 123.371, 9.000 RD: 123.540, 8.887. When the sample size is 50, the new PMs and PSDs are the following, BFH: 118.167, 6.284; LO: 117.802, 6.094; RD: 117.716, 5.948.

For comparison, different posterior densities are presented in Figures 4.1-4.4. In Figures 4.1 and 4.2, we present posterior densities of all twelve area parameters when each area, in turn, is deleted. We observe that the posterior densities are slightly different around the modes, but nothing remarkable. In Figures 4.3 and 4.4, we present posterior densities of all twelve area parameters under the FH model (unconstrained), random deletion benchmarking and deleting the last one. There are some differences among the three densities, but again these are not alarmingly different.

Finally, empirical results are presented for a simulation scenario with 99 areas, reflecting the 99 counties in Iowa. The data are generated as previously described, and the BFH model without benchmarking, with random deletion benchmarking, and with deleting the last one benchmarking is fit using 20,000 iterations for the Gibbs sampler. For each model fit, the first 10,000 iterations are used as a burn-in and every tenth iteration is kept thereafter. The BFH model fitting takes 15 seconds, while the deletion benchmarking models takes slightly less than three minutes each. For the random deletion benchmarking model parameters, the regression coefficients $\beta$ and the variance ${\sigma }^2,$ the p-values of the Geweke test are, respectively, 0.822, 0.128, 0.752 and 0.219, and the effective sample sizes are all 1,000 for the 1,000 selected iterations (i.e., an efficient Gibbs sampler). Note that the target is 12,162.93 and the sum of the PMs from the BFH model is 12,168.49, a difference of 5.56. In Figure 4.5, we present a plot of the coefficients of variation under random deletion benchmarking, deleting the last one benchmarking and BFH model versus the direct estimates by area. The differences among these models are not remarkable. Most of the points with direct CVs larger than about 0.04 fall below the ${45}^{\text{o}}$ straight line. However, some points (diamond) under the BFH model are above the ${45}^{\text{o}}$ line, four of them are noticeable, possibly shrinking too much. We conclude that it is sensible to perform the random deletion benchmarking.

• Previous
• Article main page
• Next