Small area estimation under a measurement error bivariate Fay–Herriot model

The bivariate Fay–Herriot model is an area-level linear mixed model that can be used for estimating the domain means of two correlated target variables. Under this model, the dependent variables are direct estimators calculated from survey data and the auxiliary variables are true domain means obtained from external data sources. Administrative registers do not always give good auxiliary variables, so that statisticians sometimes take them from alternative surveys and therefore they are measured with error. We introduce a variant of the bivariate Fay–Herriot model that takes into account the measurement error of the auxiliary variables and we give fitting algorithms to estimate the model parameters. Based on the new model, we introduce empirical best predictors of domain means and we propose a parametric bootstrap procedure for estimating the mean squared error. We finally give an application to estimate poverty proportions and gaps in the Spanish Living Condition Survey, with auxiliary information from the Spanish Labour Force Survey.


Introduction
Data are used as covariates in a Fay-Herriot model to estimate poverty indicators, accounting for the presence of measurement error in the covariates. Polettini and Arima (2015) introduced predictors of small area means based on area-level linear mixed models with covariates perturbed for disclosure limitation.
Adopting a Bayesian approach, Arima et al. (2015) rewrote the Ybarra-Lohr measurement error model as a hierarchical model and introduced Bayes predictors. This last work was later extended by Arima et al. (2017) proposing multivariate Fay-Herriot Bayesian predictors of small area means under functional measurement error. On the other hand, Burgard et al. (2019) followed a likelihood-based approach for extending the Ybarra-Lohr model. They proposed residual maximum likelihood (REML) estimators of the model parameters and introduced empirical best predictors and a mean squared error analytical approximation.
Concerning unit-level models, further contributions to the Bayesian SAE literature on measurement error models are Ghosh et al. (2006), Ghosh and Sinha (2007), Torabi et al. (2009), Datta et al. (2010) and Arima et al. (2012). More recently, Torabi (2013) presented an application of data cloning conducting a frequentist analysis of GLMM with covariates subject to the measurement error model. This paper introduces a three-step bivariate Fay-Herriot model by assuming that the vector of true domain means of auxiliary variables differ from the corresponding vector of direct estimators in a zero-mean multivariate normally distributed random error. The introduced functional measurement error model can be considered as a multivariate adaptation of the Ybarra-Lohr univariate model to a parametric inference setup with multivariate normal measurement errors. The proposed approach can also be considered as the non Bayesian counterpart of the statistical methodology introduced by Arima et al. (2017). Ybarra and Lohr (2008) did not assume the normality of the measurement errors and they proposed a weighted least squared approach to estimate the model parameters. Arima et al. (2017) assumed the multivariate normality of the measurement errors and they considered a Bayesian approach. They simulated the posterior distributions of the model parameters and calculated the hierarchical Bayes predictors of domain means by applying Markov Chain Monte Carlo algorithms. This paper preserves the likelihood modelling of Arima et al., but applies a non Bayesian approach. The main target is to calculate empirical best predictors of domain means and to estimate the corresponding mean squared errors (MSE).
Assuming that measurement errors have a multivariate normal distribution is a natural choice in practice, as due to the central limit theorem, the distribution of the auxiliary variable estimators has an asymptotic multivariate normal distribution. This adaption, besides giving another motivation of the empirical best predictor provided by Ybarra and Lohr (2008), has two mayor advantages. First, we derive Fisher-scoring algorithms for calculating maximum likelihood (ML) and pseudoresidual maximum likelihood (pseudo-REML) estimators of model parameters. Second, we provide a parametric bootstrap procedure for estimating the mean squared errors.
The rest of the paper is organized as follows. Section 2 introduces the measurement error bivariate Fay-Herriot model. Section 3 derives the best predictors of random effects and target domain parameters. It also calculates the MSEs of the best predictors. Section 4 presents the relative efficiency matrix of the best predictors that takes into account that variables are measured with error compared to the corresponding ones that ignores this information. Section 5 proposes a parametric bootstrap procedure for estimating the mean squared error of the empirical best predictors. Section 6 describes the pseudo-REML method for estimating the model parameters. Section 7 carries out simulation experiments to investigate the behavior of the pseudo-REML fitting algorithm, the empirical best predictors and the bootstrap estimator of the mean squared error of the empirical best predictors. Section 8 gives an application to real data where the target is the small area estimation of poverty proportions and gaps in the SLCS, with auxiliary information from the SLFS. Sections 9 summarizes some conclusions. The paper contains three appendixes that are provided as supplementary file. Appendix A gives some auxiliary results for the calculation of the best predictors and their MSEs. Appendix B shows some tables with results of the application to SLCS data. Appendix C presents the Fisher scoring algorithm for calculating the ML estimators of the model parameters.

The measurement error bivariate Fay-Herriot model
Let U be a finite population partitioned into D domains U 1 ; . . .; U D . Let l d ¼ l d1 ; l d2 ð Þ 0 be a vector of characteristics of interest in the domain d and let y d ¼ y d1 ; y d2 ð Þ 0 be a vector of direct estimators of l d calculated by using the data of the target survey sample. The measurement error bivariate Fay-Herriot (MEBFH) model is defined in three steps. The first step indicates that direct estimators are unbiased and follow the sampling model where the vectors e d ¼ ðe d1 ; e d2 Þ 0 $ N 2 0; V ed ð Þ are independent and the 2 Â 2 covariance matrices V ed are known. In most cases, V ed is taken to be the designbased covariance matrix of the direct estimator y d , d ¼ 1; . . .; D.
In the second step the true area characteristic l dk is assumed to be linearly related to p k þ q k explanatory variables, k ¼ 1; . . .;x dkp k Þ be a row vector containing the true aggregated (population) values of p k explanatory variables for l dk and letX d ¼ diag ðx 0 d1 ;x 0 d2 Þ be a 2 Â p block-diagonal matrix with p ¼ p 1 þ p 2 . Let k k ¼ ðk k1 ; . . .; k kp k Þ 0 be a column vector of size p k containing regression parameters for l dk and let k ¼ k 0 . . .; z dkq k Þ be a row vector containing the true aggregated (population) values of q k explanatory variables for l dk and let Z d ¼ diag ðz 0 d1 ; z 0 d2 Þ be a 2 Â q block-diagonal matrix with q ¼ q 1 þ q 2 . Let g k ¼ ðg k1 ; . . .; g kq k Þ 0 be a column vector of size q k containing regression parameters for l dk and let g ¼ g 0 The linking model is where the vectors u d 's are independent of the vectors e d 's. The 2 Â 2 covariance matrices V ud depend on 3 unknown parameters, This manuscript assumes that thex dk 's are unknown random vectors that are predicted from independent data sources. These data sources could be administrative registers or other surveys with larger sample sizes than the target survey. For k ¼ 1; 2, let us define the random measurement error vectors v 0 The third step considers the functional measurement error model where x 0 dk is a row vector containing the unbiased predictors of the components of x 0 dk and the vectors v d and . . .; D, are independent. In most cases, x dk is a vector of direct estimators calculated from data of a different survey and R dk 1 k 2 is taken to be the design-based covariance matrix of vectors x dk 1 and x dk 2 , k 1 ; k 2 ¼ 1; 2.
Let us also define the 2 Â p block diagonal matrices The measurement error bivariate Fay-Herriot (MEBFH) model can be expressed as a single model in the form or in the matrix form We finally assume that matrices V ud , R d and V ed are invertible and that x d , v d , u d , e d , d ¼ 1; . . .; D, are independent, but we only introduce inference procedures conditionally on X. If there are no measurement errors, then the v d 's are null vectors and the bivariate Fay-Herriot (BFH) model is obtained as a special case of (4). The MEBFH model can be considered as a multivariate generalization of the Fay-Herriot model with measurement error studied by Ybarra and Lohr (2008) or by Burgard et al. (2019). This approach was also considered by Arima et al. (2017) in a Bayesian context. Note that the MEBFH model is not a linear mixed model as the matrix B depends on the vector k of model parameters. Therefore, the MEBFH model cannot be expressed in the standard form Y ¼ Xb þ Zu þ e. It holds that Further, it holds that y d j

Best predictors under the MEBFH model
This section derives the best predictors (BP) of the random effects v d and u d and of the target parameter l d . It also calculates variances and expectations of cross products. The proofs of the following propositions are based on the properties of the multivariate normal distribution. We recall that the kernel of the n-variate normal distribution is The first two propositions deal with the BP of v d and its basic properties.
Proposition 1 Under model (4), the best predictor of v d iŝ where Proof The conditional distribution of v d , given the estimators x d and y d , is Therefore, we have We have proved that f ðv d jx d ; y d Þ is a multivariate normal distribution with parameters The following two propositions derive the BP of u d , show that it is predictively unbiased and calculate its variance.
Proposition 3 Under model (4), the best predictor of u d iŝ Proof The conditional distribution of u d , given x d and y d , is h Proposition 4 Under model (4), it holds that E½û bp hThe following two propositions give the best predictor of l d and its MSE. This section ends by defining the empirical best predictor of l d : Proposition 5 Under model (4), the best predictor (MEBFH-BP) of l d iŝ Proof As Proposition 6 Under the model (4), the MSE ofl bp d is Under model (4), the empirical best predictors (MEBFH-EBP) of v d , u d and l d are obtained from formulas (6), (7) and (8) by plugging estimatorsb,r 2 u1 ,r 2 u2 andq in the place of b, r 2 u1 , r 2 u2 and q respectively. The MEBFH-EBP of l d iŝ Section 6 introduces the maximum likelihood and the pseudo-REML of estimators of model parameters. In the application to real data, the MEBFH-EBP of l d is calculated by plugging pseudo-REML estimators.

Relative efficiency matrix of best predictors
An important question is, when using estimated auxiliary variables, how much efficiency can be gained by using the MEBFH-BP instead of the BFH-BLUP. We recall that MEBFH-BP denotes the BP of l d calculated by assuming that the MEBFH model is the true model. Similarly, BFH-BLUP denotes the best linear unbiased predictor of l d calculated by assuming that the BFH model (with no measurement errors) is the true model. The gain in efficiency is measured as the relative reduction of the MSE when using the MEBFH-BP instead of the BFH-BLUP, under the assumption that all model parameters are known and the true model is MEBFH. We first derive the BFH-BLUP of l d and its MSE under the true MEBFH model. If we equate V k to the null matrix in (8), we obtain the BFH-BLUP, i.e. where The prediction error of the BFH-BLUP iŝ Under the MEBFH model, the BFH-BLUP is predictively unbiased. It holds that Under the MEBFH model, with all model parameters known, the MSE of the BFH-BLUP is If all MEBFH model parameters are known, the relative efficiency matrix of the MEBFH-BP compared to the BFH-BLUP is where the division of the 2 Â 2 matrices MSEðl bp d jx d Þ and MSEðl bp 0 d jx d Þ is done component by component. We are mainly interested in the diagonal components RE d11 and RE d11 ; this is to say, in the relative efficiencies when predicting l d1 and l d2 respectively. We further remark that RE d does not depend on x d , g or D.
In what follows, we presents some numerical calculations of the relative efficiencies of the MEBFH model with respect to the BFH model. Let us consider a MEBFH model (4) Table 1 presents the relative efficiencies RE d11 (top) and RE d22 (bottom) for the 45 scenarios 1A,...,9E. All the scenarios take h 1 ¼ 1, h 2 ¼ 3=2. The row cases A, B, C, D, E take s ¼ 1:81; 1:41; 1:01; 0; 61; 0:21 respectively. Scenarios 1-9 take the parameters q s , c and h 3 given by the rows 2, 3 and 4 of Table 1. We observe that the relative efficiencies decrease as s increases from case E to case A. This is to say, the greater the measurement error variance of the auxiliary variables is, the greater gain of efficiency is obtained by using the best predictor based on the MEBFH model. On the other hand, if the measurement errors are negligible (s % 0), then the gain of efficiency almost null.
The greater values of RE d11 and RE d22 are in the column case 5. Therefore, the efficiency gain when using the MEBFH model is smaller when the correlations q s , c and h 3 of the components of the measurement errors v d , the sampling errors e d and the random effects u d respectively, are close to zero. In the limit q s ¼ c ¼ h 3 ¼ 0, we get two independent measurement error univariate Fay-Herriot models and it is not possible to transport information from one component to other. Figure 1 plots the relative efficiencies RE d11 and RE d22 for h 3 2 fÀ0:75; 0:85g and any q s and c. The case h 3 ¼ À0:75 covers the Scenarios 2A and 3A and any other scenario with h 3 ¼ À0:75, À0:5 q s 0:5 and À0:5 c 0:5. The case h 3 ¼ 0:85 contains the scenarios 8A and 9A as individual points in the right-hand printed surfaces. In summary, Table 1 and Figure 1 show some scenarios where the MSE of the MEBFH-BP is around a half of the MSE of the BFH-EBLUP. They also show some other scenarios where the gain of efficiency is rather small. This information is useful to decide in what situations it is more profitable to use the more complex EBP based on the MEBFH model.

Mean squared error estimation
Obtaining an approximation to the MSE of the EBP of l d under the MEBFH model requires tedious calculations. Unlike the case of no measurement errors, the obtained approximation will be quite awkward and not very useful to introduce analytic MSE estimators. This is why we propose applying a parametric bootstrap procedure, like the one introduced by González-Manteiga et al. (2008)

Estimation of model parameters
We consider two methods for estimating g, k, r 2 u1 , r 2 u2 and q under model (4): (1) pseudo-residual maximum likelihood, and (2) maximum likelihood. Both method are based on the distribution of y|X. Appendix C in the supplementary file gives the Fisher-scoring algorithm to calculate the maximum likelihood estimators of the model parameters. However, we we do not implement this last algorithm because it has a greater computational complexity. This section describes the pseudo-residual maximum likelihood method. ML estimators of model parameters have well known asymptotic properties. Under regularity conditions on the auxiliary variables, they are consistent and asymptotically normal. The Fisher-scoring algorithm maximizes the log-likelihood of y|X by solving the corresponding system of nonlinear equation, i.e. first partial derivatives equated to zero. It is a system with p þ q þ 3 equations and the Fisherscoring algorithm inverts matrices of that dimension. If p or q are big, then the algorithm speed decreases and it might have convergence problems when the number of domains D is not big enough.
The REML estimators of the parameters of a linear mixed model are quite attractive. They have similar good asymptotic properties as ML estimators, but their calculation has a lower computational cost because the REML log-likelihood involves only the variance component parameters. Nevertheless, the MEBFH model is not a linear mixed model and the REML method is thus not applicable. This is why, we introduce pseudo-REML approach by treating the components of matrix B d in (4) as known constants. In that case, the MEBFH model becomes a linear mixed model and the REML method can be applied yielding to the maximization of the derived REML log-likelihood. However, the components of B d depends on the unknown vector of parameters k and, therefore, we are not applying the REML method but a pseudo-REML approach. The small sample properties of the pseudo-REML are empirically investigated in Simulation 1.
Let us define T ¼ ½Z; (4) can be written in the form The pseudo-REML log-likelihood of model (12) is we calculate the first partial derivatives of l reml with respect to h ' , i.e.
Therefore, the score vector for ' ¼ 1; 2; 3 is For a; b ¼ 1; 2; 3, the second partial derivatives of the REML log-likelihood function are where the last equality follows from the fact that V ' is symmetric, ' ¼ 1; 2; 3. By changing the sign, taking expectations and applying PT ¼ 0, PV ¼ I À V À1 TQT 0 and the formula we get the components of the Fisher information matrix, i.e.
Therefore, the Fisher information matrix is The trace of PV a PV b can be calculated as The pseudo-REML Fisher-scoring algorithm is  The asymptotic distributions of the REML estimatorsĥ andb, can be used to construct ð1 À aÞ-level confidence intervals for the components h ' of h and b j of b, i.e.
where F À1 ðĥ;bÞ ¼ ðm ab Þ a;b¼1;2;3 , ðT 0 V À1 ðĥ;bÞTÞ À1 ¼ ðg ab Þ a;b¼1;...;pþq and z a is the a-quantile of the N(0, 1) distribution. Forb j ¼ b 0 , the p-value for testing the hypothesis We remark that we have changed the notation in (13) and (14), where b j denotes the j-th component of the vector b and not the vector of regression parameters of the j-th category.

Simulations
This section presents simulation experiments for assessing the performance of the fitting method, the EBP estimator, and the MSE estimator. The objective is to show how the new methodology works in realistic (not extreme) scenarios. In practice, some explanatory variables can be taken from an auxiliary survey that is different from the target survey and that has bigger sample sizes. This is the case of the SLCS (target survey) and the SLFS (auxiliary survey). The second one has bigger sample sizes than the first. Therefore, it is natural to choose scenarios where the variances of the measurement errors are lower than the variances of the sampling errors and lower than the variances of the random effects.
In the calculation of Sect. 4, the variances of the sampling errors are r 2 e1 ¼ r 2 e2 ¼ 1 and the variances of the random effects are r 2 u1 ¼ h 1 ¼ 1 and r 2 u2 ¼ h 2 ¼ 3=2 ¼ 0:66. Therefore, it is quite natural to choose scenario E, where the measurement error variances are s d11 ¼ s d22 ¼ 0:21 Concerning the random effect correlations, we are mainly interested in the case h 3 ¼ corr ðu d1 ; u d2 Þ ¼ 0:45 because h 3 is positive in the application to real data. This is why, we consider that the scenario 6E is the most close to the application to real data presented in Sect. 8.
For the sake of completeness, we also run simulations under the scenarios 5E and 4E with correlations h 3 ¼ 0:05 and h 3 ¼ À0:55.

Simulation 1
The target of Simulation 1 is to check the behavior of the pseudo-REML Fisherscoring algorithm for fitting the MEBFH model. The steps of Simulation 1 are 1. Generate z dk , x dk , d ¼ 1; . . .; D, k ¼ 1; 2. 2. Repeat I ¼ 1000 times (i ¼ 1; . . .; 1000)   We also note that root-MSEs decrease slowly as the number of domains increases. This is somewhat reasonable, since increasing D also increases the number of quantities l d to be predicted.

Simulation 3
Simulation 3 investigates the performance of the parametric bootstrap estimator of the MSE of the EBPs. For D ¼ 100, the steps of Simulation 3 are ÃðiÞ dk À mse dk Þ: 5. For k ¼ 1; 2, calculate the averages across domains of the relative performance measures, i.e. RB Ã k ¼ 1

Estimation of poverty proportions and gaps in Spanish provinces
This manuscript presents an application of the MEBFH model to the estimation of poverty proportions and gaps in Spanish provinces by sex. Spain is divided in 52 provinces (including the autonomous cities of Ceuta and Melilla) leading to D ¼ 104 target domains (provinces crossed by sex) of known sizes N d , d ¼ 1; . . .; D.
Let z dj be the normalized net annual income of the household where the individual j of domain d lives. Let z 0 be the poverty line, so that individuals with z dj \z 0 are considered as poor. The main goal of this section is to jointly estimate poverty proportions and gaps where y d1j ¼ Iðz dj \z 0 Þ, Iðz dj \z 0 Þ ¼ 1 if z dj \z 0 , Iðz dj \z 0 Þ ¼ 0 otherwise and y d2j ¼ z À1 0 ðz 0 À z dj ÞIðz dj \z 0 Þ. The Spanish Statistical Office calculates z dj by summing up the net annual incomes of the household members and by dividing by its normalized size. Later, the same value of the normalized net annual income of the household is assigned to all the individuals j of the household. Therefore z dj is constant within the household. The aim of normalizing the household income is to adjust for the varying size and composition of households. The definition of the total number of normalized household members is the modified OECD scale used by EUROSTAT. This scale gives a weight of 1.0 to the first adult, 0.5 to the second and each subsequent person aged 14 and over and 0.3 to each child aged under 14 in the household h. The normalized size of a household is the sum of the weights assigned to each person. So the total number of normalized household members is where H dh ! 14 is the number of people aged 14 and over and H dh\14 is the number of children aged under 14. Following the standards of the Spanish Statistical Office, the poverty threshold is fixed as the 60% of the median of the normalized net annual incomes in Spanish households. The Spanish poverty threshold (in euros) in 2008 is z 2008 ¼ 7488:65. We deal with data from the SLCS of 2008 with sample size 35967. This is an annual survey where the planned domains are the regions (autonomous communities), so that sample sizes are selected to obtain precise direct estimates at the region level. As Spain is hierarchically partitioned in regions, provinces, counties (comarcas) and municipalities, estimating province-sex poverty proportions is a small area estimation problem. The direct estimators of the size N d , the total Y dk ¼ P N d j¼1 y dkj and the mean where s d is the domain sample of size n d and the w dj 's are the official calibrated sampling weights which take non response into account.
where the case k 1 ¼ k 2 ¼ k denotes estimated variance, i.e.v ar p ð Y dir dk Þ 1 cov p ð Y dir dk ; Y dir dk Þ. The last formulas are obtained from Särndal et al. (1992), pp. 43, 185 and 391, with the simplifications w dj ¼ 1=p dj , p dj;dj ¼ p dj and p di;dj ¼ p di p dj , i 6 ¼ j in the second order inclusion probabilities.
We take data from the SLFS of 2008 to construct the data file of aggregated auxiliary variables. The SLFS is a quarterly survey with provinces as planned domains. Within each quarter, the SLFS sample sizes at the province level are fixed a priori. They are selected high enough to have precise direct estimates. By putting together the data files of the 4 quarters, the SLFS direct estimators of means at domain (province crossed by sex) level are even more precise for 2008 than using only the data of one quarter. Additionally, by doing this aggregation all the calculated SLFS direct estimators of domain means have non zero estimated variances.
The file of auxiliary variables is constructed from the SLFS data of 2008 (SLFS2008). It contains the direct estimates of the domain means by the categories of the considered auxiliary variables. It also contains the variance and covariance estimates of the calculated direct estimators. The auxiliary variables are: Nationality, with categories Spanish and Foreigner. Education, with categories LowEdu (less than secondary level completed), HighEdu (secondary or superior level completed) Age, with categories age1 ( 15), age2 (16-24), age3 (25-49), age4 (50-64), age5 ( [ 64).
Labour situation, with categories 15, Employed, Unemployed, Inactive. Table 5 presents the pseudo-REML estimates of the regression parameters of the selected MEBFH model. It also contains the corresponding p-values. We note that domains with higher proportions of Spanish or unemployed people tends to have higher poverty proportions. On the other hand, the level of education is negatively related with the poverty proportion. We note that foreigners tend to live and work in rich provinces. Therefore, the obtained results are in agreement with the Spanish socioeconomic situation of 2008. Table 6 contains the estimates of the variance component parameters and the corresponding asymptotic 95% confidence intervals. Figure 4 (left) plots the MEBFH model residuals of poverty proportions versus the corresponding EBPs. Figure 4 (right) plots the MEBFH model residuals of poverty gaps versus their EBPs.
We observe that the residuals present greater variability on the right hand side of the OX axis. This is a natural phenomenon, as the EBPs have to smooth the behavior of the direct estimators. The plots show that EBPs tend to be smaller than direct estimates when direct estimates are large. This is due to two main reasons. First, it is part of the smoothing work that EBPs have to do. Second, the variance of the directly estimated proportion is approximately ð1 À n d =N d Þ=ðn d À 1Þ Y dir dk ð1 À Y dir dk Þ, which takes its maximum at Y dir dk ¼ 0:5, for n d ; N d fixed. Allowing in this context area specific V ud is an interesting future research task, that possibly could further reduce the MSE of the EBP. Figure 5 plots the EBPs and direct estimates of poverty proportions (left) and gaps (right) for men and women. We observe that both estimators tend to coincide as the sample size increases. Figure 6 plots the estimated model-based root-MSEs of the EBPs and the estimated design-based root-MSEs of the direct estimators of poverty proportions (left) and gaps (right) for men and women. The first root-MSE is estimated by  smallest are those situated in the north and east. On the other hand, the Spanish provinces with higher poverty proportions are those situated in the center and south. Appendix B in the supplementary file contains 2 tables with basic numerical results. To give a brief summary of the obtained numerical results, we order the Spanish provinces (including the cities of Ceuta and Melilla) by sample size and we select one in five. Table B.1 presents the EBPs and direct (dir) estimates of poverty proportions for men and women respectively. Table B.2 presents the EBPs and direct (dir) estimates of poverty gaps for men and women respectively. These tables also give the estimated root-MSEs (rmse.eblup and rmse.dir), for both types of estimators. The column N d presents the estimates of the domain sizes calculated by using the data of the 4 quarters of the SLFS in 2008. The columns province and n d contain the province name and the sample size respectively.  In many applications the auxiliary information used in bivariate Fay-Herriot models is not measured exactly. Under this setting, the EBLUP based on the bivariate Fay-Herriot model is not the EBP anymore, as this model assumes the exact measurement of the auxiliary variables.
By calculating relative efficiencies, we showed that not taking into account the measurement errors may lead to predictions of target parameters with greater mean squared errors. Therefore, we extended the bivariate Fay-Herriot model by allowing for multivariate normal distributed random error on the auxiliary variables. This reflects the typical situation of estimated auxiliary information. Both, the mean squared error and the bias of the new EBP are reduced with respect to the classical EBLUP.
For fitting the new model, the pseudo-REML estimation procedure showed to be stable. Further a second fitting algorithm was introduced but not implemented. In the case of exactly measured auxiliary variables, the measurement error bivariate Fay-Herriot model reduces to the classical bivariate Fay-Herriot model and, therefore, they basically give the same predictions. Finally, we recommend to use the proposed measurement error Fay-Herriot model when the auxiliary variables are estimated.