Small area estimation of socioeconomic indicators for sampled and unsampled domains

Socioeconomic indicators play a crucial role in monitoring political actions over time and across regions. Income-based indicators such as the median income of sub-populations can provide information on the impact of measures, e.g., on poverty reduction. Regional information is usually published on an aggregated level. Due to small sample sizes, these regional aggregates are often associated with large standard errors or are missing if the region is unsampled or the estimate is simply not published. For example, if the median income of Hispanic or Latino Americans from the American Community Survey is of interest, some county-year combinations are not available. Therefore, a comparison of different counties or time-points is partly not possible. We propose a new predictor based on small area estimation techniques for aggregated data and bivariate modeling. This predictor provides empirical best predictions for the partially unavailable county-year combinations. We provide an analytical approximation to the mean squared error. The theoretical findings are backed up by a large-scale simulation study. Finally, we return to the problem of estimating the county-year estimates for the median income of Hispanic or Latino Americans and externally validate the estimates.


Introduction
Socioeconomic indicators, such as the median income of sub-populations, are key both for policy recommendations and policy evaluation. Regional indicators of income, poverty, employment, or well-being are omnipresent in current projects and research. To better reflect regional heterogeneity, their focus has shifted to deeper regional levels.
Regional information is normally published on an aggregated level. The estimation of these aggregates is usually based on survey data. Although the demand for more detailed regional social indicators has increased, the underlying surveys tend to focus on higher regional or national level so that the survey costs do not become too high. Even though survey data provide reliable direct estimates at these levels, for more detailed regional levels they usually do not. At a regionally lower level, direct estimates are usually associated with high standard errors or unavailable. Estimates can be unavailable if the region is unsampled and thus no direct estimate could be given. In addition, the publication of regional estimates can be suppressed when the associated standard errors are high. For the investigation of regional indicators a researcher is therefore usually confronted with regional aggregates which are associated with high standard errors and contain a non-negligible proportion of unpublished data.
Small area estimation (SAE) methods are increasingly used to deal with highly volatile direct estimates on regional level. For a comprehensive overview, see Rao and Molina (2015). The key idea behind SAE techniques is to borrow strength by combining regional indicators in a common model framework. Within this model, additional related information, such as registration data, can be exploited.
The model-based approach allows the introduction of best predictors (BP) that minimize the mean square errors (MSE) in the class of unbiased predictors. Since the BPs depend on the model parameters, substituting them for appropriate estimators gives the empirical BPs (EBP) that stabilize the estimates in small domains.
For aggregated univariate target data, the most prominent model is the Fay-Herriot (FH) model by Fay and Herriot (1979). Many extensions were made to the FH predictor to meet different practical problems. Inter alia, Prasad and Rao (1990) and Datta and Lahiri (2000) propose MSE estimators for the FH predictor, Li and Lahiri (2010) and Yoshimori and Lahiri (2014) introduce new adjusted maximum likelihood fitting methods, Jiang and Tang (2011) study the influence of the fitting algorithm in the empirical best prediction, Molina et al. (2015) derive preliminary testing predictors, Moura et al. (2017) modified the basic model to analyze skewed business survey data, Ybarra and Lohr (2008), Bell et al. (2019),  and  study the effect of measurement errors in the covariates, Pratesi and Salvati (2008), Gonzáalez-Manteiga et al. (2010), Articus and Burgard (2014), and Morales et al. (2015) allow for a heterogeneous dependency structure in the FH model, Esteban et al. (2012) and Marhuenda et al. (2013) estimate small area poverty proportions under temporal and spatiotemporal Fay-Herriot models, respectively.
For aggregated multivariate target data, a widely employed model is the multivariate FH model. Fay (1987) and Datta et al. (1991) investigated the gain of precision achieved by the multivariate modeling. Datta et al. (1996) employed a multivariate FH model for obtaining hierarchical Bayes predictors. González-Manteiga et al. (2008) considered a multivariate FH model with a common domain random effect for the target vector, Arima et al. (2017) and Burgard et al. (2020) study multivariate measurement errors FH models, Porter et al. (2015), Benavent and Morales (2016), Ubaidillah et al. (2019), Esteban et al. (2019) and Benavent and Morales (2021) investigate and give further applications of multivariate FH models. Many other authors have studied further variants of the FH model and the multivariate FH model adapted to different setups.
In the context of survey sampling, there are two main types of missing data. First, missingness due to non-response refers to situations where data were planned to be collected by the nature of the sampling design, but failed to be collected due to some kind of response mechanism. This may occur because individuals in the sample refuse or fail to respond, or because of processing issues. The response mechanisms are of special concern in voluntary surveys where ignoring the response mechanism could lead to biased direct estimates. The treatment of this kind of missing data is therefore of great interest for applied statisticians as shown by Matei and Ranalli (2015) and Nguyen and Zhang (2020) in their recent studies on latent modeling approaches and reweighting methods for nonresponse. More generally, the book of Longford (2005) gives an introduction to missing data modeling and imputation methods. In the present study, we do not consider any kind of non-response, but deal solely with aggregate information such as domain-specific direct estimates which might have been adjusted to nonresponse by the statistical agencies.
Second, missing data can occur from the sampling design, if the design does not allocate samples to domains of interest. These domains could, for example, be a cross-combination of small geographic units and demographic characteristics such as age classes. In this case, the domain-specific sample sizes are random and can be zero or so small that corresponding direct estimates are not published due to the estimated variances being too high. This is the kind of missing information we are considering in this paper. The problem of estimating small area indicators with missing data, i.e., with unsampled domains or simply with unavailable data, has been scarcely treated in the literature. For unit-level data, Longford (2004) did some contributions related to multivariate shrinkage estimators, where the estimation is integrated with a multiple imputation procedure. For area-level small area models, to our knowledge, the existing SAE approaches using the FH model give empirical best predictors (EBP) only for domains with observed direct estimates. For unavailable direct estimates, the existing variants of the FH model only provide synthetic estimates with vague mean squared error approximations.
As an example, we can take a look at Articus et al. (2020) where the FH model is applied to local-level rental markets based on direct estimates from the German Microcensus. We can see that in the application of FH models missing direct estimates either appear as blank spots on a map (Articus et al. 2020, Figure 6) or are filled by synthetic predictions (Articus et al. 2020, Figure 7). We attempt to fill this research gap with the introduction of a new model based on the FH model and bivariate modeling, called MBFH. It is an extension of the bivariate FH model, see, e.g., Rao and Molina (2015), Section 4.4.1, that allows for partially missing direct estimates. On the basis of this model, we derive the corresponding empirical best predictor under missing values (MBFH-EBP). In contrast to the FH or bivariate FH model, the MBFH model provides best predictors for both missing and observed values. We furthermore provide analytical mean squared error approximations of the MBFH-EBPs.
This best predictor is applicable to any indicator that is predictable by the Fay-Herriot model. For example, the dependent variables may consist of direct estimators labor market indicators (total of employed and unemployed people and unemployment rates), living conditions indicators (head count ratio, poverty gap, poverty housing, average per capita income), family budget indicators (per capita expenditure in food or housing), and others. In addition, the new predictor is applicable to setting where the target variable has missings, but a second correlated variable is observed in all domains.
The choice of auxiliary data can be made for each variable separately. We illustrate the use of the proposed method with an application to US ACS data at the county level, where direct estimates of median income for Hispanics or Latin Americans are partially missing.
The manuscript is structured as follows: Sect. 2 describes the problem of interest and the data which we use for an illustration. Section 3 provides the methodological foundation. Section 3.1 introduces the bivariate Fay-Herriot model, which is the basis for the development of the EBP theory under partially missing target estimates, the proposed MBFH predictor. Section 3.2 divides the set of domains in three groups depending on the missing structure of the direct estimates and gives the corresponding MBFH-EBPs. Section 3.3 gives an approximation to the MSE of the MBFH-EBP and proposes an explicit-formula estimator. With a model-based simulation, Sect. 4 validates the theoretical results and empirically investigates the introduced MBFH-EBPs and MSE estimators in different settings. Section 5 applies the new MBFH predictor to the publicly available US county-level ACS data on median income of Hispanic or Latino Americans. Section 6 presents a short summary and outlook. This contribution is the extension of a related working paper . Finally, the manuscript has four appendixes in the supplementary material. Section 1 gives algorithms for calculating the maximum likelihood and residual maximum likelihood estimators of the model parameters. Section 2 contains proofs of the derivation of best predictions under the new model. Section 3 contains the mathematical derivations for approximating the MSE of the MBFH-EBP. Section 4 presents a parametric bootstrap procedure for estimating the mean squared error of the MBFH-EBPs. Section 5 shows the MSE of the synthetic Fay-Herriot predictor for missing direct estimates.
Section 6 contains additional results from the model-based simulation study with 5% and 10% missing direct estimates.

3
Small area estimation of socioeconomic indicators for sampled… 2 The problem of interest

Aim
To illustrate the new MBFH predictor, we take a look at the publicly available county-level US ACS estimates of the median income of Hispanic or Latino Americans. Income is strongly related to the concepts of poverty and well-being and thereby has an essential impact on the regional distribution of resources. The US Office of Management and Budget (OMB) requires federal agencies to use at least two ethnicities in data collection and reporting: 'Hispanic or Latino' and 'not Hispanic or Latino'. Most official US survey publications, including those on income and poverty, pay extra attention to people who consider themselves as Hispanic or Latino Americans, see, e.g., Guzman (2019) and Semega et al. (2019). Therefore, we consider regional estimates of the median income of Hispanic or Latino Americans.
The ACS estimates on county-level for the years 2010 and 2011 are partly estimated with high standard errors or are not published for certain counties and years. We can use additional publicly available data from the US Census Bureau to construct bivariate FH models based on these estimates. The statistical issue is that we cannot apply the statistical methodology based on the multivariate FH model to the domains with missing direct estimates. This paper introduces a new model-based approach to address this problem. The approach makes use of the fact that some regional estimates missing in 2010 are available in 2011 and the other way around. Intuitively, the random effect correlation of the median income between two years is expected to be highly positive. As the simulation study in Sect. 4 shows, in such a situation an application of the introduced MBFH predictor is profitable for stabilizing volatile estimates and predicting missing estimates. In the following, we describe the county-level US ACS data and the collection and choice of publicly available auxiliary data.

Data description
We use the freely available US American Community Survey (ACS) data. Detailed information about the ACS is given in US Census Bureau (2014). The US Census Bureau provides aggregated county-level data of ACS 1-year estimates with associated margins of error. We take the median annual income (dollars) Hispanic or Latino origin (of any race) (HC02_EST_VC12) in 2010 and 2011 as variables of interest. 1 This is the target vector of the statistical study.
We partition our finite population, the USA, into D = 3141 ( d = 1, … , D ) US counties, our domains of interest. There are publicly available county-level data from the US Census Bureau which can be used as auxiliary data. 2 As 2010 and 2011 are temporary close, we choose the same auxiliary data for both variables, HC02_EST_VC12 in 2010 and 2011. Among the available auxiliary data, we chose the final model considering correlation patterns to the variables of interest, resulting coefficients, and model diagnostics. The chosen auxiliary variables are: Intercept, death rate in period 7/1/2010 to 6/30/2011 (RDEATH2011), and civilian labor force unemployment rate 2010 RTE (CLF040210D).
The To validate the results of the MBFH-EBP, especially the prediction of missing direct domain estimates, we need comparable county-level data. From the ACS also 5-year direct estimates are available for the variable of interest (HC02_EST_VC12). 3 ACS 5-year estimates pool data from the last 5 years, such that direct estimates in 2012 refer to ACS data from 2008-2012 and are available for more counties than the 1-year estimates. It is not recommended to directly compare overlapping ACS datasets such as ACS 1-and 5-year data. As we, however, want to evaluate whether the MBFH-EBPs of missing ACS 1-year estimates are realistic, the ACS 5-year direct estimates are chosen as benchmarks. They are available for many counties where ACS 1-year direct estimates are missing, but MBFH-EBPs can be computed.
Although ACS 5-year estimates of HC02_EST_VC12 are available for some counties where 1-year estimates are missing, they do not cover all areas either. We therefore use an additional dataset for validation, the Census-ACS 2010 estimates. Census estimates for variable HC02_EST_VC12 are not available. We therefore choose Census estimates of variable median household income in the past 12 month

3
Small area estimation of socioeconomic indicators for sampled…

The bivariate Fay-Herriot model
� be a vector of characteristics of interest in domain d and let � be a vector of direct estimates of d calculated by using the data of the target survey sample. The bivariate Fay-Herriot model is defined in two stages. The first stage indicates that direct estimators y d , ∀d ∈ {1, … D} , are unbiased and follow the sampling model where the vectors e d = (e d1 , e d2 ) � ∼ 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 design-based covariance matrix of direct estimators y d , ∀d ∈ {1, … , D} . The covariance matrices V ed are In the second stage, the true domain characteristic dk is assumed to be linearly related to p k explanatory variables, k = 1, 2 , d ∈ {1, … , D} . Let x � dk = (x dk1 , … , x dkp k ) be a row vector containing the true aggregated (population) values of p k explanatory variables for dk and let where "col" is the matrix operator stacking by columns. We finally assume that u d ,

) is a reparametrization of Model 3 introduced by Benavent and Morales (2016).
Let

Prediction with missing target values
Let us assume that some of the y dk are missing and partition the domains into three groups: If the BFH model (3.3) holds for d ∈ {1, … , D} and the missing data obey scheme {1, … , D} = 1 ∪ 2 ∪ 3 , we say that target vectors y d obey a missing data BFH (MBFH) model. If the MBFH model holds, then The supplementary material gives fitting algorithms to calculate the maximum likelihood (ML) and residual maximum likelihood (REML) estimators of the MBFH model parameters.
In a real situation where the target data follow a MBFH model, the BFH model is strictly applicable to 3 , but not to 1 or 2 .

Analytic approximation of the mean squared error of the MBFH predictor
This section gives an analytical approximation of the MSE of the MBFH-EBP for each of the three considered groups of domains. The corresponding mathematical derivations are presented in the supplementary material. Alternatively, the supplementary material also introduces a parametric bootstrap procedure for estimating the MSE of the MBFH-EBPs.

Empirical best predictors in domains of groups 1 and 2
As the estimators for groups 1 and 2 are somehow symmetric, we only present those corresponding to group 1 : where The derivatives of matrix Φ d1 ( ) with respect to , = 1, 2, 3 , are The derivatives of h d ( , ) with respect to kj and , k = 1, 2 , j = 1, … , p k , = 1, 2, 3 , arê The p k × 1 vectors containing the derivatives with respect to k , k = 1, 2 , are and the corresponding The 3 × 1 vectors containing the derivatives with respect to are and the corresponding 3 × 3 matrices are Similarly as Prasad and Rao (1990), Datta and Lahiri (2000), and Das et al. (2004), Note that when using ML instead of REML estimation an additional bias correction term has to be considered, see Datta and Lahiri (2000) for further details.

Empirical best predictors in domains of group 3
Let us consider the group 3 . We write .
(4.1) for variable 2 ∀d ∈ 1 . iv. Calculate the MSE predictor for FH-EBLUP (cf. Prasad and Rao 1990;Datta and Lahiri 2000). v. Calculate the MSE predictor for FH-SYN using the derivations in Appendix C. For variable 1 D 0 = 1 ∪ 3 , D 1 = 2 . For variable 2 D 0 = 2 ∪ 3 , D 1 = 1 . The simulation results are presented in Figs. 1, 2, 3, and 4. We show the results for 20% missing values in both variables of interest. The same plots with 5% and 10% missings look very similar and are displayed in Sect. 6 of the supplementary material. The figures differ with respect to the performance measure shown and their underlying sampling error correlation. They show the performance of the predictors MBFH-EBP, FH-EBLUP, and FH-SYN and their corresponding MSE estimates. The ML and REML methods are standard estimation methods in linear mixed models. However, REML estimators of the variance component parameters tend to have a lower bias and have a lower computational cost. This is why we obtain model parameters using the REML Fisher-scoring algorithm. The results under ML are similar and therefore not shown.

For
The performance measures are depicted separately for the three domain groups, 1 , 2 , and 3 , and the two variables of interest, resulting in six panels. The grayshaded panels indicate domain-variable combinations where direct estimates are missing. Each panel presents the results for different underlying correlations of random effects . This way the interplay of the estimators efficiency with the correlation of sampling errors and random effects can be retrieved. In each panel, a boxplot of the RRMSE or RBias over the corresponding domains is drawn for the different random effect correlations. In panels with observed direct estimates, the boxplots are shown for FH-EBLUP and MBFH-EBP, whereas in panels with missing direct estimates, i.e., in gray-shaded panels, the boxplots are shown for FH-SYN and MBFH-EBP. For facilitating the comparison, in our application in Sect. 5 we have D = 762 , D 1 = 78 , D 2 = 58 , D 3 = 626 such that about 10% of the first and 7.6% of the second variable of interest are missing. Furthermore, in the application the sampling error correlation is ed12 = 0 , and the correlation of random effects is estimated at 0.97. We first evaluate the performance of the predictors MBFH-EBP, FH-EBP, and FH-SYN via their RRMSE. Figures 1, 2, and 3 show the corresponding RRMSE for different underlying sampling error correlation ed12 ∈ {0.6, −0.3, 0} . The nonzero sampling error correlations in Figs. 1 and 2 correspond to situations where both variables stem from the same survey. A zero sampling error correlation as in Fig. 3 applies if both variables stem from independent surveys. It is visible for all three Figs. 1, 2, and 3 that in terms of RRMSE the introduced MBFH-EBP is at least as efficient as the respective FH-EBLUP or FH-SYN. This observation holds for positively, negatively, and uncorrelated sampling errors. Thus, the application of the MBFH instead of the FH estimator is profitable for various combinations of sampling error and random effect correlations.
The prediction of missing values is visible in the gray-shaded panels. There, the efficiency gain of the MBFH-EBP over the FH-SYN in terms of RRMSE is especially high for high random effect correlation, positive or negative. This finding applies to all sampling error correlations, i.e., Figs. 1, 2, and 3 with ed12 ∈ {0.6, −0.3, 0} . Thus, for the prediction of missing values the application of MBFH instead of FH is profitable as long as the absolute value of the underlying random effect correlation gets away from zero.
The performance of the MBFH-EBP for domains with no missing values is visible in the panels of group 3 . For these, in Figs. 1, 2, and 3, the efficiency gain of MBFH-EBP over FH-EBLUP is especially high when the correlation of random effects and sampling errors is of opposite sign and high magnitude. This coincides with the findings in Datta et al. (1999) for unit-level multivariate small area models without missing values.
Next, we take a look at the corresponding MSE estimates. Figure 4 presents the simulation results of the relative Bias (RBias) of the different MSE estimators under REML and sampling correlation ed12 = −0.3 for varying random effect correlations . The MSE estimates of MBFH-EBP lead good results for both missing and observed domains. The MSE estimates of FH-SYN, which derivation is shown in Appendix C, also lead good results. The relative bias of the MBFH-EBP MSE is around zero in most panels and for all random effect correlations. We can see that in group 3 there is a slight positive bias in the mse estimate for = 0.9 with ed12 = −0.3 . This bias is not visible when ed12 = 0 as in the application.
The MBFH-EBP is best for both observed and missing estimates when (1) there are many domains with both variables observed while in case a variable is missing the other variable is observed and (2) the correlation of random effects and sampling errors is of opposite sign and high magnitude. One would, for example, expect the random effects of a variable in two consecutive years to be highly positively correlated, e.g., in the application in Sect. 5, we focus on the median income of Hispanic or Latino Americans in two consecutive years). On the other hand, random effects of variables like employment and unemployment would be expected to be highly negatively correlated. The correlation of sampling errors is determined by the sampling design. As we see from the MSE estimation in Fig. 4 in case a researcher is unsure whether applying the MBFH instead of the FH estimator is improving the predictions, a comparison of the MBFH and FH MSE estimates is recommended. For the selection of variables of interest, we would therefore advise researchers to pay special attention to the missing patterns of the two variables before considering the correlation patterns.

Model results
Let y d = y d1 , y d2 � be the vector of direct estimates of the US county-level median annual incomes of Hispanic or Latino origin people in 2010 and 2011. We assume that y d follows model (3.3) where  Columns with labels beta, std.error, t-value, and p-value, respectively, contain the values of the regression parameter estimator, the standard error, the t-statistic and the p-value, respectively. The estimated coefficients are very similar for HC02_ EST_VC12 in 2010 and 2011 and highly significant. They furthermore seem plausible considering that counties with higher death rate and higher civilian labor force unemployment rate tend to have lower values of HC02_EST_VC12. Any thematic conclusions from the few freely available auxiliary data are, however, limited, e.g., due to high correlations among variables. The 2 × 2 covariance matrix V ud depends on three unknown parameters, , 2 = 2 u2 , and 3 = , i.e., As we consider the same variable in two consecutive years, we expect the correlation of random effects to be highly positive. Table 3 contains the estimates of the variance component parameters and the corresponding asymptotic 95% confidence intervals. The estimated random error correlation is highly positive.  Figure 6 plots the histograms of the standardized residuals. From the model, the residuals are expected to be normally distributed with mean zero and both figures show that the mass of the residuals is distributed near and around zero. At the same time, we detect some major outliers. In 2010 and 2011, about 2% of the standardized residuals have an absolute value greater 3. A further treatment of these outliers is necessary, but beyond the scope of this paper where we use the data for an illustration of the presented MBFH model. Unfortunately, we do not have access to new auxiliary variables that could explain the behavior of the target variables in the few domains where the model performs poorly. For further studies, an investigation of the outliers and of a robust version of the proposed model would be interesting. We refer to Sinha and Rao (2009) for general information on robust SAE, Schmid and Münnich (2014) for robust SAE including spatial effects, and Baldermann et al. (2018) for the additional consideration of spatial non-stationarity. For an implementation of robust FH models in R, see package rsae (Schoch 2014).
, ∀d ∈ {1, … , D}.  counties in which ACS 1-year direct estimates are available. Figure 7 plots the standard error of ACS 1-year direct estimates versus the difference between the direct estimates and MBFH-EBPs of HC02_EST_VC12 in 2010 and 2011. As the FH estimator is a shrinkage estimator, it is expected that the more reliable the direct estimates, and thereby the lower the ACS 1-year standard error, the closer are MBFH-EBPs and direct estimates. This shrinkage property is visible in Fig. 7 for the MBFH. Due to the shrinkage property of the MBFH, the root MSE of the MBFH-EBPs is expected to be close to the ACS standard errors for reliable direct estimates. For highly volatile direct estimates, on the other hand, it is expected to be much lower. The total number of persons HIS010210D 5 in a county is expected to be an indicator of the reliability of the direct estimates, the more persons in a county, the more reliable the resulting estimate. Figure

Validating predictions of missing direct estimates
To validate the MBFH-EBPs also for missing direct estimates, ACS 5-year estimates of the variable of interest and Census 2010 estimates of a similar variable are used. Table 4 shows the quantiles of the relative difference (in %) of the ACS 1-year direct estimates and the MBFH-EBPs to the ACS 5-year direct estimates of HC02_ EST_VC12. We compare 1-year direct estimates and EBPs in 2010 and 2011 to ACS 5-year direct estimates in 2008-2012 and 2009-2013, respectively. The results for the MBFH-EBPs are shown both for domains with available direct estimates and those without direct estimates, but were MBFH-EBPs could be computed. Due to their larger sample sizes the ACS 5-year estimates are less volatile than the ACS  1-year estimates and available for more counties (Table 4). Even though differently aggregated ACS estimates are not directly comparable, a comparison between ACS 1-year direct estimates and MBFH-EBPs to the ACS 5-year estimates should give a suitable image of the reliability of the MBFH-EBPs. As can be seen in Table 4, for domains with given direct estimate the quantiles of the relative differences of the MBFH-EBPs are smaller than those of the ACS 1-year direct estimates. Furthermore, the quantiles of the relative difference of the MBFH-EBPs of domains with no available direct estimate are close to those of the ACS 1-year direct estimates. The proximity of quantiles is more visible in 2011 than in 2010. This indicates that for both domains with available and missing ACS 1-year direct estimate, the MBFH predictor accomplishes more realistic predictions in 2011 than in 2010. Figure 9 shows the different estimates exemplary for the states Indiana and Ohio in 2010 (left plots) and 2011 (right plots). In rows one, two and three are the ACS 1-year direct estimates, the ACS 5-year direct estimates, and the MBFH-EBPs, respectively. In both years, many ACS 1-year and some ACS 5-year direct estimates are missing, indicated by the white-shaded domains. Framed counties indicate those with ACS 1-year direct estimates missing in one year, but available in the other. For these counties the MBFH model provides EBPs in contrast to the commonly used univariate FH model. Therefore, in row three the framed counties are filled by MBFH-EBPs. The plot confirms the results from Table 4  For the analysis of the county-level median income of Hispanic or Latino Americans, the MBFH estimator shows promising results. The validation indicates that the MBFH-EBPs are realistic for both observed and missing direct estimates. It should be noted that this result is already achieved under the rather small number of freely available auxiliary data. Even better results are to be expected when more detailed auxiliary data are considered. With the use of the MBFH estimator, researchers are able to improve inter-regional and temporal comparisons of subgroup-specific indicators.

Summary and outlook
Socioeconomic indicators, such as the median income of sub-populations, are key to both policy recommendations and evaluation. Taking the freely available US ACS data on county-level median income of Hispanic or Latino Americans as an 1 3 example, many county-specific direct estimates for 2010 and 2011 are associated with high standard errors or unpublished. Thereby, a comparison of different counties or time-points is partly not possible.
We present a new estimator based on small area estimation techniques and bivariate modeling, the MBFH predictor. It provides best predictors for many missing direct estimates, is easily applicable, allows for variable-specific choices of auxiliary data, and is flexible with respect to the correlation patterns and choice of the variables of interest. The MBFH predictor is a generalization of the bivariate FH model. That means in all situations where a bivariate FH model can be applied, the MBFH predictor can be applied as well, with the difference that the MBFH predictor also includes information of domains for which only one direct estimate is available. We furthermore derive an MSE estimator for the MBFH predictor and the synthetic FH predictor for observed and missing domains. Both are comparable in their relative bias. They therefore allow for a comparison of the MBFH and FH predictor for observed and missing values. The MSE estimator is convenient from a computational point of view because no resampling methods are needed to specify the uncertainty of the MBFH predictor.
A large-scale model-based Monte Carlo simulation study reveals the advantages of the MBFH predictor over the corresponding univariate FH models. For that, we consider the effects separately for the groups of domains with complete information and those where one direct estimate is missing. First, for domains with complete information the predictor can be advantageous over the FH-EBP when the variable random effects are highly correlated, positively or negatively. Thereby, it preserves the good qualities of the bivariate FH model but includes also information of domains with partially missing information in the parameter estimation. Second, for domains with one missing direct estimate we have to distinguish between the performance of the EBP for the variable with an observed and that with a missing direct estimate. For the variable with an observed direct estimate, the MBFH predictor does not give any performance gains over the corresponding FH model. On the other hand, there is also no performance loss visible in the simulation studies. For the variable with a missing direct estimate, the MBFH predictor gives significant performance gains over the corresponding FH synthetic predictor when the variable random effects are highly correlated. These are precisely the cases for which the MBFH predictor is designed and where we see the main contribution of the proposed approach.
We emphasize that we do not use imputation methods in the common way, as our goal is not to do analysis on the imputed data set. Our goal is to get best predictions 1 3 of socioeconomic indicators for the domains. Hence, our method does not stand in concurrence against imputation, as these two methods have different goals. We do not have any information on the units in the domains with missing direct estimates. Hence, the only possibility is to assume that this area follows the selected model. We do not work with missing data in the sense of non-response. We deal with unpublished direct estimates in some domains, where the sample size is typically null or too small. A basic imputation method under a selected model would be using a synthetic estimator. Instead, we propose an EBP based on the MBFH model. For an illustration, we applied the MBFH predictor to the median income of Hispanic or Latino Americans in 2010 and 2011 where publicly available direct estimates are volatile and partially missing. The validation with external data shows that the MBFH-EBP provides realistic results in practice. We therefore recommend the use of the MBFH-EBP for the estimation of regional indicators, especially in the context of unavailable direct domain estimates and unsampled domains. In the application, we detected some heavy outliers, and an investigation of a robust version of the proposed MBFH model therefore represents an interesting further area of study. In the future, we plan to publish the presented algorithm in an R package. In any case, the current R codes are available on request.