Bayesian benchmarking of the Fay-Herriot model using random deletion
Section 5. Concluding remarks

• Previous
• Next

The Bayesian Fay-Herriot (BFH) model is discussed in detail. We show that the BFH can be fit using random samples rather than a Markov chain Monte Carlo sampler. Since random samples required no monitoring, this method is beneficial because there is little time at NASS between receiving the county-level survey summary data and presenting the final estimates. In support to the BFH model, we show that the posterior density under the BFH model is proper, providing a baseline for benchmarking. The effects of benchmarking are studied in a simulation study, comparing the BFH model without benchmarking to the BFH model with two benchmarking methods.

In this study, we assume that the benchmarking constraint is of the form ${\displaystyle {\sum }_{i=1}^l{\theta }_i=a.}$ A straightforward generalization of the benchmarking methods may be developed for the constraint of the form ${\displaystyle {\sum }_{i=1}^l{w}_i{\theta }_i=a,}$ where the $w_i$ are weights. For example, this latter situation occurs for benchmarking yield, ratio of production and harvested acres.

Our major contribution is the extension of BFH model to accommodate benchmarking. Previous approaches delete the last area, giving rise to the question “Does it matter which area is deleted?”. In this paper, we develop and illustrate a method that gives each area a chance to be deleted. We show how to fit this extended BFH model using the Gibbs sampler. Because of the complexity of the joint posterior density, a sampling based method, without Markov chains, cannot be used. Using empirical studies, we show that the differences in the posterior means over no benchmarking, deleting the last county and random deletion are very small.

The effects of changing the benchmarking target are studied in a sensitivity analysis. As expected, changing the benchmarking target leads to different estimates, but, unexpectedly, the changes in the posterior standard deviations are small. Small changes in the estimates are noted for the benchmarking methods using different probabilities of deletion.

It is expected that the posterior standard deviations from deleting the last one benchmarking and random benchmarking be larger than those from the BFH model because of the jittering effect from benchmarking. However, in the empirical studies we present, deleting the last one benchmarking and random benchmarking have about the same posterior standard deviations with a small reduction when random benchmarking is used. The key strength of the random benchmarking approach is that there is no preferential treatment for any area/county.

Disclaimer and acknowledgements

The Findings and Conclusions in This Preliminary Publication Have Not Been Formally Disseminated by the U.S. Department of Agriculture and Should Not Be Construed to Represent Any Agency Determination or Policy. This research was supported in part by the intramural research program of the U.S. Department of Agriculture, National Agriculture Statistics Service.

Dr. Nandram’s work was supported by a grant from the Simons Foundation (353953, Balgobin Nandram). The authors thank the Associate Editor and the referees for their comments and suggestions. The work of Erciulescu was completed as a Research Associate at the National Institute of Statistical Sciences (NISS) working on NASS projects.

Appendix A

Exemplification of the sensitivity of deletion

Let $y_i\stackrel{\text{ind}}{{\displaystyle \sim }}\text{Normal}\left({\mu }_i,{\sigma }_i^2\right),i=1,2,$ such that $y_1+y_2=a,$ and $\lambda ={\sigma }_2^2/\left({\sigma }_1^2+{\sigma }_2^2\right).$ Then, if we start by deleting $y_2,$ the joint density of $\left(y_1,y_2\right)$ is

$$f\left(y_1,y_2|\varphi =0\right)={\delta }_{y2}\left(a-y_1\right)\text{Normal}\left\{{\lambda }_{{\mu }_1}+\left(1-\lambda \right)\left(a-{\mu }_2\right),\left(1-\lambda \right){\sigma }_2^2\right\},$$

where ${\delta }_a\left(b\right)=1$ if $a=b$ and ${\delta }_a\left(b\right)=0$ if $a\ne b.$ However, if we start by deleting $y_1,$ the joint density of $\left(y_1,y_2\right)$ is

$$g\left(y_1,y_2|\varphi =0\right)={\delta }_{y1}\left(a-y_2\right)\text{Normal}\left\{\lambda \left(a-{\mu }_1\right)+\left(1-\lambda \right){\mu }_2,\left(1-\lambda \right){\sigma }_2^2\right\}.$$

It will matter in the estimation procedure which variable is deleted because the two joint distributions are different. Note that the two distributions are the same if and only if

$$\lambda {\mu }_1+\left(1-\lambda \right)\left(a-{\mu }_2\right)=\lambda \left(a-{\mu }_1\right)+\left(1-\lambda \right){\mu }_2,$$

which gives

$$\lambda \left({\mu }_1-a/2\right)=\left(1-\lambda \right)\left({\mu }_2-a/2\right).$$

Even if we assume that ${\sigma }_1^2={\sigma }_2^2,$ the two distributions are different. However, under this assumption, $\lambda =1/2,$ and the condition for the two distributions to be the same is that ${\mu }_1={\mu }_2.$ That is, overall the condition for the two joint distributions to be the same is that ${\mu }_1={\mu }_2$ and ${\sigma }_1={\sigma }_2,$ thereby making $y_1$ and $y_2$ exchangeable. However, this is a very restricted situation.

One way out of this diffculty is to actually delete both $y_1$ and $y_2$ in the following way. Let $z=1$ if $y_1$ is deleted and let $z=0$ if $y_2$ is deleted. Then,

$$p\left(y_1,y_2,z|\varphi =0\right)={\left[pg\left(y_1,y_2|\varphi =0\right)\right]}^z{\left[\left(1-p\right)f\left(y_1,y_2|\varphi =0\right)\right]}^{1-z},$$

where we have taken $z\sim \text{Bernoulli}\left(p\right)$ and, because $z$ is not really identifiable, we will take $p=1/2$ (i.e., we randomly delete one or the other). However, note that

$$z|y_1,y_2,\varphi =0\sim \text{Bernoulli}\left\{\frac{pg\left(y_1,y_2|\varphi =0\right)}{\left[pg\left(y_1,y_2|\varphi =0\right)\right]+\left[\left(1-p\right)f\left(y_1,y_2|\varphi =0\right)\right]}\right\}.$$

Appendix B

Fitting the Bayesian Fay-Herriot model

The Bayesian Fay-Herriot (BFH) model is given in (2.1) and the joint posterior density under the BFH model is given in (2.3), which for convenience we state here,

$$\pi \left(\theta ,\beta ,{\sigma }^2|\widehat{\theta }\right)\propto \frac{1}{{\left(1+{\sigma }^2\right)}^2}{\left(\frac{1}{{\sigma }^2}\right)}^{l/2}{\displaystyle \prod _{i=1}^l\left\{\exp \left[-\frac{1}{2}\left\{\frac{1}{s_i^2}{\left({\widehat{\theta }}_i-{\theta }_i\right)}^2+\frac{1}{{\sigma }^2}{\left({\theta }_i-x_i^{\prime }\beta \right)}^2\right\}\right]\right\}.}(\text{B}\text{.1})$$

We show how to fit the joint posterior density of the parameters using random samples (not even a Gibbs sampler), thereby avoiding any monitoring. We will use the multiplication rule to write

$$\pi \left(\theta \mathrm{,}\beta \mathrm{,}{\sigma }^2|\widehat{\theta }\right)\mathrm{<mo>=</mo>}{\pi }_1\left(\theta |\beta \mathrm{,}{\sigma }^2\mathrm{,}\widehat{\theta }\right){\pi }_2\left(\beta |{\sigma }^2\mathrm{,}\widehat{\theta }\right){\pi }_3\left({\sigma }^2|\widehat{\theta }\right)\mathrm{,}$$

where ${\pi }_1\left(\theta |\beta ,{\sigma }^2,\widehat{\theta }\right)$ and ${\pi }_2\left(\beta |{\sigma }^2,\widehat{\theta }\right)$ have standard forms and ${\pi }_3\left({\sigma }^2|\widehat{\theta }\right)$ is nonstandard but it is density of a single parameter.

Momentarily, we will drop the term, $\frac{1}{{\left(1+{\sigma }^2\right)}^2}{\left(\frac{1}{{\sigma }^2}\right)}^{l/2},$ because it only affects the posterior density of ${\sigma }^2.$ That is,

$${\pi }_1\left(\theta |\beta \mathrm{,}{\sigma }^2\mathrm{,}\widehat{\theta }\right)\propto {\displaystyle \prod _{i\mathrm{<mo>=</mo>\; 1}}^l}\left\{\exp \left[-\frac{1}{2}\left\{\frac{1}{s_i^2}{\left({\widehat{\theta }}_i-{\theta }_i\right)}^2+\frac{1}{{\sigma }^2}{\left({\theta }_i-x_i^{\prime }\beta \right)}^2\right\}\right]\right\}\mathrm{.}$$

Standard calculations reduce the argument (without $-1/2)$ of the exponential term to

$$\frac{1}{\left(1-{\lambda }_i\right){\sigma }^2}{\left\{{\theta }_i-\left({\lambda }_i{\widehat{\theta }}_i+\left(1-{\lambda }_i\right)x_i^{\prime }\beta \right)\right\}}^2+\frac{{\lambda }_i}{{\sigma }^2}{\left({\widehat{\theta }}_i-x_i^{\prime }\beta \right)}^2,{\lambda }_i=\frac{{\sigma }^2}{s_i^2+{\sigma }^2},i=1,\dots ,l.$$

Hence, for ${\pi }_1\left(\theta |\beta ,{\sigma }^2,\widehat{\theta }\right),$

$${\theta }_i|\beta \mathrm{,}{\sigma }^2\mathrm{,}\widehat{\theta }\stackrel{\text{ind}}{{\displaystyle \sim }}\text{Normal}\left\{{\lambda }_i{\widehat{\theta }}_i+\left(1-{\lambda }_i\right)x_i^{\prime }\beta \mathrm{,}\left(1-{\lambda }_i\right){\sigma }^2\right\}\mathrm{,}i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l\mathrm{.}(\text{B}\text{.2})$$

Momentarily, we will drop the term, ${\displaystyle {\prod }_{i\mathrm{<mo>=</mo>\; 1}}^l}{\left[\left(1-{\lambda }_i\right){\sigma }^2\right]}^{1/2}.$ Then, integrating out the ${\theta }_i,$ we get

$${\pi }_2\left(\beta |{\sigma }^2\mathrm{,}\widehat{\theta }\right)\propto \exp \left\{-\frac{1}{2}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^l}\frac{{\lambda }_i}{{\sigma }^2}{\left({\widehat{\theta }}_i-x_i^{\prime }\beta \right)}^2\right\}\mathrm{.}$$

Hence, the exponent (without $-1/2)$ can be written as,

$${\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^l}\frac{{\lambda }_i}{{\sigma }^2}{\left({\widehat{\theta }}_i-x_i^{\prime }\widehat{\beta }\right)}^2+{\left(\beta -\widehat{\beta }\right)}^{\prime }{\widehat{\Sigma }}^{-1}\left(\beta -\widehat{\beta }\right)\mathrm{,}$$

where

$$\widehat{\beta }\mathrm{<mo>=</mo>}\widehat{\Sigma }{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^l}\frac{{\widehat{\theta }}_i{x}_i}{s_i^2+{\sigma }^2}\text{and}{\widehat{\Sigma }}^{-1}\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^l}\frac{x_i{x}_i^{\prime }}{s_i^2+{\sigma }^2}\mathrm{.}$$

It is worth noting that $\widehat{\beta }$ and $\widehat{\Sigma }$ are well defined for all ${\sigma }^2$ provided that the design matrix, $X,$ where $X^{\prime }=\left(x_1,\dots ,x_l\right)$ is full rank. Then,

$${\pi }_2\left(\beta |{\sigma }^2,\widehat{\theta }\right)\propto \exp \left\{-\frac{1}{2}{\displaystyle \sum _{i=1}^l\frac{{\lambda }_i}{{\sigma }^2}{\left({\widehat{\theta }}_i-x_i^{\prime }\widehat{\beta }\right)}^2-\frac{1}{2}{\left(\beta -\widehat{\beta }\right)}^{\prime }{\widehat{\Sigma }}^{-1}\left(\beta -\widehat{\beta }\right)}\right\}.$$

That is,

$$\beta |{\sigma }^2,\widehat{\theta }\sim \text{Normal}\left(\widehat{\beta },\widehat{\Sigma }\right).(\text{B}.3)$$

Now, integrating out $\beta$ and incorporating the terms in ${\sigma }^2,$ which were dropped, we have

$${\pi }_3\left({\sigma }^2|\widehat{\theta }\right)\propto Q\left({\sigma }^2\right)\frac{1}{{\left(1+{\sigma }^2\right)}^2}\mathrm{,}(\text{B}\text{.4})$$

where

$$Q\left({\sigma }^2\right)\mathrm{<mo>=</mo>}{\left|\widehat{\Sigma }\right|}^{1/2}{\displaystyle \prod _{i\mathrm{<mo>=</mo>\; 1}}^l}\frac{1}{{\left(s_i^2+{\sigma }^2\right)}^{1/2}}\exp \left\{-\frac{1}{2}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^l}\frac{1}{s_i^2+{\sigma }^2}{\left({\widehat{\theta }}_i-x_i^{\prime }\widehat{\beta }\right)}^2\right\}\mathrm{.}$$

To obtain a random sample from (B.1), we sample ${\sigma }^2$ from (B.4), $\beta$ from (B.3) and the ${\theta }_i$ independently from (B.2). The conditional posterior density in (B.4) is nonstandard, and to draw a sample from it, we use a grid method (e.g., Nandram and Yin, 2016). First, we transform ${\sigma }^2$ to $\varphi ={\sigma }^2/\left(1+{\sigma }^2\right)$ so that $0<\varphi <1.$ Then, we divide (0, 1) into 100 grids. Actually, we have located the range of $\varphi$ in (0, 1) and we have divided this interval into 100 grids. This gives us a probability mass function that we sample. Jittering is used in the selected grid to get deviates, which are different with probability one; see Nandram and Yin (2016) for more details.

Appendix C

Proof of Theorem 2

It is convenient to make the following transformations,

$${\theta }_i={\theta }_i,i=1,\dots ,l-1,\varphi ={\displaystyle \sum _{i=1}^l{\theta }_i-a}.$$

Here, $\varphi$ is a dummy variable, which holds the benchmarking constraint, and it ensures a non-singular transformation. The Jacobian is unity and the inverse transformation is

$${\theta }_i={\theta }_i,i=1,\dots ,l-1,{\theta }_l=\varphi +a-{\displaystyle \sum _{i=1}^{l-1}{\theta }_i}.$$

The transformed density is

$$\tilde{\pi }\left({\theta }_1,\dots ,{\theta }_{l-1},\varphi \right).$$

Then, the density that holds the benchmarking constraint exactly is

$$\tilde{\pi }\left({\theta }_1,\dots ,{\theta }_{l-1}|\varphi =0\right),{\theta }_l=a-{\displaystyle \sum _{i=1}^{l-1}{\theta }_i.}$$

Therefore,

$$\pi \left({\theta }_1,\dots ,{\theta }_{l-1}|\varphi =0\right)\propto \exp \left\{-\frac{1}{2{\delta }^2}\left[{\displaystyle \sum _{i=1}^{l-1}{\left({\theta }_i-u_i^{\prime }\beta \right)}^2+}{\left\{{\displaystyle \sum _{i=1}^{l-1}{\theta }_i-\left(a-u_i^{\prime }\beta \right)}\right\}}^2\right]\right\}.$$

Dropping terms that do not involve ${\theta }_{\left(l\right)}={\left({\theta }_1,\dots ,{\theta }_{l-1}\right)}^{\prime },$ it is easy to show that the exponent is

$$\frac{1}{2{\delta }^2}\left\{{\theta }_{\left(l\right)}^{\prime }\left(I+J\right){\theta }_{\left(l\right)}-2\left[{\left(u_1-u_l\right)}^{\prime }\beta ,\dots ,{\left(u_{l-1}-u_l\right)}^{\prime }\beta +a{j}^{\prime }\right]{\theta }_{\left(l\right)}\right\}.$$

Then, using the properties of a multivariate normal density, we have

$${\theta }_{\left(l\right)}|\varphi =0\sim \text{Normal}\left({\left(I+J\right)}^{-1}{\left(a{j}^{\prime }+{\left(u_1-u_l\right)}^{\prime }\beta ,\dots ,{\left(u_{l-1}-u_l\right)}^{\prime }\beta \right)}^{\prime },{\delta }^2{\left(I+J\right)}^{-1}\right).$$

Finally, using the Sherman-Morrison formula, ${\left(I+J\right)}^{-1}=I-\frac{1}{l}J,$ we have

$${\theta }_{\left(l\right)}|\varphi =0\sim \text{Normal}\left\{\left(I-\frac{1}{l}J\right){\left(a{j}^{\prime }+{\left(u_1-u_l\right)}^{\prime }\beta ,\dots ,{\left(u_{l-1}-u_l\right)}^{\prime }\beta \right)}^{\prime },{\delta }^2\left(I-\frac{1}{l}J\right)\right\}.$$

It is worth noting that the matrix determinant lemma gives $\det \left(I+J\right)=l$ and so $\det \left({\delta }^2\left(I-\frac{1}{l}J\right)\right)=\frac{1}{l}{\left({\delta }^2\right)}^{l-1}.$

References

Battese, G., Harter, R. and Fuller, W. (1988). An error-components model for prediction of county crop areas using survey and satellite data. Journal of the American Statistical Association, 83(401), 28-36. doi:10.2307/2288915.

Bell, W.R., Datta, G.S. and Ghosh, M. (2013). Benchmarking small area estimators. Biometrika, 100(1), 189-202.

Cruze, N.B., Erciulescu, A.L., Nandram, B., Barboza, W. and Young, L.J. (2019). Producing official county-level agricultural estimates in the United States: Needs and challenges. Statistical Science. To appear.

Datta, G.S., Ghosh, M., Steorts, R. and Maples, J. (2011). Bayesian benchmarking with applications to small area estimation. Test, 20(3), 574-588.

Erciulescu, A.L., Cruze, N.B. and Nandram, B. (2019). Model-based county level crop estimates incorporating auxiliary sources of information. Journal of Royal Statistical Society, Series A, 182, 283-303. doi:10.1111/rssa.12390.

Fay, R.E., and Herriot, R.A. (1979). Estimates of income for small places: An application of James-Stein procedures to census data. Journal of the American Statistical Association, 74(366), 269-277.

Ghosh, M., and Steorts, R. (2013). Two stage Bayesian benchmarking as applied to small area estimation. Test, 22(4), 670-687.

Janicki, R., and Vesper, A. (2017). Benchmarking techniques for reconciling Bayesian small area models at distinct geographic levels. Statistical Methods and Applications, 26, 4, 557-581.

Jiang, J., and Lahiri, P. (2006). Mixed model prediction and small area estimation. Test, 15(1), 1-96.

Nandram, B., and Sayit, H. (2011). A Bayesian analysis of small area probabilities under a constraint. Survey Methodology, 37, 2, 137-152. Paper available at https://www150.statcan.gc.ca/n1/en/pub/12-001-x/2011002/article/11603-eng.pdf.

Nandram, B., and Toto, M.C.S. (2010). Bayesian predictive inference for benchmarking crop production for Iowa counties. Journal of the Indian Society of Agricultural Statistics, 64 (2), 191-207.

Nandram, B., and Yin, J. (2016). A nonparametric Bayesian prediction interval for a finite population mean. Journal of Statistical Computation and Simulation, 86 (16), 3141-3157.

Nandram, B., Berg, E. and Barboza, W. (2014). A hierarchical Bayesian model for forecasting state-level corn yield. Journal of Environmental and Ecological Statistics, 21, 507-530.

Nandram, B., Toto, M.C.S. and Choi, J.W. (2011). A Bayesian benchmarking of the Scott-Smith model for small areas. Journal of Statistical Computation and Simulation, 81 (11), 1593-1608.

Pfeffermann, D. (2013). New important developments in small area estimation. Statistical Science, 28(1), 40-68.

Pfeffermann, D., Sikov, A. and Tiller, R. (2014). Single-and two-stage cross-sectional and time series benchmarking procedures for small area estimation. Test, 23(4), 631-666.

Pfeffermann, D., and Tiller, R. (2006). Small area estimation with state-space models subject to benchmark constraints. Journal of the American Statistical Association, 101 (476), 1387-1397.

Rao, J.N.K., and Molina, I. (2015). Small Area Estimation, Wiley Series in Survey Methodology.

Toto, M.C.S., and Nandram, B. (2010). A Bayesian predictive inference for small area means incorporating covariates and sampling weights. Journal of Statistical Planning and Inference, 140, 2963-2979.

Wang, J., Fuller, W.A. and Qu, Y. (2008). Small area estimation under a restriction. Survey Methodology, 34, 1, 29-36. Paper available at https://www150.statcan.gc.ca/n1/en/pub/12-001-x/2008001/article/10619-eng.pdf.

You, Y., and Rao, J.N.K. (2002). Small area estimation using unmatched sampling and linking models. The Canadian Journal of Statistics, 30, 3-15.

You, Y., Rao, J.N.K. and Dick, J.P. (2004). Benchmarking hierarchical Bayes small area estimators in the Canadian census under coverage estimation. Statistics in Transition, 6, 631-640.

• Previous
• Article main page
• Next