Hierarchical Bayes Small Area Estimation under a Unit Level Model with Applications in Agriculture

To studied Bayesian aspect of small area estimation using Unit level model. In this paper we proposed and evaluated new prior distribution for the ratio ) ( / 2 2     e v of variance components in unit level model rather than uniform prior. To approximate the posterior moments of small area means, Laplace approximation method is applied. This choice of prior avoids the extreme skewness, usually present in the posterior distribution of variance components. This property leads to more accurate Laplace approximation. We apply the proposed model to the analysis of horticultural data and results from the model are compared with frequestist approach and with Bayesian model of uniform prior in terms of average relative bias, average squared relative bias and average absolute bias. The numerical results obtained highlighted the superiority of using the proposed prior over the uniform prior. Thus Bayes estimators (with new prior) of small area means have good frequentist properties such as MSE and ARB as compared to other traditional methods viz., Direct, Synthetic and Composite estimators.


Introduction
Model-based small area estimation methods have been widely used in practice due to the increasing demand for precise estimates for local regions and various small areas.It is now generally accepted that the indirect estimates should be based on explicit models that provide links to related areas through the use of supplementary data such as census counts or administrative records; see, for example, (Jiang and Lahiri 2006) and (Rao 2003) for more discussion on model-based small area methods.Also, (Adam et. Al., 2013) summarise the main methodological approaches to SAE and their linkages.[Jiango et. Al., 2013] investigate two new approaches: one relying on the work of Chambers, and the second using the concept of conditional bias to measure the influence of units in the population.[Chambers et. al., 2014] proposed two different analytical mean-squared error estimators for the ensuing bias-corrected outlier robust estimators.[Rao et al., 2013] relaxed the assumption of linear regression for the fixed part of the model and replace it by a weaker assumption of a semi-parametric regression.The model-based estimates are obtained to improve the direct design-based estimates in terms of precision and reliability, i.e., smaller coefficients of variation (CVs).There are two broad classifications for small area models: area level models and unit level models.The basic unit level model is based on unit level auxiliary variables.In this paper we focus on Unit level models that are related to the unit level values of response through a nested error linear regression model, under the assumption that the nested error and the model error are independent of each other and normally distributed with common mean zero and common or different variances.The nested error unit level regression model was first used to model county crop areas in USA (Battese et al., 1988), they have used the normally distributed common errors variance assumption and revealed that based on the fitting-of-constants method the estimates of errors variances are slightly different from each other.Techniques for validating their model on the basis of unit level auxiliary variables are also considered.This type of model is appropriate for continuous value response variables.Various extensions of this type of model have been proposed to handle binary responses, two-stage sampling within areas, multivariate responses and others (Rao, 2003(Rao, , 2003a)).
The objective of this paper is to consider new improved prior on hyperparameters of variance component and obtain HB inference for small area parameters through the MCMC method using Laplace approximation.In section 4, we apply the proposed model to the analysis of small area data from the Horticultural Survey.We compare the performance the proposed model i.e hierarchical Bayes estimate with the proposed prior (Proposed HB) with the hierarchical Bayes estimate with uniform prior (HB(WINBUG)) and EBLUP estimates to investigate the effects of incorporating new prior on the variance component .Bayesian model comparison and model fit analysis are also provided.Finally in section 5, we offer some concluding remarks.

Unit level model
The basic unit level model is based on unit level auxiliary variables.This type of model can be represented by the following mathematical equation where ij y is the value of the study variable for the th j unit belonging to the th i area, ij x is the unit level vector of unknown regression coefficients, i v is the random small area effect, and ij e is the error term.We

EBLUP estimator of unit level model
For the basic unit level model defined in section (2) the BLUP estimator of i  is given as: The BLUP estimator (2.1.1)depends on the variance ratio where The second order MSE approximation i.e EBLUP estimate of  Rao (2003), we assume that the sample values also follow the model (3.4), that is, we do not address the situation involving informative sampling.The assumption that the sample values also obey the model (3.5.1) holds true for simple random sampling from each area.For a proof of this absence of selection bias (Rao, 2003).
We can also write the small area mean where is the finite population correction, is y and ins y are the means of the sampled and non-sampled units, respectively, for the th i area.Now applying a hierarchical Bayesian approach to estimate the small area means (3.1) for finite populations.A hierarchical Bayesian version of the Unit level model is given by Our aim is to find the Bayesian estimator of small area means of unit level model and its measure of uncertainty, which are given by ) respectively, where y is the vector of sampled values of y.Whether we are in Bayesian or frequentist paradigm depends on whether we assume a prior on the hyperparameters ) , , ( at the third level.To obtain the Bayesian summary statistics for Yi, we proceed as follows (assuming that the sample values follow the same model above, as discussed in Section 2).First we write the likelihood using level 1 and level 2 of the hierarchical Bayesian model as: Now, following Ghosh and Lahiri (1989), Datta and Ghosh (1991), we specify prior on 2 e  and  at level 3 as: 3), we consider a flat prior (uniform prior) on 2 e  i.e., we consider Typically, enough data will be available to estimate 2 e  precisely.Any reasonable non- informative prior usually works well for 2 e  (Gelman, 2006).
After collecting terms involving i v from the exponent of (3.7.3), we can say that: 3), using matrix notation, we write the joint distribution (without the normalizing constant) of where The subscript  in both   ˆand   indicates the dependence of the terms on the ratio of variance component From the equation (3.7), it follows that where IG(a,b) is the inverse gamma distribution with shape parameter "a" and scale parameter "b".A random variable Z has an inverse Gamma distribution IG(a,b), if its pdf is given by 0 Using the property of inverse gamma distribution, we express the conditional posterior mean and variance of 2 Integrating out 2 e  from (3.7.7), we can write the posterior distribution of a single variance component  as To present the results in a unified way, we restate (3.11) using matrix notation as where This formulation is presented for the ease of writing codes in simulation and data analysis.

Posterior moments of finite population mean when  known
We derive in two steps.First, we find analytical expressions for these quantities when  is known.Then using the Laplace approximation, we obtain the final results.Using the iterative expectation and variance technique and the results given in (3..4), (3. 6), and (3.9), it is not very difficult to get Common mean model, Balanced case, and 0  i f Here, besides the common mean assumption, we assume that . Then the posterior distribution of  is given by: are the usual definition of within sum of square and between sum of square used in the balanced one-way ANOVA.After some modification and simplification of (3.7.13) and (3.7.14), we find

2 Choice of proir on 
To obtain the posterior mean and variance of i Y from (3.1.1)and (3.1.2),we need to perform one-dimensional integral with respect to the posterior distribution of  .For that we need to assume a prior distribution on  .The uniform prior on  is non-informative and yields a posterior distribution of  for which the mode is identical to the residual maximum likelihood (REML) estimator of  , following the arguments discussed above.The posterior mode of  with uniform prior or equivalently, the REML estimator of  can be zero for a particular application.In many practical applications, the maximum likelihood (ML) or restricted maximum likelihood (REML) estimates of hyperparameters occur at the boundary point.we match the posterior distribution of  given in (3. 11) to the adjusted profile likelihood function of  to obtain an appropriate prior distribution of  .This prior leads to a posterior distribution of  for which the mode is always positive and results in estimators of the shrinkage factor and small area mean that have good frequentist properties.By profile likelihood, we mean the likelihood of  that does not account for the loss of degrees of freedom due to the estimation of regression coefficient  .This is given by Note that, to obtain the profile likelihood (3.2.1), we consider  as the only nuisance parameter.We plug in the ML estimate of  in (3.5) [without the term


] to obtain the joint profile likelihood of  and 2 e  .Then integrating out 2 e  , we derive the marginal profile likelihood of  .This leads to the prior 0 To simplify the prior for some special cases, we proceed as follows: . Thus, we can write For the common mean model, balanced case (i.e.
where SSW and SSB is defined earlier.Hence, simplifying it further we can say

3 Laplace approximation
The posterior moments of i Y , using the prior in (3.2.2) are not in a closed-form, but can be obtained either by numerical integration or by the Monte Carlo Markov Chain (MCMC) method (Spiegelhalter et al., 1997).We have written a program, using the R2JAGS package in R [R Development Core Team, 2008], that allows the hierarchical Bayes analysis for our new prior using the MCMC method.But the slow computation speed of the MCMC method does not permit its evaluation by repeated use in simulation.Thus, for convenient implementation and evaluation of our hierarchical Bayes method, we approximate the posterior moments of i Y using Laplace's method.This method has been applied by many authors in the context of Bayesian analysis.See [Butar and Lahiri 2002], [Kass Staffey 1989], [Tierney andKadane1986] [Tierney et al., 1989].Following [Kass Staffey 1989] the first order approximation to the posterior mean and variance of i  under the area level model described above in section-2 is given byThe first order Laplace approximation of the posterior moments of i Y can be given by ) Below we derive the analytical expressions of the posterior moments of

Numerical study
This section compares the performance of the approximate hierarchical Bayesian (HB) approach using our new/proposed prior on  with several other existing choices.The HB approach is approximate as it is implemented through Laplace approximation.Moreover, we consider one classical approach in estimating small area means that uses the REML method to estimate the variance component  .Here we study the frequentist properties of the resulting Bayes procedure by a Monte Carlo simulation study under the assumption of a balanced set-up, common mean ( )  model with 0  i f .In the tables, we denote the hierarchical Bayes methods resulting from new proposed prior for variance component  by Proposed HB.For inference about the small area means i  we compare our proposed Bayesian estimator with empirical best linear unbiased predictors (EBLUP) and Heirarchical Bayes (HB(WINBUG)) estimators of  .
Data collected through pilot survey conducted by the Division of Agri-Statistics on estimation of area and yield of apple in District Baramulla has been used for the purpose of our proposed small area estimation.The district Baramulla comprises of 12 blocks viz., Zanigeer, Boniyar, Tangmarg, Wagoora, Sopore, Baramulla, Uri, Pattan, Rohama, Singphora, Rafiabad and Kunzer.Each block consists of different number of villages.A fixed number of five villages were selected at random from each block by simple random sampling.The data set was named apple-1 for analysis and modeling in R/SAS software"s.It has 61 rows and 7 columns.The columns names are Blocks, N, n, Yield, Area, Trees, Actual Yield for names of blocks, total number of villages in each block, number of villages selected from each block, yield of apple from each selected village in metric tons, area under apple orchards and total number of apple trees in each of the selected village, actual yield obtained as per census records.This data set has been used for unit level estimation.

Comparison of different estimators
The performance of different estimators is examined from the accuracy of the point estimates.This is considered through the relative bias and absolute relative bias of different estimators.The different estimators are compared according to four different criteria recommended by the panel on small area estimates of population and income set up by the United States committee on National Statistics (1978), [Datta et al., 2002] [ Ghosh et al., 1996], [ Datta et al., 2012.],[ Pfeffermann 2013.].We compare different estimators on the basis of average relative bias, average squared relative bias, average absolute bias and average squared deviation.Suppose act est m ASD Now using the above four criteria on the apple data set already discussed above.A Comparison of HB estimates using proposed prior with EBLUP estimates and HB estimates with uniform prior is made using above discussed four different criteria and is reported in Table-1.It is clear from the Table-1 that the proposed prior for  performed significantly better than the uniform gamma prior and EBLUP estimates in terms of all the four criterion.This motivates us to argue that our proposed prior is superior even in the cases where the standard hierarchical Bayes estimate with uniform prior seems more appropriate.The value of % ARB is 3.9 per cent in HB (will proposed prior) as compared to 5.61 and 8.73 per cent in HB (with uniform prior) and EBLUP respectively.Now An Empirical comparison of proposed prior for  with REML and WINBUG estimates for all the small areas separately using Percent Absolute Relative Error(ARE) and Absolute Error(AE) has been carried out and is depicted in the following Table 2.       Table-5 reports the relative contribution of the three terms to the posterior variance of i  obtained by using the new/proposed prior.The three columns T 1 , T 2 , T 3 exhibits the relative contribution of term1, term 2, term 3 respectively.From the values in the table we conclude that the contribution of the term which accounts for the uncertainty in estimating the variance component is substantially small relative to the first term which accounts for the uncertainty in the model in estimating the small area means.

Conclusion
In this paper Bayesian implementation of Unit level model is carried out and new prior is proposed on the variance component  .Besides being simple, this prior has two main advantages.It removes the possibility of yielding zero estimates for the variance component; the popular choice of uniform prior on  suffers from this drawback if posterior mode is considered as an estimator.This prior also enjoys good small sample frequentist properties; real agricultural study results justify this conclusion.Also, in order to have closed form expressions of the posterior mean and variance of the true small area mean, Laplace approximation to ratio of integrals, following [Kass and Staffey 1989] is being used.To illustrate the method numerically a real data set on apple has been used and the results showed that the Bayes estimators (with new prior) of small area means have good frequentist properties such as MSE and ARB as compared to other methods viz., EBLUP and HB(WINBUG) estimators.
hierarchical Bayesian model.Combining the likelihood (3. 2) and the prior at level 3, we write the joint distribution of on both  and 2 e  .Integrating out  from (3.5), we find

Fig. 1 Fig. 1 :
Fig.1 shows the comparison of the values of percent relative bias and absolute bias for EBLUP, hierarchical Bayes and proposed Heirarchical Bayes ( LL HB )

Fig. 2
Fig. 2 plots the point estimates of i  for EBLUP, HB (with Uniform prior), HB (with new/proposed prior) against the small areas and also provides a comparison of these values with the actual value of yield obtained in each of the small area.

i  2
Uncertainty in estimating coefficient  3 Uncertainty in model in estimating the variance component A.
of Unit Level model and illustrate the usefulness of this models through an application to horticultural survey data.The paper is organized as follows.In section 2, we first study Unit level models including EBLUP estimators of unit level model.Then in section 3 we propose Bayesian formulation of urea level model with new prior for variance component

Table - 2: Empirical comparison of estimators
From the Table-2 it is evident that HB model with proposed prior exhibits smaller errors and a lower incidence of extreme error than either of the HB (with uniform prior) and EBLUP estimates.

Average Relative Bias Comparison for EBLUP, Heirarchical Bayes(HB) and Proposed HB
1 plots the values of % ARB and AB of EBLUP, hierarchical Bayes and proposed Heirarchical Bayes against small areas.significant disparity is observed among the three estimators.Heirarchical Bayes with the new/proposed prior performed better by

Average Relative Bias Comparison for EBLUP, Heirarchical Bayes(HB) and Proposed HB
the lowest value of both % ARB and AB for each of the small areas compared to EBLUP and hierarchical Bayes (with uniform prior).Mean Square Error (MSE) of Estimators of Variance components for 12 small areas are reported below in providing

Table - 3
reports the different estimates of A for each of the 12 small areas.From the table it is clear that in term of MSE the performance of proposed prior is better than rest of two techniques.

Table - 4: EBLUP and HB estimates and associated Standard Errors Estimators Small Areas EBLUP HB (uniform prior) HB (proposed prior) Actual value
Table-4 reports the EBLUP, HB (uniform prior) and proposed prior estimates and their associated standard errors for all the 12 small areas separately.It is clear from the table that the new/proposed prior for  provides better estimates than EBLUP and HB (uniform prior).In terms of standard error also new/proposed prior provides better estimates than the EBLUP and HB (uniform prior).

Table 5 : Relative contribution of the three terms to the posterior variance of i  using new/proposed prior on 
1 Uncertainty in the model in estimating