# Unifying area and unit-level small area estimation through calibration

William Acero∗
Department of Statistics and Operations Research
Complutense University of Madrid
and
Isabel Molina
Institute of Interdisciplinary Mathematics (IMI)
Department of Statistics and Operations Research
Complutense University of Madrid  This work was supported by the Ministry of
Science and Innovation of Spain PID2020-115598RB-I0.

###### Abstract

When estimating area means, direct estimators based on area-specific data, are
usually consistent under the sampling design without model assumptions.
However, they are inefficient if the area sample size is small. In small area
estimation, model assumptions linking the areas are used to “borrow strength”
from other areas. The basic area-level model provides design-consistent
estimators but error variances are assumed to be known. In practice, they are
estimated with the (scarce) area-specific data. These estimators are
inefficient, and their error is not accounted for in the associated mean
squared error estimators. Unit-level models do not require to know the error
variances but do not account for the survey design. Here we describe a unified
estimator of an area mean that may be obtained both from an area-level model
or a unit-level model and based on consistent estimators of the model error
variances as the number of areas increases. We propose bootstrap mean squared
error estimators that account for the uncertainty due to the estimation of the
error variances. We show a better performance of the new small area estimators
and our bootstrap estimators of the mean squared error. We apply the results
to education data from Colombia.

Keywords: Mean square error, Variance estimation, Calibration, Area-level
model, Unit-level model.

##  1 Introduction

When estimating linear characteristics of a finite population, survey
practitioners rely on expansion estimators, obtained as weighted sums of the
sample observations using the survey weights, which are often previously
calibrated. Design-based inference provides estimators with good properties
under the sampling replication mechanism. A good property of common design-
based estimators, such as the usual expansion estimators, is design
consistency as the sample size increases, which holds without any model
assumption for the values of the target variable. Since the survey sample size
is typically planned to achieve a desired level of precision when estimating
at the population level or for certain subpopulations, design-based estimators
perform well for those planned domains.

Often, after the survey is conducted, disaggregated estimates are required for
other subpopulations or areas, for which the survey was not originally
planned. In this case, the expansion estimators obtained with the area-
specific survey data, called direct estimates, may lack precision due to small
area sample sizes. Small area estimation (SAE) methods address the problem of
data scarcity by making model assumptions about the study variable. These
assumptions often link all the target areas through common parameters. These
common parameters are estimated with the survey data from all the areas,
leading to “indirect” estimators of the target area characteristic that
“borrow strength” from related areas. Model-based estimators typically improve
the precision of direct estimators, provided that model assumptions hold.

When estimating small area means, two well-known models are often used. The
first one is the area-level model proposed by Fay & Herriot (1979), referred
to as the FH model, which uses only aggregated data at the area level.
Specifically, this model linearly relates the direct estimators of the target
area characteristics to a vector of true area-level means of auxiliary
variables, and includes area effects and model errors. Based on this model,
empirical best linear unbiased predictors (EBLUPs) are obtained. Direct
estimates are previously calculated using the unit-level data on the target
variable from a survey, whereas the true area means of the auxiliary variables
are typically obtained from census data.

The second model is the unit-level model proposed by Battese et al. (1988),
referred to as the BHF model, which assumes a linear regression model for the
individual values of the target variable in terms of the available unit-level
covariates and includes random area effects. Linear characteristics are then
estimated with the EBLUPs obtained under this model, which also requires the
population means of the covariates for the area, in addition to the unit-level
survey data on the target variable and the covariates. Both FH and BHF models
belong to the class of linear mixed models. Their main difference is that FH
model is fit using only area-level data, whereas BHF model is fit with the
unit-level survey data.

An important advantage of FH model over BHF one is that the resulting EBLUPs
based on FH model may use the survey weights through the direct estimates,
which act as response variables in the FH model. As consequence, the EBLUPs
based on FH model are design consistent as the area sample size becomes large,
unlike EBLUPs of area means based on BHF model, which ignore the sampling
design.

You & Rao (2002) obtained a Pseudo EBLUP for a small area mean under BHF model
that is design consistent as the area sample size increases. They provided an
analytical approximation for the mean square error (MSE) of the Pseudo EBLUP,
as well as an estimator of the MSE. Moreover, if survey weights for the units
in each area are calibrated to add up to the known population count of that
area, the Pseudo EBLUPs of the area totals automatically satisfy the
convenient self-benchmarking property of adding up to the survey regression
estimator of the population total.

A well-known issue of the FH model is that the variances of the model errors
are assumed to be known. Since these variances vary across areas, they cannot
be estimated from the available area-level data. This problem is typically
solved by replacing these variances with some estimates of the sampling
variances of the direct estimators, obtained from the unit-level area-specific
survey data. These estimates are then treated as the true error variances in
the FH model. However, they are also subject to high sampling variability for
areas with small sample sizes, and this variability is typically ignored.

A common approach to avoid the problem of the high variability of the (direct)
estimators of the sampling variances is to smooth them using a previous model
for the estimated variances in terms of the same or different predictors.
However, the choice of predictors, which tunes the smoothing level of the
final variance estimates, obtained as predicted values in this preliminary
model, is often unclear and arbitrary. Note that the purpose of this
preliminary model is not to find the estimated variances that best fit the
direct variance estimates, as that would eliminate the smoothing effect.

Prasad & Rao (1990) obtained a second-order approximation for the MSE of the
EBLUP of a small area mean under the FH model, as the number of areas tends to
infinity. They also provided an estimator of the MSE with second-order bias,
as the number of areas increases. However, this estimator again assumes that
the true error variances in the FH model are known. When these true error
variances are replaced by the estimated sampling variances of direct
estimators, the resulting Prasad-Rao MSE estimator might underestimate the
true MSE. Consequently, gains of the usual EBLUPs based on FH model with
respect to direct estimators might be overstated in the areas with small
sample sizes.

Bell (2008) and Rivest & Vandal (2003) studied the effect of replacing model
error variances with estimated sampling variances of direct estimators in the
MSE of the EBLUP based on FH model. They considered simple random sampling
within the areas and assumed that the other FH model parameters were known.
Simulation results of Rivest & Vandal (2003) indicate that the Prasad-Rao MSE
estimator can underestimate the true MSE for areas with small sample sizes.
Rivest & Vandal (2003) proposed an estimator of the MSE that accounts for the
uncertainty due to the estimation of these variances. Wang & Fuller (2003)
provided another estimator of the MSE when all the model parameters are
estimated, using certain moment estimators, but further research is needed to
extend their results to other sampling designs and fitting methods. Fay (2018)
proposed to consider the survey-regression estimator as direct estimator in
the FH model, when unit-level data are available, and used a smoothed version
of its estimated variance as the true error variance. He found a better
performance of the resulting EBLUP for some simulation scenarios in contrast
to using those using customary direct estimators but, again, the smoothing
procedure is not clear. Rao et al. (2023) studied the relative efficiency of
EBLUPs obtained under Fay’s proposed model (referred to as hybrid model)
compared to the EBLUPs based on the BHF unit-level model. They obtained MSE
ratios of EBLUPs under BHF model by those under the hybrid model close to one,
showing no major gain in efficiency.

In this paper, we demonstrate how to obtain the area-level FH model through
aggregation in the unit-level BHF model, utilizing the survey weights and
calibration at the area totals of the covariates. The aggregation process
provides a correct specification for the error variances in the aggregated FH
model in terms of the variance of the unit-level model errors in the BHF
model. This implies that the FH model parameters, including the
heteroscedastic error variances, can be estimated using the unit-level data on
the target variable and the covariates from the survey. Utilizing unit-level
data leads to increased precision due to the larger number of observations
used to fit the model.

In cases where no access is provided to the unit-level survey data, the FH
model with the same specification of error variances can be also fitted with
the area-level data. In both cases, we propose bootstrap MSE estimators that
account for the uncertainty due to the estimation of error variances,
applicable to general fitting methods and sampling designs.

We compare the usual EBLUPs based on the FH model, which use estimated
sampling variances as the true error variances, with the EBLUPs obtained by
fitting the unit-level BHF model that leads to the FH model through
aggregation. Subsequently, we illustrate that the customary EBLUPs based on FH
model result in greater MSEs, especially for areas with the smallest sample
sizes. Finally, we investigate the underestimation of the MSE when using the
Prasad-Rao MSE estimator, replacing true error variances with estimated
sampling variances. We show that our bootstrap MSE estimators provide better
approximations to the true MSEs, even for very small area sample sizes.

##  2 Area-level FH model

Let $U$ be a finite population of size $N$, composed of $D$ non-overlapping
areas or domains $U_{1},\ldots,U_{D}$, each with sizes $N_{1},\ldots,N_{D}$.
Let $y_{di}$ represent the value of the study variable for the $i$th unit
within the $d$th area, $i=1,\ldots,N_{d}$, and $d=1,\ldots,D$. The objective
is to estimate the area means $\mu_{d}=N_{d}^{-1}\sum_{i=1}^{N_{d}}y_{di}$,
for $d=1,\ldots,D$.

To estimate $\mu_{d}$ for each $d$, we use data from a survey conducted on the
entire population, measuring the variable of interest. In this survey, a
sample $s_{d}\subseteq U_{d}$ of size $n_{d}>0$ with $n_{d}\leq N_{d}$, is
assumed to be drawn independently from each area $d$, for $d=1,\ldots,D$,
according to a given sampling design. The overall sample is denoted as
$s=s_{1}\cup\cdots\cup s_{D}$ with a total sample size of
$n=\sum_{d=1}^{D}n_{d}$. From the $n_{d}$ survey observations in area $d$, a
direct estimator $\hat{\mu}_{d}^{DIR}$ of $\mu_{d}$ is obtained for
$d=1,\ldots,D$. The problem arises when direct estimators
$\hat{\mu}_{d}^{DIR}$ for some areas exhibit unacceptable sampling errors, due
to small $n_{d}$.

Fay & Herriot (1979) addressed this problem by considering the model known as
FH model and defined in two stages. In the first stage, the true area mean
$\mu_{d}$ were assumed to vary linearly with respect to a vector
$\bar{\boldsymbol{X}}_{d}$ with the area means of $p<n$ covariates:

$\mu_{d}=\bar{\boldsymbol{X}}_{d}^{\prime}\boldsymbol{\beta}+u_{d},\quad
u_{d}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\sigma_{u}^{2})\quad
d=1,\ldots,D.$ (1)

Here, $\boldsymbol{\beta}$ is a $p\times 1$ vector of unknown regression
coefficients, and $u_{d}$, $d=1,\ldots,D$, are known as area effects, with
$\sigma_{u}^{2}>0$ being unknown. This model, known as the linking model,
allows to “borrow strength” by linking all the areas through the common
$\boldsymbol{\beta}$.

In model (1), the true means $\mu_{d}$, are not observable but we consider
that direct estimators $\hat{\mu}_{d}^{DIR}$ of $\mu_{d}$ are available from a
survey. These direct estimators have non-negligible sampling variances, which
vary by areas because the area sample sizes $n_{d}$ are generally different.
To account for the sampling error of the direct estimators, at the second
stage, Fay & Herriot (1979) assumed the following model, known as the sampling
model,

$\hat{\mu}_{d}^{DIR}=\mu_{d}+e_{d},\quad e_{d}\stackrel{{\scriptstyle
ind}}{{\sim}}N(0,\psi_{d}),\quad d=1,\ldots,D.$ (2)

Here, the $e_{d}$’s represent the sampling errors, assumed to be independent
of the area effects $u_{d}$, and with error variances $\psi_{d}>0$ assumed to
be known for $d=1,\ldots,D$. By combining (1) and (2), the FH model can be
expressed as a single linear mixed model,

$\hat{\mu}_{d}^{DIR}=\bar{\boldsymbol{X}}_{d}^{\prime}\boldsymbol{\beta}+u_{d}+e_{d},\quad
d=1,\ldots,D.$ (3)

Let $\boldsymbol{\psi}=(\psi_{1},\ldots,\psi_{D})^{\prime}$ be the vector of
error variances. The best predictor (BP) of the mixed effect
$\mu_{d}=\bar{\boldsymbol{X}}_{d}^{\prime}\boldsymbol{\beta}+u_{d}$ under the
FH model (3) is the predictor $\tilde{\mu}_{d}$ obtained in terms of the
available data
$\boldsymbol{y}=(\hat{\mu}_{1}^{DIR},\ldots,\hat{\mu}_{D}^{DIR})^{\prime}$
that minimizes the MSE under the FH model (3),
$\mbox{MSE}(\tilde{\mu}_{d})=E[(\tilde{\mu}_{d}-\mu_{d})^{2}]$. Under
normality of $u_{d}$ and $e_{d}$, the BP of
$\mu_{d}=\bar{\boldsymbol{X}}_{d}^{\prime}\boldsymbol{\beta}+u_{d}$ is

$\tilde{\mu}_{d}(\boldsymbol{\beta},\sigma_{u}^{2},\boldsymbol{\psi})=\gamma_{d}(\sigma_{u}^{2},\boldsymbol{\psi})\hat{\mu}_{d}^{DIR}+(1-\gamma_{d}(\sigma_{u}^{2},\boldsymbol{\psi}))\bar{\boldsymbol{X}}_{d}^{\prime}\boldsymbol{\beta},\quad\gamma_{d}(\sigma_{u}^{2},\boldsymbol{\psi})=\frac{\sigma_{u}^{2}}{\sigma_{u}^{2}+\psi_{d}}.$
(4)

According to (4), the BP of $\mu_{d}$ under the FH model is a weighted average
of the direct estimator $\hat{\mu}_{d}^{DIR}$ and the regression-synthetic
estimator $\bar{\boldsymbol{X}}_{d}^{\prime}\boldsymbol{\beta}$, tending to
$\hat{\mu}_{d}^{DIR}$ as $\psi_{d}$ goes to zero (equivalently, as $n_{d}$
grows) and to the regression-synthetic estimator otherwise.

The MSE of the BP
$\tilde{\mu}_{d}(\boldsymbol{\beta},\sigma_{u}^{2},\boldsymbol{\psi})$
obtained with the true values of $\boldsymbol{\beta}$, $\sigma_{u}^{2}$ and
$\boldsymbol{\psi}$ satisfies

$\mbox{MSE}[\tilde{\mu}_{d}(\boldsymbol{\beta},\sigma_{u}^{2},\boldsymbol{\psi})]=\gamma_{d}(\sigma_{u}^{2},\boldsymbol{\psi})\psi_{d}<\psi_{d}=\mbox{\mbox{Var}}(\hat{\mu}_{d}^{DIR}|\mu_{d});$

hence, the BP
$\tilde{\mu}_{d}(\boldsymbol{\beta},\sigma_{u}^{2},\boldsymbol{\psi})$ is more
efficient than the direct estimator $\hat{\mu}_{d}^{DIR}$, for all
$d=1,\ldots,D$. However, in practice, $\boldsymbol{\beta}$, $\sigma_{u}^{2}$
and the error variances $\psi_{d}$, $d=1,\ldots,D$, are unknown.

Assuming that $\sigma_{u}^{2}$ and $\psi_{d}$, $d=1,\ldots,D$, are known, but
avoiding normality in (1) and (2), Fay & Herriot (1979) obtained the best
linear unbiased predictor (BLUP) of the mixed effect
$\mu_{d}=\bar{\boldsymbol{X}}_{d}^{\prime}\boldsymbol{\beta}+u_{d}$, which is
the optimal predictor in terms of minimizing the MSE under the FH model (3),
among the class of predictors $\tilde{\mu}_{d}$ that are linear on
$\boldsymbol{y}=(\hat{\mu}_{1}^{DIR},\ldots,\hat{\mu}_{D}^{DIR})^{\prime}$ and
unbiased in the sense $E[\tilde{\mu}_{d}-\mu_{d}]$. The BLUP of
$\mu_{d}=\bar{\boldsymbol{X}}_{d}^{\prime}\boldsymbol{\beta}+u_{d}$ equals the
BP $\tilde{\mu}_{d}(\boldsymbol{\beta},\sigma_{u}^{2},\boldsymbol{\psi})$,
when the true $\boldsymbol{\beta}$ is replaced by the weighted least squares
(WLS) estimator obtained from area-level data
$\\{(\hat{\mu}_{d}^{DIR},\bar{\boldsymbol{X}}_{d}),\ d=1,\ldots,D\\}$ and
given by

$\tilde{\boldsymbol{\beta}}_{FH}(\sigma_{u}^{2},\boldsymbol{\psi})=\left(\sum_{d=1}^{D}\gamma_{d}(\sigma_{u}^{2},\boldsymbol{\psi})\bar{\boldsymbol{X}}_{d}\bar{\boldsymbol{X}}_{d}^{\prime}\right)^{-1}\sum_{d=1}^{D}\gamma_{d}(\sigma_{u}^{2},\boldsymbol{\psi})\bar{\boldsymbol{X}}_{d}\hat{\mu}_{d}^{DIR}.$
(5)

Therefore, for $\sigma_{u}^{2}$ and $\boldsymbol{\psi}$ known, the BLUP of
$\mu_{d}$ under the FH model (3) is given by

$\tilde{\mu}_{d}^{FH}(\sigma_{u}^{2},\boldsymbol{\psi})=\tilde{\mu}_{d}\left(\tilde{\boldsymbol{\beta}}_{FH}(\sigma_{u}^{2},\boldsymbol{\psi}),\sigma_{u}^{2},\boldsymbol{\psi}\right)=\gamma_{d}(\sigma_{u}^{2},\boldsymbol{\psi})\hat{\mu}_{d}^{DIR}+\left(1-\gamma_{d}(\sigma_{u}^{2},\boldsymbol{\psi})\right)\bar{\boldsymbol{X}}_{d}^{\prime}\tilde{\boldsymbol{\beta}}_{FH}(\sigma_{u}^{2},\boldsymbol{\psi}).$
(6)

For known $\psi_{d}$, $d=1,\ldots,D$, a consistent estimator of
$\sigma_{u}^{2}$ may be obtained by maximum likelihood (ML), restricted ML
(REML) or by a moments method, such as the one proposed by Fay & Herriot
(1979), known as FH fitting method. In terms of relative efficiency, these
three estimation methods are equivalent (see Rao & Molina, 2015, sec. 6.1.3).

Replacing a consistent estimator $\hat{\sigma}_{u}^{2}(\boldsymbol{\psi})$ of
$\sigma_{u}^{2}$ as $D\to\infty$ in $\eqref{blupFayModel}$, we obtain the
empirical BLUP (EBLUP) of $\mu_{d}$. Denoting
$\hat{\boldsymbol{\beta}}_{FH}(\boldsymbol{\psi})=\tilde{\boldsymbol{\beta}}_{FH}\left(\hat{\sigma}_{u}^{2}(\boldsymbol{\psi}),\boldsymbol{\psi}\right)$
and
$\hat{\gamma}_{d}(\boldsymbol{\psi})=\gamma_{d}(\hat{\sigma}_{u}^{2}(\boldsymbol{\psi}),\boldsymbol{\psi})$,
the EBLUP of $\mu_{d}$ based on the FH model (3) with error variances
$\psi_{d}$, $d=1,\ldots,D$, reads

$\hat{\mu}_{d}^{FH}(\boldsymbol{\psi})=\tilde{\mu}_{d}^{FH}\left(\hat{\sigma}_{u}^{2}(\boldsymbol{\psi}),\boldsymbol{\psi}\right)=\hat{\gamma}_{d}(\boldsymbol{\psi})\hat{\mu}_{d}^{DIR}+\left(1-\hat{\gamma}_{d}(\boldsymbol{\psi})\right)\bar{\boldsymbol{X}}_{d}^{\prime}\hat{\boldsymbol{\beta}}_{FH}(\boldsymbol{\psi}),\
\hat{\gamma}_{d}(\boldsymbol{\psi})=\frac{\hat{\sigma}_{u}^{2}(\boldsymbol{\psi})}{\hat{\sigma}_{u}^{2}(\boldsymbol{\psi})+\psi_{d}}.$
(7)

Since, for $\boldsymbol{\psi}$ known,
$\hat{\sigma}_{u}^{2}(\boldsymbol{\psi})$ and
$\hat{\boldsymbol{\beta}}_{FH}(\boldsymbol{\psi})$ are consistent estimators
of $\sigma_{u}^{2}$ and $\boldsymbol{\beta}$, respectively, as $D\to\infty$,
we expect $\hat{\mu}_{d}^{FH}(\boldsymbol{\psi})$ to be close to the BP
$\tilde{\mu}_{d}(\boldsymbol{\beta},\sigma_{u}^{2},\boldsymbol{\psi})$ for $D$
sufficiently large, which is more efficient than the direct estimator
$\hat{\mu}_{d}^{DIR}$.

Assuming normality of $u_{d}$ and $e_{d}$, known and uniformly bounded error
variances $\psi_{d}$, $d=1,\ldots,D$, and $\sup_{1\leq d\leq
D}h_{dd}=O(D^{-1})$, for
$h_{dd}=\bar{\boldsymbol{X}}_{d}^{\prime}\left(\sum_{d=1}^{D}\bar{\boldsymbol{X}}_{d}\bar{\boldsymbol{X}}_{d}^{\prime}\right)^{-1}\bar{\boldsymbol{X}}_{d}$,
$d=1,\ldots,D$, Prasad & Rao (1990) approximated the MSE of
$\hat{\mu}_{d}^{FH}(\boldsymbol{\psi})$ in (7) with

$\displaystyle\mbox{MSE}_{PR}[\hat{\mu}_{d}^{FH}(\boldsymbol{\psi})]$
$\displaystyle=$ $\displaystyle
g_{1d}(\sigma_{u}^{2},\boldsymbol{\psi})+g_{2d}(\sigma_{u}^{2},\boldsymbol{\psi})+g_{3d}(\sigma_{u}^{2},\boldsymbol{\psi}),$
(8)

where

$\displaystyle g_{1d}(\sigma_{u}^{2},\boldsymbol{\psi})$ $\displaystyle=$
$\displaystyle\gamma_{d}(\sigma_{u}^{2},\boldsymbol{\psi})\psi_{d},$
$\displaystyle g_{2d}(\sigma_{u}^{2},\boldsymbol{\psi})$ $\displaystyle=$
$\displaystyle\left(1-\gamma_{d}\sigma_{u}^{2},\boldsymbol{\psi})\right)^{2}\bar{\boldsymbol{X}}_{d}^{\prime}\left(\sum_{d=1}^{D}\frac{\bar{\boldsymbol{X}}_{d}\bar{\boldsymbol{X}}_{d}^{\prime}}{\psi_{d}+\sigma_{u}^{2}}\right)^{-1}\bar{\boldsymbol{X}}_{d},$
$\displaystyle g_{3d}(\sigma_{u}^{2},\boldsymbol{\psi})$ $\displaystyle=$
$\displaystyle\left(1-\gamma_{d}(\sigma_{u}^{2},\boldsymbol{\psi})\right)^{2}\frac{\overline{\mbox{var}}(\hat{\sigma}_{u}^{2})}{\sigma^{2}_{u}+\psi_{d}}.$

Here, $\overline{\mbox{var}}(\hat{\sigma}_{u}^{2})$ denotes the asymptotic
variance of $\hat{\sigma}_{u}^{2}(\boldsymbol{\psi})$, which depends on the
procedure employed to estimate $\sigma_{u}^{2}$. The term
$g_{1d}(\sigma_{u}^{2},\boldsymbol{\psi})$ accounts for the uncertainty due to
the prediction of the area effect $u_{d}$ and is $O(1)$ as $D\to\infty$,
whereas $g_{2d}(\sigma_{u}^{2},\boldsymbol{\psi})$ and
$g_{3d}(\sigma_{u}^{2},\boldsymbol{\psi})$ are the uncertainty due to the
estimation of $\boldsymbol{\beta}$ and $\sigma_{u}^{2}$, respectively, and are
both $O(D^{-1})$ as $D\to\infty$. The approximation (8) is second-order
correct, meaning that

$\mbox{MSE}[\hat{\mu}_{d}^{FH}(\boldsymbol{\psi})]=\mbox{MSE}_{PR}[\hat{\mu}_{d}^{FH}(\boldsymbol{\psi})]+o(D^{-1}).$

Based on (8), the authors provided a second-order unbiased estimator of the
MSE of $\hat{\mu}_{d}^{FH}(\boldsymbol{\psi})$, that is, with $o(D^{-1})$
bias, which again depends on the estimation method used for $\sigma_{u}^{2}$.
If $\hat{\sigma}_{u}^{2}(\boldsymbol{\psi})$ is the REML estimator of
$\sigma_{u}^{2}$, the second-order unbiased MSE estimator is

$\mbox{mse}_{PR}[\hat{\mu}_{d}^{FH}(\boldsymbol{\psi})]=g_{1d}(\hat{\sigma}_{u}^{2}(\boldsymbol{\psi}),\boldsymbol{\psi})+g_{2d}(\hat{\sigma}_{u}^{2}(\boldsymbol{\psi}),\boldsymbol{\psi})+2g_{3d}(\hat{\sigma}_{u}^{2}(\boldsymbol{\psi}),\boldsymbol{\psi}).$
(9)

In practice, the error variances in
$\boldsymbol{\psi}=(\psi_{1},\ldots,\psi_{D})^{\prime}$ are unknown. Hence, we
need to estimate them to obtain the EBLUPs based on FH model,
$\hat{\mu}_{d}^{FH}(\boldsymbol{\psi})$. A common approach is to replace
$\psi_{d}$ with the estimated design-based variance of the direct estimator,
$\psi_{d0}=\widehat{\mbox{Var}}_{\pi}(\hat{\mu}_{d}^{DIR}|\mu_{d})$, which is
obtained from the $n_{d}$ unit-level survey observations from area $d$ and
hence it is also “direct”. Defining
$\boldsymbol{\psi}_{0}=(\psi_{10},\ldots,\psi_{D0})^{\prime}$, the FH
estimator obtained by setting $\boldsymbol{\psi}=\boldsymbol{\psi}_{0}$ is
denoted here as
$\hat{\mu}_{d}^{FHD}=\hat{\mu}_{d}^{FH}(\boldsymbol{\psi}_{0})$.

Smoothed versions of $\psi_{d0}$, $d=1,\ldots,D$, obtained as fitted values
from a preliminary regression model for them are often employed to avoid the
high instability of $\psi_{d0}$. However, no clear procedure for selecting the
predictors in that model seems to be given, and this determines the level of
smoothing of the variance estimators.

Moreover, using the approximation $\mbox{MSE}_{PR}(\hat{\mu}_{d}^{FHD})$
obtained by setting $\boldsymbol{\psi}=\boldsymbol{\psi}_{0}$ in (8) is likely
to understate the true MSE of
$\hat{\mu}_{d}^{FHD}=\hat{\mu}_{d}^{FH}(\boldsymbol{\psi}_{0})$ for small
$n_{d}$, a problem that inherits the associated MSE estimator,
$\mbox{mse}_{PR}(\hat{\mu}_{d}^{FHD})$, obtained by setting
$\boldsymbol{\psi}=\boldsymbol{\psi}_{0}$ in (9). They are not accounting for
the uncertainty due to the estimation of $\psi_{d}$ by $\psi_{d0}$,
$d=1,\ldots,D$, which is not negligible for $n_{d}$ small, as will be
illustrated in the simulation experiments of Section 8. Actually, as will be
also shown in our simulations, in areas with a small sample size $n_{d}$,
since the true MSE of $\hat{\mu}_{d}^{FHD}$ is larger than expected,
$\hat{\mu}_{d}^{FHD}$ might be even less efficient than the corresponding
direct estimator $\hat{\mu}_{d}^{DIR}$.

##  3 Unit-level BHF model

Battese et al. (1988) proposed a nested error linear regression model for data
at the unit level, in the context of estimating mean acreage under a crop for
counties (small areas) in Iowa. The model, assumed for each population unit
and known as BHF model, is given by

$y_{di}=\boldsymbol{x}_{di}^{\prime}\boldsymbol{\beta}+u_{d}+e_{di},\quad
i=1,\ldots,N_{d},\ d=1,\ldots,D,$ (10)

where $\boldsymbol{x}_{di}=(x_{di1},\ldots,x_{dip})^{\prime}$ contains the
unit-specific values of $p<n$ auxiliary variables, $\boldsymbol{\beta}$ is a
$p\times 1$ vector of unknown regression coefficients, $u_{d}$ is the random
effect of area $d$, assumed to satisfy $u_{d}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,\sigma_{u}^{2})$, $d=1,\ldots,D$, and $e_{di}$ are the unit
model errors, satisfying $e_{di}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,\sigma_{e}^{2})$, $i=1,\ldots,N_{d}$, and all independent of
the area effects. If sample selection bias is absent, the sample measurements
$(y_{di},\boldsymbol{x}_{di})$, $i\in s_{d}$, $d=1,\ldots,D$, follow the same
population model given in (10).

Battese et al. (1988) obtained the EBLUP of
$\mu_{d}=\bar{\boldsymbol{X}}_{d}^{\prime}\boldsymbol{\beta}+u_{d}$ under the
model (10) and Prasad & Rao (1990) provided a second-order approximation for
the MSE of the EBLUP, together with a second-order unbiased MSE estimator.
However, the EBLUP of $\mu_{d}$ under the BHF model (10) ignores the sampling
design. As a consequence, it is not design consistent as the area sample size
$n_{d}$ grows, unlike the EBLUP based on the FH model (3).

Prasad & Rao (1999) introduced a Pseudo EBLUP for a small area mean that uses
the survey weights and is design-consistent as the area sample size $n_{d}$
increases. Let $w_{di}>0$ be the survey weight for unit $i$ in the sample from
area $d$ and $w_{d\cdot}=\sum_{i\in s_{d}}w_{di}$. The survey-weighted direct
estimator $\bar{y}_{dw}=w_{d\cdot}^{-1}\sum_{i\in s_{d}}w_{di}y_{di}$ is
design-consistent as $n_{d}$ increases. They noted that taking weighted
average in BHF model (10) over the sample units in area $d$ results in the
following aggregated area-level model,

$\bar{y}_{dw}=\bar{\boldsymbol{x}}_{dw}^{\prime}\boldsymbol{\beta}+u_{d}+\bar{e}_{dw},\quad
d=1,\ldots,D,$ (11)

where $\bar{\boldsymbol{x}}_{dw}=w_{d\cdot}^{-1}\sum_{i\in
s_{d}}w_{di}\boldsymbol{x}_{di}$ and $\bar{e}_{dw}=w_{d\cdot}^{-1}\sum_{i\in
s_{d}}w_{di}e_{di}$. Moreover, under the unit-level model (10), the model
errors in (11) satisfy $\bar{e}_{dw}\stackrel{{\scriptstyle
ind}}{{\sim}}N(0,\psi_{d}(\sigma_{e}^{2}))$, where

$\psi_{d}(\sigma_{e}^{2})=\sigma_{e}^{2}w_{d\cdot}^{-2}\sum_{i\in
s_{d}}w_{di}^{2},\quad d=1,\ldots,D.$ (12)

The aggregated model (11) has an important advantage over the usual FH model,
which is that all of the error variances $\psi_{d}(\sigma_{e}^{2})$,
$d=1,\ldots,D$, depend on a common parameter $\sigma_{e}^{2}$, for which a
consistent estimator as $D\to\infty$ may be obtained from the available area-
or unit-level data for the $D$ areas. In contrast, the usual FH model assumes
that the error variances $\psi_{d}$, $d=1,\ldots,D$, are known, and they are
often replaced with “direct” estimators
$\psi_{d0}=\widehat{\mbox{Var}}_{\pi}(\hat{\mu}_{d}^{DIR}|\mu_{d})$,
$d=1,\ldots,D$, obtained from the $n_{d}$ unit-level survey observations from
area $d$. Note that fitting the aggregated model (11) requires the same data
sources as those used to fit the usual FH model.

Prasad & Rao (1999) defined the Pseudo BLUP of
$\mu_{d}=\bar{\boldsymbol{X}}_{d}^{\prime}\boldsymbol{\beta}+u_{d}$ as the
BLUP of $\mu_{d}$ under the aggregated model (11). The Pseudo BP of
$\mu_{d}=\bar{\boldsymbol{X}}_{d}^{\prime}\boldsymbol{\beta}+u_{d}$ under
normality is

$\tilde{\mu}_{d}^{YR}(\boldsymbol{\beta},\sigma_{u}^{2},\sigma_{e}^{2})=\gamma_{d}(\sigma_{u}^{2},\sigma_{e}^{2})\left[\bar{y}_{dw}+(\bar{\boldsymbol{X}}_{d}-\bar{\boldsymbol{x}}_{dw})^{\prime}\boldsymbol{\beta}\right]+\left(1-\gamma_{d}(\sigma_{u}^{2},\sigma_{e}^{2})\right)\bar{\boldsymbol{X}}_{d}^{\prime}\boldsymbol{\beta},$
(13)

where now
$\gamma_{d}(\sigma_{u}^{2},\sigma_{e}^{2})=\sigma_{u}^{2}/(\sigma_{u}^{2}+\psi_{d}(\sigma_{e}^{2}))$,
for $\psi_{d}(\sigma_{e}^{2})$ given in (12).

Prasad & Rao (1999) proposed to estimate $\boldsymbol{\beta}$ in the
aggregated area-level model (11) using the aggregated data
$\\{(\bar{y}_{dw},\bar{\boldsymbol{x}}_{dw});\ d=1,\ldots,D\\}$. However,
$\sigma_{u}^{2}$ and $\sigma_{e}^{2}$ were estimated using a moments method
based on the unit-level data. Then Prasad & Rao (1999) defined the Pseudo
EBLUP of $\mu_{d}$ by replacing the unknown $\boldsymbol{\beta}$,
$\sigma_{u}^{2}$ and $\sigma_{e}^{2}$, with the resulting estimates in (13).
Later, You & Rao (2002) noted that, since $n\geq D$, a more efficient Pseudo
EBLUP is obtained by estimating $\boldsymbol{\beta}$ in the original unit-
level model (10) using the unit-level data
$\\{(y_{di},\boldsymbol{x}_{di},w_{di});\ i\in s_{d},\,d=1,\ldots,D\\}$ as
follows,

$\tilde{\boldsymbol{\beta}}_{U}(\sigma_{u}^{2},\sigma_{e}^{2})=\left[\sum_{d=1}^{D}\sum_{i\in
s_{d}}w_{di}\boldsymbol{x}_{di}\left(\boldsymbol{x}_{di}-\gamma_{d}(\sigma_{u}^{2},\sigma_{e}^{2})\bar{\boldsymbol{x}}_{dw}\right)^{\prime}\right]^{-1}\sum_{d=1}^{D}\sum_{i\in
s_{d}}w_{di}\left(\boldsymbol{x}_{di}-\gamma_{d}(\sigma_{u}^{2},\sigma_{e}^{2})\bar{\boldsymbol{x}}_{dw}\right)y_{di}.$
(14)

Estimators $\hat{\sigma}_{u,U}^{2}$ and $\hat{\sigma}_{e,U}^{2}$ of
$\sigma_{u}^{2}$ and $\sigma_{e}^{2}$, respectively, were also obtained using
the method of Henderson III (H3) based on the unit-level data. You & Rao
(2002) then defined another Pseudo EBLUP of $\mu_{d}$ as
$\hat{\mu}_{d}^{YR}=\tilde{\mu}_{d}^{YR}(\hat{\boldsymbol{\beta}}_{U},\hat{\sigma}_{u,U}^{2},\hat{\sigma}_{e,U}^{2})$,
where
$\hat{\boldsymbol{\beta}}_{U}=\tilde{\boldsymbol{\beta}}_{U}(\hat{\sigma}_{u,U}^{2},\hat{\sigma}_{u,U}^{2})$.

You & Rao (2002) noted that, if survey weights $w_{di}$, $i\in s_{d}$, are
calibrated so that $\sum_{i\in s_{d}}w_{di}=N_{d}$, the Pseudo EBLUP
$\hat{\mu}_{d}^{YR}$ obtained considering intercept (i.e., with $x_{di1}=1$)
automatically satisfies a benchmarking property. More specifically, when
adding the Pseudo EBLUPs of the area totals over all the areas of the
population, we obtain the survey regression estimator of the population total
$Y=\sum_{d=1}^{D}\sum_{i=1}^{N_{d}}y_{di}$, that is,

$\sum_{d=1}^{D}N_{d}\hat{\mu}_{d}^{YR}=\hat{Y}_{w}+(\boldsymbol{X}-\hat{\boldsymbol{X}}_{w})^{\prime}\hat{\boldsymbol{\beta}}_{U},$
(15)

where $\hat{Y}_{w}=\sum_{d=1}^{D}\sum_{i\in s_{d}}w_{di}y_{di}$ and
$\hat{\boldsymbol{X}}_{w}=\sum_{i\in s_{d}}w_{di}\boldsymbol{x}_{di}$ are the
survey-weighted direct estimators of the population totals $Y$ and
$\boldsymbol{X}=\sum_{d=1}^{D}\sum_{i=1}^{N_{d}}\boldsymbol{x}_{di}$,
respectively.

You & Rao (2002) approximated the MSE of the Pseudo EBLUP $\hat{\mu}_{d}^{YR}$
when $\sigma_{u}^{2}$ and $\sigma_{e}^{2}$ are estimated by H3 method, with

$\mbox{MSE}_{YR}(\hat{\mu}_{d}^{YR})=g_{1dw}(\sigma_{u}^{2},\sigma_{e}^{2})+g_{2dw}(\sigma_{u}^{2},\sigma_{e}^{2})+g_{3dw}(\sigma_{u}^{2},\sigma_{e}^{2}),$
(16)

where

$g_{1dw}(\sigma_{u}^{2},\sigma_{e}^{2})=\left(1-\gamma_{d}(\sigma_{u}^{2},\sigma_{e}^{2})\right)\sigma_{u}^{2},$
$g_{2dw}(\sigma_{u}^{2},\sigma_{e}^{2})=\left(\bar{\boldsymbol{X}}_{d}-\gamma_{d}(\sigma_{u}^{2},\sigma_{e}^{2})\bar{\boldsymbol{x}}_{dw}\right)^{\prime}\Phi_{w}\left(\bar{\boldsymbol{X}}_{d}-\gamma_{d}(\sigma_{u}^{2},\sigma_{e}^{2})\bar{\boldsymbol{x}}_{dw}\right),$
$g_{3dw}(\sigma_{u}^{2},\sigma_{e}^{2})=\gamma_{d}(\sigma_{u}^{2},\sigma_{e}^{2})\left(1-\gamma_{d}(\sigma_{u}^{2},\sigma_{e}^{2})\right)^{2}\sigma_{e}^{-4}\sigma_{u}^{-2}h(\sigma_{u}^{2},\sigma_{e}^{2}),$
(17)

where $h(\sigma_{u}^{2},\sigma_{e}^{2})$ and $\Phi_{w}$ and are respectively
given in (7) and (16) of You & Rao (2002). For the same fitting method, they
provided the following second-order unbiased MSE estimator,

$\mbox{mse}_{YR}(\hat{\mu}_{d}^{YR})=g_{1dw}(\hat{\sigma}_{u,U}^{2},\hat{\sigma}_{e,U}^{2})+g_{2dw}(\hat{\sigma}_{u,U}^{2},\hat{\sigma}_{e,U}^{2})+2g_{3dw}(\hat{\sigma}_{u,U}^{2},\hat{\sigma}_{e,U}^{2}).$
(18)

##  4 Unification of unit-level and area-level models

This section shows that aggregation of BHF unit-level model using the survey
weights leads to a FH model, provided that the survey weights are calibrated
to the known population totals of the auxiliary variables at the area-level.

Calibration estimators were defined by Deville & Särndal (1992) as expansion
estimators based on new weights $w_{di}^{C}$, $i\in s_{d}$, obtained by
minimizing $\sum_{i\in s_{d}}d(w_{di},w_{di}^{C})$, where $d(\cdot,\cdot)$ is
a pseudo-distance, subject to certain calibration restrictions. In this
section, we consider that $w_{di}^{C}$, $i\in s_{d}$, are obtained by
calibration of the original survey weights $w_{di}$, $i\in s_{d}$, to the
known population totals of the auxiliary variables
$\boldsymbol{X}_{d}=\sum_{i=1}^{N_{d}}\boldsymbol{x}_{di}$ in area $d$, that
is, satisfying the calibration restrictions $\sum_{i\in
s_{d}}w_{di}^{C}\boldsymbol{x}_{di}=\boldsymbol{X}_{d}$. Taking the first
element of $\boldsymbol{x}_{di}$ as one, the first calibration restriction is
actually $\sum_{i\in s_{d}}w_{di}^{C}=N_{d}$.

The most common type of calibration is linear calibration, obtained using the
chi-squared pseudo-distance
$d(w_{di},w_{di}^{C})=(w_{di}-w_{di}^{C})^{2}/w_{di}$. Actually, linear
calibration of the area survey weights $w_{di}$, $i\in s_{d}$, to the area
totals $\boldsymbol{X}_{d}$ leads to the well-known GREG estimator of the
population total for area $d$, $\hat{Y}_{d}^{GREG}=\sum_{i\in
s_{d}}w_{di}^{C}y_{di}$.

Let $\bar{\boldsymbol{x}}_{dw}^{C}=N_{d}^{-1}\sum_{i\in
s_{d}}w_{di}^{C}\boldsymbol{x}_{di}$ and
$\bar{y}_{dw}^{C}=N_{d}^{-1}\sum_{i\in s_{d}}w_{di}^{C}y_{di}$ be the weighted
area means of the response variable and the auxiliary variables, respectively,
using the calibrated weights. Note that the calibration restrictions are
$\bar{\boldsymbol{x}}_{dw}^{C}=\bar{\boldsymbol{X}}_{d}$. Consequently, the
model (11), obtained by aggregation of the unit-level BHF model (10), but
using now the calibrated survey weights $w_{di}^{C}$, is a FH model as in (3),
based on the direct estimator $\hat{\mu}_{d}^{DIR}=\bar{y}_{dw}^{C}$. More
specifically, the aggregated model (11) based on the calibrated weights
$w_{di}^{C}$ reads

$\bar{y}_{dw}^{C}=\bar{\boldsymbol{X}}_{d}^{\prime}\boldsymbol{\beta}+u_{d}+\bar{e}_{dw}^{C},\quad\bar{e}_{dw}^{C}\stackrel{{\scriptstyle
ind}}{{\sim}}N(0,\psi_{d}^{C}(\sigma_{e}^{2})),\quad d=1,\ldots,D.$ (19)

with error variances now given by

$\psi_{d}^{C}(\sigma_{e}^{2})=\sigma_{e}^{2}N_{d}^{-2}\sum_{i\in
s_{d}}(w_{di}^{C})^{2},\quad d=1,\ldots,D.$ (20)

We have seen that survey-weighted aggregation of the unit-level BHF model
using the calibrated weights $w_{di}^{C}$, $i\in s_{d}$, leads to a FH model.
However, as already said, the aggregated FH model (19), obtained assuming that
the unit-level BHF model holds, has an advantage over the usual FH model in
(3). In the usual FH model, error variances $\psi_{d}$, $d=1,\ldots,D$ are
assumed to be known but, in practice, they need to be estimated and are often
replaced with the unstable estimators
$\psi_{d0}=\widehat{\mbox{Var}}_{\pi}(\hat{\mu}_{d}^{DIR}|\mu_{d})$ ,
$d=1,\ldots,D$. Moreover, the uncertainty due to the estimation of $\psi_{d}$
is not accounted for in the usual MSE estimators of $\hat{\mu}_{d}^{FHD}$.

In contrast, in the aggregated area-level FH model (19) obtained with the
calibrated weights $w_{di}^{C}$, error variances
$\psi_{d}^{C}(\sigma_{e}^{2})$, $d=1,\ldots,D$, given in (20), are specified
in terms of a single parameter $\sigma_{e}^{2}$. This parameter
$\sigma_{e}^{2}$ may be estimated, similarly as in the Pseudo EBLUP, either
with the unit-level survey data by fitting the original BHF model, or by
fitting the aggregated FH model (19) to the corresponding area aggregates,
when the analyst has no access to the original unit-level survey data.

It turns out that, when using the calibrated weights $w_{di}^{C}$, $i\in
s_{d}$, the Pseudo BP
$\tilde{\mu}_{d}^{YR}(\boldsymbol{\beta},\sigma_{u}^{2},\sigma_{e}^{2})$ of
$\mu_{d}$ given in (13) based on the unit-level BHF model, equals the BP
$\tilde{\mu}_{d}(\boldsymbol{\beta},\sigma_{u}^{2},\boldsymbol{\psi})$ of
$\mu_{d}$ given in (4), based on the FH model (19). This predictor, obtained
by unified area-level and unit-level models and referred to as “unified”
predictor, is given by

$\tilde{\mu}_{d}^{U}(\boldsymbol{\beta},\sigma_{u}^{2},\sigma_{e}^{2})=\gamma_{d}^{C}(\sigma_{u}^{2},\sigma_{e}^{2})\bar{y}_{dw}^{C}+\left(1-\gamma_{d}^{C}(\sigma_{u}^{2},\sigma_{e}^{2})\right)\bar{\boldsymbol{X}}_{d}^{\prime}\boldsymbol{\beta},\quad\gamma_{d}^{C}(\sigma_{u}^{2},\sigma_{e}^{2})=\frac{\sigma_{u}^{2}}{\sigma_{u}^{2}+\psi_{d}^{C}(\sigma_{e}^{2})}.$
(21)

It holds that $\gamma_{d}^{C}(\sigma_{u}^{2},\sigma_{e}^{2})\to 1$ as
$n_{d}\to\infty$, provided that $w_{di}^{C}>0$, for all $i\in s_{d}$, and
$\max_{i\in s_{d}}(w_{di}^{C}/N_{d})=O(n_{d}^{-1})$, as proven in the
Appendix. Hence, in that case, under the same conditions under which
$\bar{y}_{dw}^{C}$ is design-consistent as $n_{d}\to\infty$, for known
$\boldsymbol{\beta}$, $\sigma_{u}^{2}$ and $\sigma_{e}^{2}$, the unified
predictor
$\tilde{\mu}_{d}^{U}(\boldsymbol{\beta},\sigma_{u}^{2},\sigma_{e}^{2})$ given
in (21) is also design consistent as $n_{d}\to\infty$.

Empirical versions of the unified predictor in (21) may be obtained, similarly
as the Pseudo EBLUP, either by fitting the original BHF model using the
available unit-level survey data, or by fitting the aggregated FH model (19)
using the aggregated area-level data. In the first case, we can use the
estimator of $\boldsymbol{\beta}$ given in (14), but using the calibrated
weights $w_{di}^{C}$ rather than the original weights $w_{di}$, denoted
$\tilde{\boldsymbol{\beta}}_{w}^{C}(\sigma_{u}^{2},\sigma_{e}^{2})$.
Consistent estimators $\hat{\sigma}_{u,U}^{2}$ and $\hat{\sigma}_{e,U}^{2}$ of
$\sigma_{u}^{2}$ and $\sigma_{e}^{2}$, respectively, as $D\to\infty$, may be
obtained by ML, REML or H3 methods, applied to the original BHF model. We
define an empirical unified (EU) predictor as
$\hat{\mu}_{d}^{U}=\tilde{\mu}_{d}^{U}(\hat{\boldsymbol{\beta}}_{U}^{C},\hat{\sigma}_{u,U}^{2},\hat{\sigma}_{e,U}^{2})$,
for
$\hat{\boldsymbol{\beta}}_{U}^{C}=\tilde{\boldsymbol{\beta}}_{U}^{C}(\hat{\sigma}_{u,U}^{2},\hat{\sigma}_{e,U}^{2})$.

As proven in the Appendix, the EU predictors of the area totals add up to the
expansion estimator of the population total $Y$ based on the calibrated
weights, $\hat{Y}_{w}^{C}=\sum_{d=1}^{D}\sum_{i\in s_{d}}w_{di}^{C}y_{di}$,
that is,

$\sum_{d=1}^{D}N_{d}\hat{\mu}_{d}^{U}=\hat{Y}_{w}^{C}.$ (22)

Now consider that the analyst is not given access to the unit-level survey
data, only to the aggregated area-level data
$\\{(\bar{y}_{dw}^{C},\bar{\boldsymbol{X}}_{d},W_{2d});\ d=1,\ldots,D\\}$,
where $W_{2d}=\sum_{i\in s_{d}}(w_{di}^{C})^{2}$. The regression coefficients
$\boldsymbol{\beta}$ under the aggregated FH model (19) may be estimated with
the usual WLS estimator, as follows

$\tilde{\boldsymbol{\beta}}_{A}^{C}(\sigma_{u}^{2},\sigma_{e}^{2})=\left(\sum_{d=1}^{D}\gamma_{d}^{C}(\sigma_{u}^{2},\sigma_{e}^{2})\bar{\boldsymbol{X}}_{d}\bar{\boldsymbol{X}}_{d}^{\prime}\right)^{-1}\sum_{d=1}^{D}\gamma_{d}^{C}(\sigma_{u}^{2},\sigma_{e}^{2})\bar{\boldsymbol{X}}_{d}\bar{y}_{dw}^{C},$
(23)

where the subscript $A$ now means that only area-level data are used.
Consistent estimators $\hat{\sigma}_{u,A}^{2}$ and $\hat{\sigma}_{e,A}^{2}$ of
$\sigma_{u}^{2}$ and $\sigma_{e}^{2}$ respectively, as $D\to\infty$, may be
also obtained based on the same area-level data, either by ML, REML or a
method of moments. Then, a second EU predictor of $\mu_{d}$ is obtained by
replacing the model parameters in the unified predictor with these estimates,
that is, as
$\hat{\mu}_{d}^{UA}=\tilde{\mu}_{d}^{U}(\hat{\boldsymbol{\beta}}_{A}^{C},\hat{\sigma}_{u,A}^{2},\hat{\sigma}_{e,A}^{2})$,
for
$\hat{\boldsymbol{\beta}}_{A}^{C}=\tilde{\boldsymbol{\beta}}_{A}^{C}(\hat{\sigma}_{u,A}^{2},\hat{\sigma}_{e,A}^{2})$.

Advantages of the Pseudo EBLUP $\hat{\mu}_{d}^{YR}$ and the two EU predictors
$\hat{\mu}_{d}^{U}$ or $\hat{\mu}_{d}^{UA}$ over the usual EBLUPs based on the
FH model is that the specification of error variances in terms of a single
parameter $\sigma_{e}^{2}$ allows us to find consistent estimators
$\hat{\psi}_{d}^{C}=\hat{\sigma}_{e}^{2}N_{d}^{-2}\sum_{i\in
s_{d}}(w_{di}^{C})^{2}$, as $D\to\infty$, for the model error variances
$\psi_{d}^{C}(\sigma_{e}^{2})$, $d=1,\ldots,D$. This, in turn, allows us to
propose consistent parametric bootstrap MSE estimators that incorporate the
uncertainty due to the estimation of the error variances, for general fitting
procedures providing consistent estimators of the model parameters under the
unit-level BHF model (10).

##  5 FH model with new specification of error variances

The agreement of the unit-level BHF and the area-level FH models holds only
when considering direct estimators $\hat{\mu}_{d}^{DIR}=\bar{y}_{dw}^{C}$ that
use the calibrated survey weights $w_{di}^{C}$, $i\in s_{d}$. On the other
hand, the FH model is not tied to this particular direct estimator. However,
the only theoretical justification of the FH model assumptions that we find is
actually through aggregation of the unit-level model. Still, calibration at
the area level is not compulsory and one may consider the direct estimator
$\hat{\mu}_{d}^{DIR}=\bar{y}_{dw}$ based on the original weights $w_{di}$,
$i\in s_{d}$, which might be previously calibrated at a larger aggregation
level.

We can still assume that the unit-level BHF model (10) holds. Then, according
to the aggregated model (11), the variance of
$\hat{\mu}_{d}^{DIR}=\bar{y}_{dw}$ given
$\mu_{d}=\bar{\boldsymbol{X}}_{d}^{\prime}\boldsymbol{\beta}+u_{d}$, is given
by

$\mbox{Var}(\hat{\mu}_{d}^{DIR}|\mu_{d})=\sigma_{e}^{2}w_{d\cdot}^{-2}\sum_{i\in
s_{d}}w_{di}^{2}\triangleq\psi_{d}(\sigma_{e}^{2}),\quad d=1,\ldots,D.$ (24)

In this section, we consider the original FH model (3) using
$\hat{\mu}_{d}^{DIR}=\bar{y}_{dw}$, but with error variances
$\psi_{d}=\psi_{d}(\sigma_{e}^{2})$ as specified in (24), depending on a
single parameter $\sigma_{e}^{2}$. This parameter can be estimated using the
area-level data from the $D$ areas. This approach accounts for the sampling
design through $\psi_{d}(\sigma_{e}^{2})$, but avoids assuming that $\psi_{d}$
are known and the application of smoothing procedures for $\psi_{d0}$,
$d=1,\ldots,D$.

Using the specification of error variances in (24) in the FH model, the BP of
$\mu_{d}=\bar{\boldsymbol{X}}_{d}^{\prime}\boldsymbol{\beta}+u_{d}$ under that
model with normality is the same as in (4), but we now write it making
explicit the dependence on a single variance parameter $\sigma_{e}^{2}$, as

$\tilde{\mu}_{d}(\boldsymbol{\beta},\sigma_{u}^{2},\sigma_{e}^{2})=\gamma_{d}(\sigma_{u}^{2},\sigma_{e}^{2})\bar{y}_{dw}+\left(1-\gamma_{d}(\sigma_{u}^{2},\sigma_{e}^{2})\right)\bar{\boldsymbol{X}}_{d}^{\prime}\boldsymbol{\beta}.$
(25)

Considering only the aggregated data, the BLUP of
$\mu_{d}=\bar{\boldsymbol{X}}_{d}^{\prime}\boldsymbol{\beta}+u_{d}$ is
obtained by replacing the true $\boldsymbol{\beta}$ in the BP
$\tilde{\mu}_{d}(\boldsymbol{\beta},\sigma_{u}^{2},\sigma_{e}^{2})$ by the WLS
estimator as in (5), given now by

$\tilde{\boldsymbol{\beta}}_{FHA}(\sigma_{u}^{2},\sigma_{e}^{2})=\left(\sum_{d=1}^{D}\gamma_{d}(\sigma_{u}^{2},\sigma_{e}^{2})\bar{\boldsymbol{X}}_{d}\bar{\boldsymbol{X}}_{d}^{\prime}\right)^{-1}\sum_{d=1}^{D}\gamma_{d}(\sigma_{u}^{2},\sigma_{e}^{2})\bar{\boldsymbol{X}}_{d}\bar{y}_{dw}.$
(26)

Consistent estimators $\hat{\sigma}_{u,A}^{2}$ and $\hat{\sigma}_{e,A}^{2}$ of
$\sigma_{u}^{2}$ and $\sigma_{e}^{2}$, respectively, as $D\to\infty$, may be
obtained by ML, REML or H3 methods applied to the FH model (3) with
$\psi_{d}=\psi_{d}(\sigma_{e}^{2})$ given in (24). Based on the resulting
estimators, we define an empirical best predictor as
$\hat{\mu}_{d}^{FHA}=\tilde{\mu}_{d}(\hat{\boldsymbol{\beta}}_{FHA},\hat{\sigma}_{u,A}^{2},\hat{\sigma}_{e,A}^{2})$,
where
$\hat{\boldsymbol{\beta}}_{FHA}=\tilde{\boldsymbol{\beta}}_{FHA}(\hat{\sigma}_{u,A}^{2},\hat{\sigma}_{e,A}^{2})$.

An asymptotic approximation to the MSE of $\hat{\mu}_{d}^{FHA}$ and its
corresponding analytical estimator may be obtained following similar arguments
as those used by Prasad & Rao (1990). In Section 7 we propose a parametric
bootstrap MSE estimator.

##  6 MSE of PEBLUP and EU predictor

When using H3 method to estimate the variance components, the MSE of PEBLUP
can be estimated by the analytical second-order unbiased estimator given in
(18). Similarly, using the same estimation method, an estimator for the MSE of
the EU predictor $\hat{\mu}_{d}^{U}$ is obtained by replacing
$\gamma_{d}(\sigma_{u}^{2},\sigma_{e}^{2})$ with
$\gamma_{d}^{C}(\sigma_{u}^{2},\sigma_{e}^{2})$ and
$\bar{\boldsymbol{x}}_{dw}$ with $\bar{\boldsymbol{x}}_{dw}^{C}$ in (18). This
estimator does account for the uncertainty due to the estimation of all the
unknown model parameters, but is specific for H3 estimators of
$\sigma_{u}^{2}$ and $\sigma_{e}^{2}$.

Below, we describe a parametric bootstrap procedure that provides alternative
estimators for the MSE of $\hat{\mu}_{d}^{YR}$ or $\hat{\mu}_{d}^{U}$,
applicable to more general fitting methods. The resulting bootstrap MSE
estimator is consistent for the true MSE as $D\to\infty$, at the same rate and
under the same conditions as the consistency of the estimators of the
parameters in the unit-level BHF model (10). This result can be proven using
the method of imitation, similar to the approach followed by González-Manteiga
et al. (2008) for H3. The bootstrap procedure is outlined below for
$\hat{\mu}_{d}^{YR}$. Replacing in this procedure $\hat{\mu}_{d}^{YR}$ with
$\hat{\mu}_{d}^{U}$, based on the expansion estimator
$\hat{\mu}_{d}^{DIR}=\bar{y}_{dw}^{C}$ that uses the calibrated survey weights
$w_{di}^{C}$ rather than the original weights $w_{di}$, it yields a parametric
bootstrap MSE estimator for $\hat{\mu}_{d}^{U}$.
Parametric bootstrap for MSE of PEBLUP and EU predictor

  1. 1.

Fit BHF model (10) to the unit-level survey data
$\\{(y_{di},\boldsymbol{x}_{di});\ i\in s_{d},\,d=1,\ldots,D\\}$, obtaining
consistent estimators $\hat{\boldsymbol{\beta}}_{U}$, $\hat{\sigma}_{u,U}^{2}$
and $\hat{\sigma}_{e,U}^{2}$ of $\boldsymbol{\beta}$, $\sigma_{u}^{2}$ and
$\sigma_{e}^{2}$, respectively.

  2. 2.

For $b=1,\ldots,B$ with $B$ large:

    1. (i)

Generate area effects and bootstrap errors for sampled units as

$u_{d}^{*(b)}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,\hat{\sigma}_{u,U}^{2}),\quad
e_{di}^{*(b)}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,\hat{\sigma}_{e,U}^{2}),\quad i\in s_{d},\ d=1,\ldots,D.$

    2. (ii)

With the generated area effects $u_{d}^{*(b)}$ from (i) and the fitted model
parameters obtained in step 1, obtain the true bootstrap area means as

$\mu_{d}^{*(b)}=\bar{\boldsymbol{X}}_{d}^{\prime}\hat{\boldsymbol{\beta}}_{U}+u_{d}^{*(b)},\quad
d=1,\ldots,D.$

    3. (iii)

Calculate the bootstrap observations for the sample units as

$y_{di}^{*(b)}=\boldsymbol{x}_{di}^{\prime}\hat{\boldsymbol{\beta}}+u_{d}^{*(b)}+e_{di}^{*(b)},\quad
i\in s_{d},\ d=1,\ldots,D.$ (27)

    4. (iv)

Fit the bootstrap unit-level BHF model (27), obtaining bootstrap model
parameter estimates $\hat{\sigma}_{u,U}^{2*(b)}$, $\hat{\sigma}_{e,U}^{2*(b)}$
and $\hat{\boldsymbol{\beta}}_{U}^{*(b)}$, and from them obtain
$\hat{\mu}_{d}^{YR*(b)}$.

  3. 3.

A parametric bootstrap estimator of $\mbox{MSE}(\hat{\mu}_{d}^{YR})$ is then
given by

$\mbox{mse}_{PB}(\hat{\mu}_{d}^{YR})=\frac{1}{B}\sum_{b=1}^{B}\left(\hat{\mu}_{d}^{YR*(b)}-\mu_{d}^{*(b)}\right)^{2}.$
(28)

For the customary EBLUP under the FH model, $\hat{\mu}_{d}^{FHD}$, obtained
with the direct estimator $\hat{\mu}_{d}^{DIR}=\bar{y}_{dw}^{C}$ based on the
calibrated weights $w_{di}^{C}$, we propose in Section 7 a consistent
parametric bootstrap MSE estimator which, unlike
$\mbox{mse}_{PR}(\hat{\mu}_{d}^{FHD})$, incorporates the uncertainty due to
the estimation of $\psi_{d}$ with
$\psi_{d0}=\widehat{\mbox{Var}}_{\pi}(\bar{y}_{dw}^{C}|\mu_{d})$,
$d=1,\ldots,D$.

##  7 MSE of EBLUP under FH model

The use of the direct estimator $\hat{\mu}_{d}^{DIR}=\bar{y}_{dw}^{C}$ based
on the calibrated survey weights $w_{di}^{C}$ in FH model provides a
specification of the error variances in terms of a common parameter
$\sigma_{e}^{2}$ for all the areas. This specification allows us to estimate
these variances consistently as $D\to\infty$, similarly to the other model
parameters. The consistency of the parameter estimators used to generate new
data under a bootstrap procedure yields consistency of the parametric
bootstrap MSE estimator.

We propose here two different parametric bootstrap estimators for the MSE of
the EBLUP $\hat{\mu}_{d}^{FHD}$ of $\mu_{d}$ based on the usual FH model (3).
The second MSE estimator uses the Prasad-Rao analytical MSE estimator given in
(9), correcting this formula for the missing uncertainty based on parametric
bootstrap.
Parametric bootstrap for MSE of EBLUP under FH with calibrated weights

  1. 1.

Obtain consistent estimators, as $D\to\infty$, $\hat{\boldsymbol{\beta}}$,
$\hat{\sigma}_{u}^{2}$ and $\hat{\sigma}_{e}^{2}$ of $\boldsymbol{\beta}$,
$\sigma_{u}^{2}$ and $\sigma_{e}^{2}$ respectively. They can be obtained
either by fitting the unit-level BHF model (10) to the unit-level survey data
$\\{(y_{di},\boldsymbol{x}_{di});\ i\in s_{d},\,d=1,\ldots,D\\}$, or by
fitting the aggregated model (19) using the area-level data
$\\{(\bar{y}_{dw}^{C},\bar{\boldsymbol{X}}_{d});\ d=1,\ldots,D\\}$.

  2. 2.

For $b=1,\ldots,B$ with $B$ large:

    1. (i)

Taking as true error variances
$\hat{\psi}_{d}=\hat{\sigma}_{e}^{2}N_{d}^{-2}\sum_{i\in
s_{d}}(w_{di}^{C})^{2}$, $d=1,\ldots,D$, generate, independently, bootstrap
area effects and sampling errors as

$u_{d}^{*(b)}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,\hat{\sigma}_{u}^{2}),\quad
e_{d}^{*(b)}\stackrel{{\scriptstyle ind}}{{\sim}}N(0,\hat{\psi}_{d}),\quad
d=1,\ldots,D.$

    2. (ii)

With $u_{d}^{*(b)}$ generated in (i), obtain the true bootstrap area means as

$\mu_{d}^{*(b)}=\bar{\boldsymbol{X}}_{d}^{\prime}\hat{\boldsymbol{\beta}}+u_{d}^{*(b)},\quad
d=1,\ldots,D.$

    3. (iii)

Compute the bootstrap direct estimators

$\hat{\mu}_{d}^{DIR*(b)}=\mu_{d}^{*(b)}+e_{d}^{*(b)},d=1,\ldots,D.$

    4. (iv)

Taking $\hat{\mu}_{d}^{DIR}=\bar{y}_{dw}^{C}$ and
$\psi_{d0}=\widehat{\mbox{\mbox{Var}}}_{\pi}(\bar{y}_{dw}^{C}|\mu_{d})$, fit
the FH model (3) to the area-level data
$\\{(\bar{y}_{dw}^{C},\bar{\boldsymbol{X}}_{d},\psi_{d0});\ d=1,\ldots,D\\}$,
and calculate $\hat{\mu}_{d}^{FHD*(b)}$.

  3. 3.

A parametric bootstrap estimator of $\mbox{MSE}(\hat{\mu}_{d}^{FHD})$ is then

$\mbox{mse}_{PB1}(\hat{\mu}_{d}^{FHD})=\frac{1}{B}\sum_{b=1}^{B}\left(\hat{\mu}_{d}^{FHD*(b)}-\mu_{d}^{*(b)}\right)^{2}.$
(29)

We now propose an alternative MSE estimator based on correcting through a
bootstrap procedure the Prasad-Rao analytical MSE estimator
$\mbox{mse}_{PR}(\hat{\mu}_{d}^{FHD})$ given in (9), by incorporating the
uncertainty due to replacing $\psi_{d}$ with
$\psi_{d0}=\widehat{\mbox{Var}}_{\pi}(\bar{y}_{dw}^{C}|\mu_{d})$. The
bootstrap procedure is the same as the previous one, but additionally, in step
2(iv), calculates the EBLUP based on the “true” bootstrap error variances,
$\hat{\boldsymbol{\psi}}=(\hat{\psi}_{1},\ldots,\hat{\psi}_{D})^{\prime}$,
$d=1,\ldots,D$, denoted as
$\hat{\mu}_{d}^{FHT*(b)}=\hat{\mu}_{d}^{FH*(b)}(\hat{\boldsymbol{\psi}})$.
Then, in step 3, it additionally calculates the parametric bootstrap estimator
of the MSE of the EBLUP under FH model based on the true error variances,
$\hat{\mu}_{d}^{FHT}=\hat{\mu}_{d}^{FH}(\boldsymbol{\psi})$, as follows

$\mbox{mse}_{PB}(\hat{\mu}_{d}^{FHT})=\frac{1}{B}\sum_{b=1}^{B}\left(\hat{\mu}_{d}^{FHT*(b)}-\mu_{d}^{*(b)}\right)^{2}.$
(30)

The MSE estimator (30) does not account for the uncertainty due to the
estimation of the error variances. Therefore, the difference between
$\mbox{mse}_{PB1}(\hat{\mu}_{d}^{FHD})$ and
$\mbox{mse}_{PB}(\hat{\mu}_{d}^{FHT})$ represents the uncertainty that we need
to add to the PR estimator $\mbox{mse}_{PR}(\hat{\mu}_{d}^{FHD})$. Defining
$\Delta_{d}=\max\left\\{0,\mbox{mse}_{PB1}(\hat{\mu}_{d}^{FHD})-\mbox{mse}_{PB}(\hat{\mu}_{d}^{FHT})\right\\}$,
another estimator of $\mbox{MSE}(\hat{\mu}_{d}^{FHD})$ is then

$\mbox{mse}_{PB2}(\hat{\mu}_{d}^{FHD})=\mbox{mse}_{PR}(\hat{\mu}_{d}^{FHD})+\Delta_{d}.$
(31)

When there is no access to the unit-level survey data, we may consider also
the EU predictor with model parameters estimated based on the area-level data,
$\hat{\mu}_{d}^{UA}=\tilde{\mu}_{d}^{U}(\hat{\boldsymbol{\beta}}_{A}^{C},\hat{\sigma}_{u,A}^{2},\hat{\sigma}_{e,A}^{2})$.
A consistent MSE estimator $\mbox{mse}_{PB1}(\hat{\mu}_{d}^{UA})$ may be
obtained with the same parametric bootstrap procedure, but employing the area-
level data in step 1 to obtain the model parameter estimators
$\hat{\boldsymbol{\beta}}_{A}^{C}$, $\hat{\sigma}_{u,A}^{2}$ and
$\hat{\sigma}_{e,A}^{2}$. Then, in step 2(iv), we calculate
$\hat{\mu}_{d}^{UA*(b)}$ with estimated error variances
$\hat{\psi}_{d,A}^{C}=\hat{\sigma}_{e,A}^{2}N_{d}^{-2}\sum_{i\in
s_{d}}(w_{di}^{C})^{2}$, $d=1,\ldots,D$, instead of $\hat{\mu}_{d}^{FHD*(b)}$.

Similarly, if the direct estimator $\hat{\mu}_{d}^{DIR}=\bar{y}_{dw}$ based on
the original weights is used, an estimator of the MSE of
$\hat{\mu}_{d}^{FHA}$, denoted as $\mbox{mse}(\hat{\mu}_{d}^{FHA})$, could be
obtained using the same bootstrap procedure as above, by replacing
$\hat{\mu}_{d}^{FHD}$ with $\hat{\mu}_{d}^{FHA}$ in step 2 (iv).

##  8 Simulation experiments

This section describes Monte Carlo (MC) simulation experiments designed to
compare the properties (bias, MSE) of the alternative estimators of the area
means $\mu_{d}$ considered here and to study the performance of the proposed
MSE estimators.

The first experiment compares the estimators of $\mu_{d}$ obtained using the
calibrated weights: 1) Direct estimator
$\hat{\mu}_{d}^{DIR}=\bar{y}_{dw}^{C}$; 2) EBLUP based on the FH model,
$\hat{\mu}_{d}^{FHD}$, obtained setting $\hat{\mu}_{d}^{DIR}=\bar{y}_{dw}^{C}$
and $\psi_{d}=\widehat{\mbox{Var}}_{\pi}(\bar{y}_{dw}^{C}|\mu_{d})$; 3)
Empirical unified predictor, $\hat{\mu}_{d}^{UA}$, where error variances are
given by $\psi_{d}^{C}(\sigma_{e}^{2})=\sigma_{e}^{2}N_{d}^{-2}\sum_{i\in
s_{d}}(w_{di}^{C})^{2}$ and $\sigma_{e}^{2}$ is estimated, along with
$\sigma_{u}^{2}$ and $\boldsymbol{\beta}$, with the area-level data; 4)
Empirical unified predictor, $\hat{\mu}_{d}^{U}$, where error variances are
also given by $\psi_{d}^{C}(\sigma_{e}^{2})=\sigma_{e}^{2}N_{d}^{-2}\sum_{i\in
s_{d}}(w_{di}^{C})^{2}$, but estimating $\sigma_{e}^{2}$, $\sigma_{u}^{2}$ and
$\boldsymbol{\beta}$ with the survey data at the unit level.

In this simulation experiment, we consider a population of $N=250\,000$ units,
distributed into $D=25$ areas, with $N_{d}=10\,000$ units in each area
$d=1\ldots,D$. We consider two continuous auxiliary variables, $x_{q}$,
$q=1,2$, whose values are generated as $x_{q,di}\sim\mbox{Gamma}(k_{qd},1)$
with $k_{1d}=5+3d/D$, $k_{2d}=2$, $i=1\ldots,N_{d}$, $d=1\ldots,D$. The vector
of true regression coefficients is taken as
$\boldsymbol{\beta}=(4,0\mbox{.}5,-0\mbox{.}4)^{\prime}$ and the variances of
the area effects and errors are taken, respectively, as
$\sigma_{u}^{2}=0\mbox{.}1^{2}$ and $\sigma_{e}^{2}=0\mbox{.}3^{2}$.

The sample $s$ is drawn by simple random sampling without replacement
(SRSWOR), independently from each area $d$, with area sample sizes $n_{d}$ are
taken from $\\{3,5,10,15,50\\}$, where each different value is repeated for 5
consecutive areas. Then, the original survey weights $w_{di}=N_{d}/n_{d}$,
$i\in s_{d}$, are calibrated for each area $d$ to
$\boldsymbol{X}_{d}=(N_{d},\bar{X}_{1d},\bar{X}_{2d})$, $d=1,\ldots,D$, where
$\bar{X}_{qd}=\sum_{i=1}^{N_{d}}x_{q,di}$, $q=1,2$, using linear calibration.

The simulation experiment proceeds as follows. For each MC simulation
$\ell=1,\ldots,L$, with $L=1\,000$, we generate a population of
$y_{di}^{(\ell)}$ values, $i=1,\ldots,N_{d}$, $d=1,\ldots,D$, from the unit-
level BHF model in (10) with the above sizes, auxiliary variables and model
parameters. From the generated population values, we calculate the true area
means $\mu_{d}^{(\ell)}=\sum_{i=1}^{N_{d}}y_{di}^{(\ell)}$, $d=1,\ldots,D$.
Then, taking the values $y_{di}^{(\ell)}$ for the sample units, $i\in s_{d}$,
$d=1,\ldots,D$, we compute the four estimators of interest, namely
$\hat{\mu}_{d}^{DIR(\ell)}$, $\hat{\mu}_{d}^{FHD(\ell)}$,
$\hat{\mu}_{d}^{UA(\ell)}$ and $\hat{\mu}_{d}^{U(\ell)}$. The sample units,
$s$, the auxiliary information
$\boldsymbol{x}_{di}=(1,x_{1,di},x_{2,di})^{\prime}$, $i=1,\ldots,N_{d}$,
$d=1,\ldots,D$, as well as the calibrated weights, are held fixed across MC
simulations.

For a generic estimator $\hat{\mu}_{d}$, we evaluate its performance in terms
of relative bias (RB) and relative root MSE (RRMSE), approximated empirically
as

$\mbox{RB}(\hat{\mu}_{d})=\frac{L^{-1}\sum_{\ell=1}^{A}(\hat{\mu}_{d}^{(\ell)}-\mu_{d}^{(\ell)})}{L^{-1}\sum_{\ell=1}^{L}\mu_{d}^{(\ell)}},\quad\mbox{RRMSE}(\hat{\mu}_{d})=\frac{\sqrt{L^{-1}\sum_{\ell=1}^{L}(\hat{\mu}_{d}^{(\ell)}-\mu_{d}^{(\ell)})^{2}}}{L^{-1}\sum_{\ell=1}^{L}\mu_{d}^{(\ell)}}.$

Average across domains of absolute RB (ARB) and RRMSE are also calculated as

$\overline{\mbox{ARB}}=D^{-1}\sum_{d=1}^{D}|\mbox{RB}(\hat{\mu}_{d})|,\quad\overline{\mbox{RRMSE}}=D^{-1}\sum_{d=1}^{D}\mbox{RRMSE}(\hat{\mu}_{d}).$

Figure 1 shows the percent RB of the four estimators, $\hat{\mu}_{d}^{DIR}$
(labelled DIR), $\hat{\mu}_{d}^{FHD}$ (labelled FHD), $\hat{\mu}_{d}^{UA}$
(labelled UA) and $\hat{\mu}_{d}^{U}$ (labelled U), for each domain
$d=1,\ldots,D$ in the $x$-axis, with area sample sizes $n_{d}$ within
parenthesis. This plot shows that the two EU predictors have similar RB. The
same occurs for DIR and FHD, which have a larger absolute RB than the EU
predictors, for small area sample sizes ($n_{d}<5$).

Figure 2 displays the corresponding percent RRMSEs. This figure shows again a
similar performance for the two EU predictors, although the one fitted using
the more detailed unit-level data has slightly smaller RRMSEs. We can see that
the RRMSE of FHD grows considerably for $n_{d}\leq 10$ and does not show much
improvement compared with the direct estimator. Actually, it becomes as
inefficient as the direct estimator for $n_{d}\leq 3$. This occurs because the
variance estimators used in FHD are “direct”, based on really small area
sample sizes $n_{d}$. When using these inefficient variance estimators, the
EBLUPs based on FH model become also inefficient. In contrast, the two EU
predictors use consistent estimators of the error variances as $D\to\infty$.
As expected, increasing the area sample size yields similar performance of all
the estimators.

Figure 1: % RB of FHD, UA and U estimators for each area. Sample sizes in
parenthesis. Figure 2: % RRMSE of FHD, UA and U estimators for each area.
Sample sizes in parenthesis.

Table 1 reports averages across the areas with the same sample size of ARB and
RRMSE, for the three model-based estimators. In this table, FHD has an average
RRMSE that is about twice that of UA for $n_{d}=3$, but this difference
reduces as $n_{d}$ increases. U has the smallest average RRMSE for $n_{d}<50$,
but for $n_{d}=50$, the average RRMSE and ARB values are very similar for the
three estimators.

Table 1: Averages across areas with the same sample size of % ARB and RRMSE, for FHD, UA and U estimators of $\mu_{d}$. | $n_{d}=3$ | $n_{d}=5$ | $n_{d}=10$ | $n_{d}=15$ | $n_{d}=50$
---|---|---|---|---|---
Estimator | $\overline{\mbox{ARB}}$ | $\overline{\mbox{RRMSE}}$ | $\overline{\mbox{ARB}}$ | $\overline{\mbox{RRMSE}}$ | $\overline{\mbox{ARB}}$ | $\overline{\mbox{RRMSE}}$ | $\overline{\mbox{ARB}}$ | $\overline{\mbox{RRMSE}}$ | $\overline{\mbox{ARB}}$ | $\overline{\mbox{RRMSE}}$
$\hat{\mu}_{d}^{FHD}$ | 0.11 | 5.43 | 0.06 | 3.21 | 0.06 | 1.67 | 0.02 | 1.13 | 0.02 | 0.59
$\hat{\mu}_{d}^{UA}$ | 0.05 | 2.53 | 0.04 | 2.07 | 0.06 | 1.58 | 0.02 | 1.13 | 0.02 | 0.59
$\hat{\mu}_{d}^{U}$ | 0.06 | 1.93 | 0.03 | 1.71 | 0.05 | 1.36 | 0.02 | 1.06 | 0.02 | 0.60

The next MC simulation experiment evaluates the performance of the two
available MSE estimators of the EU predictor $\hat{\mu}_{d}^{U}$, for finite
number of areas ($D=25$): 1) The analytical second-order unbiased estimator
$\mbox{mse}_{YR}(\hat{\mu}_{d}^{U})$, given by (18) with
$\gamma_{d}(\sigma_{u}^{2},\sigma_{e}^{2})$ replaced by
$\gamma_{d}^{C}(\sigma_{u}^{2},\sigma_{e}^{2})$ and
$\bar{\boldsymbol{x}}_{dw}$ by $\bar{\boldsymbol{x}}_{dw}^{C}$; 2) The
parametric bootstrap MSE estimator $\mbox{mse}_{PB}(\hat{\mu}_{d}^{U})$
proposed in Section 6. In this experiment, we consider area sample sizes
$n_{d}$ in the sequence $\\{5,10,15,20,25\\}$, with each value repeated for
five consecutive areas.

The true MSEs of the EU predictors were previously approximated as described
above, but now with $L=100\,000$ MC replicates. Then, a new simulation is
performed with $L=1\,000$ MC replicates, calculating in each replicate the
analytical MSE estimator $\mbox{mse}_{YR}(\hat{\mu}_{d}^{U})$ and the
parametric bootstrap estimator $\mbox{mse}_{PB}(\hat{\mu}_{d}^{U})$ given in
(28) with $B=500$ bootstrap replicates.

Figure 3 shows the true MSEs and the MC averages of the two MSE estimators,
labelled YR and PB respectively, for each area. This plot shows that YR
analytical MSE estimator overestimates the true MSEs for all the areas, but
especially for those with $n_{d}\leq 20$. On the other hand, PB tracks very
well the true MSEs for all the areas, performing clearly better than the
analytical YR estimator, even for the areas with the largest sample sizes.

Figure 3: True MSEs of $\hat{\mu}^{U}_{d}$ and MC averages of YR and PB
estimators for each area.

Next we study the alternative estimators of the MSE of FHD. The true MSE of
FHD was again previously approximated with $L=100\,000$ replicates. Figure 4
shows the true MSEs, the means across MC simulations of the Prasad-Rao
analytical MSE estimator $\mbox{mse}_{PR}(\hat{\mu}_{d}^{FHD})$ given in (9)
(labelled PR mse) and of the two alternative parametric bootstrap MSE
estimators $\mbox{mse}_{PB1}(\hat{\mu}_{d}^{FHD})$ and
$\mbox{mse}_{PB2}(\hat{\mu}_{d}^{FHD})$ given in (29) and (31) (labelled PB1
mse and PB2 mse respectively). This figure shows that PR analytical MSE
estimator, which was not designed to account for the uncertainty due to the
estimation of the error variances, underestimates seriously the true MSE of
FHD estimator for the areas with $n_{d}\leq 10$, while the two proposed
bootstrap estimators perform practically without error. For the areas with
sample sizes greater than 10, the three MSE estimators perform similarly,
becoming nearly equal for large $n_{d}$. Actually, PR mse is the analytical
version of the parametric bootstrap MSE estimator of FHT,
$\mbox{mse}_{PB}(\hat{\mu}_{d}^{FHT})$ given in (30), which assumes that the
error variances $\psi_{d}$, $d=1,\ldots,D$ are known. None of them incorporate
the error that arises from the estimation of these variances.

Figure 4: True MSEs of $\hat{\mu}^{FHD}_{d}$ and MC averages of PR, PB1 and
PB2 estimators, for each area.

We now investigate the performance of the MSE estimators for the EU predictor
$\hat{\mu}_{d}^{UA}$ based on fitting the model with the area-level data.
Figure 5 shows the true MSEs, and the MC averages of the Prasad-Rao analytical
MSE estimator $\mbox{mse}_{PR}(\hat{\mu}_{d}^{UA})$ given in (9) (labelled PR
mse) and of the parametric bootstrap MSE estimator
$\mbox{mse}_{PB}(\hat{\mu}_{d}^{UA})$ obtained with the bootstrap procedure
described in Section 7 (labelled PB1 mse). This plot shows a better
performance of PB1 estimator of the MSE, with PR mse clearly underestimating
the true MSEs for $n_{d}\leq 20$. For $n_{d}\geq 20$, both MSE estimators are
close to the true values.

Figure 5: True MSEs of $\hat{\mu}_{d}^{UA}$ and MC averages of PR and PB1
estimators for each area.

We have obtained clear efficiency gains for the EU predictors of $\mu_{d}$,
obtained when survey weights are calibrated to the known area counts and area
totals of the auxiliary variables. Of course, this type of calibration is not
compulsory when applying the FH model. However, the FH model has the problem
of estimation of $\psi_{d}$, $d=1,\ldots,D$. Even though the equivalency
between the EBLUP based on FH model and the Pseudo EBLUP based on BHF model
holds only when using the calibrated survey weights, we also studied whether
specifying the error variances $\psi_{d}$ in FH model to those variances
obtained by aggregation of the unit-level BHF model still leads to efficiency
gains compared to the usual FHD estimator obtained with the original survey
weights. Hence, in an additional simulation experiment, we compared the
performance of the alternative estimators of $\mu_{d}$ obtained with the
original weights $w_{di}$, $i\in s_{d}$: 1) Direct estimator
$\hat{\mu}_{d}^{DIR}=\bar{y}_{dw}$; 2) EBLUP based on the FH model
$\hat{\mu}_{d}^{FHD}$ that uses the direct estimator
$\hat{\mu}_{d}^{DIR}=\bar{y}_{dw}$; 3) EBLUP based on the FH model
$\hat{\mu}_{d}^{FHA}$ with error variances
$\psi_{d}=\psi_{d}(\sigma_{e}^{2})=\sigma_{e}^{2}w_{d\cdot}^{-2}\sum_{i\in
s_{d}}w_{di}^{2}$, where $\sigma_{e}^{2}$ is estimated together with
$\sigma_{u}^{2}$ and $\boldsymbol{\beta}$ with area-level data; 4) Pseudo
EBLUP (PEBLUP) $\hat{\mu}_{d}^{YR}$.

This simulation experiment was conducted under the same setup as before, with
samples in each area drawn by SRSWOR, leading to constant survey weights
$w_{di}=N_{d}/n_{d}$ for all $i\in s_{d}$, $d=1,\ldots,D$. Here the survey-
weighted direct estimator $\hat{\mu}_{d}^{DIR}=\bar{y}_{dw}$ reduces to the
area sample mean, $\bar{y}_{d}=n_{d}^{-1}\sum_{i\in s_{d}}y_{di}$, and the
error variances resulting from aggregation in the unit-level model are given
by $\psi_{d}(\sigma_{e}^{2})=\sigma_{e}^{2}/n_{d}$, $d=1,\ldots,D$.

Figure 8 in the Appendix shows RB and RRMSE of the four estimators, labelled
as DIR, FHD, FHA and PEBLUP, respectively. This figure shows that, even using
the original survey weights $w_{di}$ and fitting the FH model with area-level
data only, FHA estimator obtained with the new specification of the error
variances performs substantially better than FHD estimator. Moreover, the loss
of efficiency with respect to PEBLUP, based on fitting the model with unit-
level data, is small. Interestingly, FHD is even less efficient than the
direct estimator in several areas.

##  9 Application

In this section, we consider data from the Colombian Saber 11 test, an annual
test conducted by the Colombian Institute for the evaluation of education to
measure the quality of formal high school education. We select the test for
Calendar B students in 2018. This test evaluates Critical Reading,
Mathematics, Social and Human Sciences, Natural Sciences and English. The data
represent a census of the population of students and, from them, official
education statistics are published. The data set includes also information
about high schools, such as whether they are public or private or their
geographical location, as well as socio-economic conditions of students.

Colombia is divided into $D=33$ administrative and political divisions, 32
departments and one capital district, hereafter called departments for
simplicity. The department population sizes $N_{d}$ (number of students taking
the test in that department) range from 268 to 87 556. For this application,
we select a SRSWOR of 1.5% from each department. The total sample size is
$n=8\,250$, where $n_{d}$ varies between 4 and 1 313.

We are interested in estimating the mean $\mu_{d}$ of the Maths scores
$y_{di}$, $i=1,\ldots,N_{d}$, for each of the $D=33$ departments, using as
auxiliary information the number of rooms of the student’s household,
$x_{1,di}$. We will compare the two alternative EU predictors
$\hat{\mu}_{d}^{U}$ and $\hat{\mu}_{d}^{UA}$, with FHD estimator,
$\hat{\mu}_{d}^{FHD}$. The MSE of $\hat{\mu}_{d}^{U}$ will be estimated with
the parametric bootstrap estimator, $\mbox{mse}_{PB}(\hat{\mu}_{d}^{U})$,
described in Section 6. For $\hat{\mu}_{d}^{FHD}$, we consider the PR
analytical MSE estimator, $\mbox{mse}_{PR}(\hat{\mu}_{d}^{FHD})$, given in
Section 7, and for $\hat{\mu}_{d}^{UA}$ we obtain the parametric bootstrap
estimator obtained following the approach of Section 7,
$\mbox{mse}_{PB1}(\hat{\mu}_{d}^{UA})$.

First of all, for each area $d=1,\ldots,33$, we calibrate the survey weights
$w_{di}$, $i\in s_{d}$, to the known department count $N_{d}$ and the total
$X_{1,d}=\sum_{i=1}^{N_{d}}x_{1,di}$, using linear calibration. Now, since
both EU predictors assume that BHF model holds for the Maths scores $y_{di}$,
we check the model assumptions. Figure 9 in the Appendix shows a histogram and
a normal Q-Q plot of unit-level residuals obtained from a REML fit of the
model, given by
$\hat{e}_{di}=y_{di}-\boldsymbol{x}_{di}^{\prime}\hat{\boldsymbol{\beta}}-\hat{u}_{d}$,
for
$\hat{u}_{d}=\hat{\gamma}_{d}(\bar{y}_{d}-\bar{\boldsymbol{x}}_{d}^{\prime}\hat{\boldsymbol{\beta}})$
with
$\hat{\gamma}_{d}=\hat{\sigma}_{u}^{2}/(\hat{\sigma}_{u}^{2}+\hat{\sigma}_{e}^{2}/n_{d})$.
This figure indicates that the normal distribution is tenable, showing only a
mild deviation in the tails.

Given that no serious departures from the unit-level BHF model are found, we
depict in Figure 6 the values of the model-based estimators
$\hat{\mu}_{d}^{FHD}$, $\hat{\mu}_{d}^{UA}$ and $\hat{\mu}_{d}^{U}$ of the
Maths department means, $\mu_{d}$, with departments sorted in the ascending
order of sample size. In this application, $\hat{\mu}_{d}^{UA}$ and
$\hat{\mu}_{d}^{U}$ take practically the same values, which differ from
$\hat{\mu}_{d}^{FHD}$.

Figure 6: FHD, UA and U estimators of the means of Mathematics $\mu_{d}$ by
department, in ascending order of sample size.

Table 2 in the Appendix reports summaries of MSE estimates for FHD, UA and U.
We can see the larger MSE estimates obtained by FHD compared with UA and U.
The maximum estimated MSE for U is almost half of that for UA, and that of UA
is also clearly smaller than for FHD. Estimated MSEs of the two EU estimators
are shown for each department in Figure 7, in the ascending order of sample
size. FHD is not displayed here because its much greater MSE estimates distort
the plot. We can see that the estimated MSEs of UA are larger than those of U
for the departments with the smaller samples sizes, but they get close as the
sample size increases.

Figure 7: MSE estimates for UA and U estimators of the mean Maths scores
$\mu_{d}$ by department, in the ascending order of sample size.

##  10 Conclusions

This paper unifies the small area estimators of area means $\mu_{d}$ based on
the most common unit-level and area-level models, by previously performing
calibration of the survey weights to the area counts and totals of the
auxiliary variables. The unified model can be fitted with data at the area or
at the unit levels, requiring exactly the same data sources as the area-level
FH model. The unified model allows us to estimate consistently the error
variances in the FH model, and this in turn allows us to construct consistent
MSE estimators that take into account the uncertainty due to the estimation of
these error variances.

Our simulation experiments show that the empirical unified predictors obtained
with consistently estimated error variances, either using area-level or unit-
level data, $\hat{\mu}_{d}^{UA}$ and $\hat{\mu}_{d}^{U}$, perform clearly
better, in terms of RB and RRMSE, than the usual FH estimator
$\hat{\mu}_{d}^{FHD}$ based on direct estimators of the error variances. Our
simulations illustrate that $\hat{\mu}_{d}^{FHD}$ might be even less efficient
than the direct estimator $\hat{\mu}_{d}^{DIR}$ for some areas.

Our bootstrap estimators of $\mbox{MSE}(\hat{\mu}_{d}^{U})$ perform clearly
better than the analytical MSE estimator that ignores the uncertainty due to
the estimation of the error variances. The latter overestimates the true MSE,
specially in areas with small sample sizes. In contrast, our MSE estimators
track well the corresponding true MSEs of the empirical unified predictors
$\hat{\mu}_{d}^{UA}$ and $\hat{\mu}_{d}^{U}$ for all the areas, including
those with very small sample sizes.

## References

  * (1)
  * Battese et al. (1988) Battese, G. E., Harter, R. M. & Fuller, W. A. (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.
  * Bell (2008) Bell, W. R. (2008), Examining sensitivity of small area inferences to uncertainty about sampling error variances, in ‘Proceedings of the American Statistical Association, Survey Research Methods Section’, Vol. 327, Citeseer, p. 334.
  * Deville & Särndal (1992) Deville, J. C. & Särndal, C.-E. (1992), ‘Calibration estimation in survey sampling’, Journal of the American Statistical Association 87, 376–382.
  * Fay (2018) Fay, R. E. (2018), Further comparisons of unit-and area-level small area estimators, in ‘Survey Research Methods Section’, p. 14.
  * Fay & Herriot (1979) Fay, R. E. & 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(366a), 269–277.
  * González-Manteiga et al. (2008) González-Manteiga, W., Lombardía, M. J., Molina, I., Morales, D. & Santamaría, L. (2008), ‘Bootstrap mean squared error of a small-area eblup’, Journal of Statistical Computation and Simulation 78(5), 443–462.
  * Prasad & Rao (1990) Prasad, N. N. & Rao, J. (1990), ‘The estimation of the mean squared error of small-area estimators’, Journal of the American statistical association 85(409), 163–171.
  * Prasad & Rao (1999) Prasad, N. & Rao, J. (1999), ‘On robust small area estimation using a simple random effects model’, Survey Methodology 25, 67–72.
  * Rao et al. (2023) Rao, J., Estevao, V., Beaumont, J. F. & Bosa, K. (2023), Relative efficiency of area level and unit level small area estimators when unit level data are available, Presented in Annual Meeting of the Statistical Society of Canada, Ottawa, Canada.
  * Rao & Molina (2015) Rao, J. N. & Molina, I. (2015), Small area estimation, John Wiley & Sons, Inc., Hoboken.
  * Rivest & Vandal (2003) Rivest, L.-P. & Vandal, N. (2003), Mean squared error estimation for small areas when the small area variances are estimated, in ‘Proceedings of the International Conference on Recent Advances in Survey Sampling’, Laboratory for Research in Statistics and Probability, Carleton University.
  * Wang & Fuller (2003) Wang, J. & Fuller, W. A. (2003), ‘The mean squared error of small area predictors constructed with estimated area variances’, Journal of the American Statistical Association 98(463), 716–723.
  * You & Rao (2002) You, Y. & Rao, J. (2002), ‘A pseudo-empirical best linear unbiased prediction approach to small area estimation using survey weights’, Canadian Journal of Statistics 30(3), 431–439.

##  Appendix A Proof of design consistency of $\hat{\mu}^{U}_{d}$

We proof that, if $w_{di}^{C}>0$ for all $i\in s_{d}$ and $\max_{i\in
s_{d}}(w_{di}^{C}/N_{d})=O(n_{d}^{-1})$, then
$\gamma_{d}^{C}(\sigma_{u}^{2},\sigma_{e}^{2})\to 1$ as $n_{d}\to\infty$. We
have that

$\displaystyle\sum\limits_{i\in
s_{d}}\left(\frac{w_{di}^{C}}{N_{d}}\right)^{2}$ $\displaystyle\leq
n_{d}\max\limits_{i\in s_{d}}\left(\frac{w_{di}^{C}}{N_{d}}\right)^{2}$
$\displaystyle=n_{d}\left\\{\max\limits_{i\in
s_{d}}\left(\frac{w_{di}^{C}}{N_{d}}\right)\right\\}^{2}$
$\displaystyle=n_{d}O(n_{d}^{-2})$ $\displaystyle=O(n_{d}^{-1}).$

As a consequence,
$\psi_{d}^{C}(\sigma_{e}^{2})=\sigma_{e}^{2}\sum\limits_{i\in
s_{d}}(w_{di}^{C}/N_{d})^{2}\to 0$ as $n_{d}\rightarrow\infty$, since
$\sigma_{e}^{2}\in(0,\infty)$ is fixed. Therefore,
$\gamma_{d}^{C}(\sigma_{u}^{2},\sigma_{e}^{2})=\sigma_{u}^{2}/(\sigma_{u}^{2}+\psi_{d}^{C}(\sigma_{e}^{2}))\to
1$ as $n_{d}\rightarrow\infty$, since $\sigma_{u}^{2}\in(0,\infty)$ is also
fixed.

##  Appendix B Proof of selfbenchmarking property of $\hat{\mu}^{U}_{d}$

For given $\sigma_{u}^{2}$ and $\sigma_{e}^{2}$,
$\tilde{\boldsymbol{\beta}}_{U}^{C}(\sigma_{u}^{2},\sigma_{e}^{2})$ is the
solution of the following equation system for $\boldsymbol{\beta}$:

$\sum_{d=1}^{D}\sum_{i\in
s_{d}}w_{di}^{C}\boldsymbol{x}_{di}[y_{di}-\boldsymbol{x}_{di}^{\prime}\boldsymbol{\beta}-\tilde{u}_{dw}^{C}(\boldsymbol{\beta},\sigma_{u}^{2},\sigma_{e}^{2})]=\mathbf{0}_{p},$
(32)

where
$\tilde{u}_{dw}^{C}(\boldsymbol{\beta},\sigma_{u}^{2},\sigma_{e}^{2})=\gamma_{d}^{C}(\sigma_{u}^{2},\sigma_{e}^{2})(\bar{y}_{dw}^{C}-(\bar{{\boldsymbol{x}}}_{dw}^{C})^{\prime}\boldsymbol{\beta})$
and $\boldsymbol{x}_{di}=(x_{di1},\ldots,x_{dip})^{\prime}$ with $x_{di1}=1$.
Taking the first equation in (32), corresponding to $x_{di1}=1$, and noting
that, by the calibration restrictions, $\sum_{d=1}^{D}\sum_{i\in
s_{d}}w_{di}^{C}\boldsymbol{x}_{di}=\sum_{d=1}^{D}\boldsymbol{X}_{d}=\boldsymbol{X}$,
it holds

$\displaystyle\sum_{d=1}^{D}\sum_{i\in
s_{d}}w_{di}^{C}y_{di}-\sum_{d=1}^{D}\sum_{i\in
s_{d}}w_{di}^{C}\boldsymbol{x}_{di}^{\prime}\tilde{\boldsymbol{\beta}}_{U}^{C}(\sigma_{u}^{2},\sigma_{e}^{2})$
$\displaystyle=\sum_{d=1}^{D}\sum_{i\in
s_{d}}w_{di}^{C}\tilde{u}_{dw}^{C}(\tilde{\boldsymbol{\beta}}_{U}^{C}(\sigma_{u}^{2},\sigma_{e}^{2}),\sigma_{u}^{2},\sigma_{e}^{2})$
$\displaystyle\Leftrightarrow\hat{Y}_{w}^{C}-\boldsymbol{X}^{\prime}\tilde{\boldsymbol{\beta}}_{U}^{C}(\sigma_{u}^{2},\sigma_{e}^{2})$
$\displaystyle=\sum_{d=1}^{D}N_{d}\tilde{u}_{dw}^{C}(\tilde{\boldsymbol{\beta}}_{U}^{C}(\sigma_{u}^{2},\sigma_{e}^{2}),\sigma_{u}^{2},\sigma_{e}^{2}).$

This is also valid for $\hat{\sigma}_{u,U}^{2}$ and $\hat{\sigma}_{e,U}^{2}$,
and denoting
$\hat{\gamma}_{d}^{C}=\gamma_{d}^{C}(\hat{\sigma}_{u,U}^{2},\hat{\sigma}_{e,U}^{2})$,
$\hat{\boldsymbol{\beta}}_{U}^{C}=\tilde{\boldsymbol{\beta}}_{U}^{C}(\hat{\sigma}_{u,U}^{2},\hat{\sigma}_{e,U}^{2})$
satisfies

$\hat{Y}_{w}^{C}-\boldsymbol{X}^{\prime}\hat{\boldsymbol{\beta}}_{U}^{C}=\sum_{d=1}^{D}N_{d}\hat{u}_{dw}^{C},$

where
$\hat{u}_{dw}^{C}=\tilde{u}_{dw}^{C}(\hat{\boldsymbol{\beta}}_{U}^{C},\hat{\sigma}_{u,U}^{2},\hat{\sigma}_{e,U}^{2})$.
Then, we have

$\displaystyle\sum_{d=1}^{D}N_{d}\hat{\mu}_{d}^{U}$
$\displaystyle=\sum_{d=1}^{N_{d}}N_{d}\left[\bar{\boldsymbol{X}}_{d}^{\prime}\hat{\boldsymbol{\beta}}_{U}^{C}+\hat{\gamma}_{d}^{C}(\bar{y}_{dw}^{C}-\boldsymbol{X}_{d}^{\prime}\hat{\boldsymbol{\beta}}_{U}^{C})\right]$
$\displaystyle=\sum_{d=1}^{D}N_{d}\bar{\boldsymbol{X}}_{d}^{\prime}\hat{\boldsymbol{\beta}}^{C}_{w}+\sum_{d=1}^{D}N_{d}\hat{u}_{dw}^{C}$
$\displaystyle=\boldsymbol{X}^{\prime}\hat{\boldsymbol{\beta}}_{U}^{C}+\hat{Y}_{w}^{C}-\boldsymbol{X}^{\prime}\hat{\boldsymbol{\beta}}_{U}^{C}$
$\displaystyle=\hat{Y}_{w}^{C}.$

##  Appendix C Additional simulation results

Figure 8: % RRMSE of DIR, FHA, FHD and U estimators with original survey
weights, for each area, with area sample sizes in parentheses.

##  Appendix D Additional application results

Figure 9: Histogram (A) and normal Q-Q plot (B) of unit-level residuals from BHF model. Table 2: Summary of MSE estimates for FHD, UA and U estimators of Maths mean score. Estimator | Min. | Q1 | Median | Mean | Q3 | Max.
---|---|---|---|---|---|---
$\mbox{mse}_{PR}(\hat{\mu}_{d}^{FHD})$ | 6.83 | 7.23 | 7.58 | 8.79 | 9.57 | 20.60
$\mbox{mse}_{PB1}(\hat{\mu}_{d}^{UA})$ | 0.08 | 0.39 | 0.56 | 2.17 | 1.58 | 15.80
$\mbox{mse}_{PB}(\hat{\mu}_{d}^{U})$ | 0.11 | 0.48 | 0.72 | 1.77 | 1.56 | 8.41


ilesRGB.icc

# UltraDark.jl: A Julia package for simulation of cosmological scalar fields

Nathan Musoke  0000-0001-9839-9256

††margin:  Preprint: Prepared for submission to Journal of Open Source
Software Software • Repository License
Authors of papers retain copyright and release the work under a Creative
Commons Attribution 4.0 International License (CC BY 4.0).

## Summary

UltraDark.jl is a Julia package for the simulation of cosmological scalar
fields. Scalar fields are proposed solutions to two of the fundamental
questions in cosmology: the nature of dark matter and the universe’s initial
conditions. Modeling their dynamics requires solving the Gross-Pitaevskii-
Poisson equations, which is analytically challenging. This makes simulations
essential to understanding the dynamics of cosmological scalar fields.
UltraDark.jl is an open, performant and user friendly option for solving these
equations numerically.

Figure 1: Growth of density perturbations in a scalar field in an expanding
background, as the universe’s scale factor growths from $a=1.0$ to $a=100.0$.
The density is measured in units of the universe’s critical density
$\rho_{\mathrm{crit}}$. This scenario is similar to that in (Musoke et al.,
2020), in which the scalar field is an inflaton fragmenting after the end of
primordial inflation.

## Statement of need

Scalar fields are ubiquitous in physics, as solutions to partial differential
equations describing the spatial variation of physical quantities. As dark
matter candidates, scalar fields including axion-like particles (ALPs) would
explain the nature of the missing 85% of the universe’s matter (Adams et al.,
2022; Ade & others, 2016; Hu et al., 2000; Matos et al., 2000). As inflaton
candidates, scalar fields are proposed to cause a phase of accelerated
expansion that sets the stage for big bang nucleosynthesis (Albrecht &
Steinhardt, 1982; Amin et al., 2012; Guth, 1981; Linde, 1982). In each case, a
scalar field $\psi(t,\mathbf{x})$ represents the density of particles as a
function of space and time.

Subject to reasonable conditions, a cosmological scalar field $\psi$ whose
particles have mass $m$ obeys the Gross-Pitaevskii equation for the scalar
field,

$i\hbar\frac{\partial\psi}{\partial
t}=-\frac{\hbar^{2}}{2ma(t)^{2}}\nabla^{2}\psi+m\Phi\psi\;,$ (1)

coupled to the Poisson equation for the gravitational potential $\Phi$

$\nabla^{2}\Phi=\frac{4\pi G}{a(t)}m{|\psi|}^{2}\;,$ (2)

where $a(t)$ is the scale factor characterising the expansion of the universe.

These equations are difficult to solve analytically – even static equilibrium
solutions do not have a closed form – and necessitate the use of computer
simulations. There are codes which solving Equation 1 and Equation 2 with
different methods and in different domains, including PyUltraLight (Edwards et
al., 2018), a code written in Chapel (Padmanabhan et al., 2020), AxioNyx
(Schwabe et al., 2020), SCALAR (Mina et al., 2020), i-SPin 2 (Jain et al.,
2023; Jain & Amin, 2023); see Zhang et al. (2018) for an overview.

UltraDark.jl solves Equation 1 and Equation 2 with a pseudo-spectral
symmetrized split-step method, in which each time step consists of four sub-
steps:

$\displaystyle\psi\to\exp\left(-i\frac{h}{2}\Phi\right)\psi$ (3)
$\displaystyle\psi\to\mathcal{F}^{-1}\left\\{\exp\left(-ih\frac{k^{2}}{2}\right)\mathcal{F}\left\\{\psi\right\\}\right\\}$
(4)
$\displaystyle\Phi=\mathcal{F}^{-1}\left\\{-\frac{4\pi}{ak^{2}}\mathcal{F}\left\\{|\psi|^{2}\right\\}\right\\}$
(5) $\displaystyle\psi\to\exp\left(-i\frac{h}{2}\Phi\right)\psi,$ (6)

where $\mathcal{F}$ is a Fourier transform, $k$ are the corresponding
frequencies, and $h$ is the time step. UltraDark.jl has adaptive time steps
which allow it to accelerate simulations while preserving numerical
convergence. This is particularly useful in an expanding universe, where the
time step is roughly $h\propto a^{2}$. Such time steps result in orders of
magnitude speedups when simulating collapse of an inflaton field in the early
universe (Musoke et al., 2020).

Figure 2: Wall time for a single time step, as a function of number of CPUs.
The points represent measured times and the lines represent theoretical
$1/\text{\\#CPU}$ scalings. The circles and solid line are for grids
constructed from Arrays and the triangles and dashed lines are for MPI-
distributed PencilArrays.

Julia (Bezanson et al., 2017) has seen increasing use in scientific computing;
see for example Eschle et al. (2023) and Roesch et al. (2021) for overviews of
its use in high energy physics and biology. The use of Julia is one of the
choices that separates UltraDark.jl from similar codes. UltraDark.jl uses
Julia’s rich parallelism capabilities. The Threads.@threads macro provides
simple parallelisation of for loops. Folds.jl enables simple parallelisation
of reduction operations (Arakaki, 2020). In a cluster environment,
PencilArrays.jl and PencilFFTs.jl enable straightforward cross-node
parallelism, a capability that is challenging to reproduce in Python (Polanco,
2021a, 2021b). The scaling of this two approaches to parallelism is
demonstrated in Figure 2.

The features described above have allowed collaborators and I to produce
results presented in 4 publications, each exploring an aspect of the small
scale structure of ultralight dark matter. We have used UltraDark.jl to
explore tidal disruption in dark matter halos comprised of self-interacting
ALPs (Glennon et al., 2022), the effects of self-interactions on dynamical
friction (Glennon et al., 2024), perform the first simulations of multi-
species ALPs with intra- and inter-species interactions (Glennon, Musoke, &
Prescod-Weinstein, 2023b) and discover a novel mechanism for vortex
stabilisation in scalar dark matter (Glennon, Mirasola, Musoke, et al.,
2023a). More sample output can be found in Glennon, Mirasola, Musoke, et al.
(2023b), Glennon, Musoke, & Prescod-Weinstein (2023a), and Glennon, Musoke,
Nadler, et al. (2023).

## Acknowledgements

I thank Arka Banerjee, Ben Chandran, Richard Easther, Mateja Gosenca, Shaun
Hotchkiss, Emily Kendall, Anthony Mirasola, Ethan Nadler, Mark Neyrinck,
Chanda Prescod-Weinstein, Yourong Frank Wang, Risa Wechsler, Luna Zagorac for
discussion and comments on various stages of the development of UltraDark.jl.
Particular thanks to Noah Glennon as the first user of UltraDark.jl.

Computations were performed on Marvin, a Cray CS500 supercomputer at UNH
supported by the NSF MRI program under grant AGS-1919310. The author wishes to
acknowledge the use of New Zealand eScience Infrastructure (NeSI) high
performance computing facilities, consulting support and/or training services
as part of this research. New Zealand’s national facilities are provided by
NeSI and funded jointly by NeSI’s collaborator institutions and through the
Ministry of Business, Innovation & Employment’s Research Infrastructure
programme. This work was performed in part at Aspen Center for Physics, which
is supported by National Science Foundation under Grant No. PHY-1607611. This
work was partially supported by a grant from the Sloan Foundation.

## References

reAdams, C. B., Aggarwal, N., Agrawal, A., Balafendiev, R., Bartram, C.,
Baryakhtar, M., Bekker, H., Belov, P., Berggren, K. K., Berlin, A., Boutan,
C., Bowring, D., Budker, D., Caldwell, A., Carenza, P., Carosi, G., Cervantes,
R., Chakrabarty, S. S., Chaudhuri, S., … Zhou, T. (2022, March). Axion dark
matter. _Snowmass 2021_. https://doi.org/10.48550/ARXIV.2203.14923

preAde, P. A. R., & others. (2016). Planck 2015 results. XIII. Cosmological
parameters. _Astron. Astrophys._ , _594_ , A13.
https://doi.org/10.1051/0004-6361/201525830

preAlbrecht, A., & Steinhardt, P. J. (1982). Cosmology for Grand Unified
Theories with Radiatively Induced Symmetry Breaking. _Phys. Rev. Lett._ , _48_
, 1220–1223. https://doi.org/10.1103/PhysRevLett.48.1220

preAmin, M. A., Easther, R., Finkel, H., Flauger, R., & Hertzberg, M. P.
(2012). Oscillons After Inflation. _Phys. Rev. Lett._ , _108_ , 241302.
https://doi.org/10.1103/PhysRevLett.108.241302

preArakaki, T. (2020). Folds: Sequential, threaded, and distributed fold
interface for julia. In _GitHub repository_. GitHub.
https://github.com/JuliaFolds/Folds.jl

preBezanson, J., Edelman, A., Karpinski, S., & Shah, V. B. (2017). Julia: A
fresh approach to numerical computing. _SIAM Review_ , _59_(1), 65–98.
https://doi.org/10.1137/141000671

preEdwards, F., Kendall, E., Hotchkiss, S., & Easther, R. (2018).
PyUltraLight: A pseudo-spectral solver for ultralight dark matter dynamics.
_JCAP_ , _10_ , 027. https://doi.org/10.1088/1475-7516/2018/10/027

preEschle, J., Gal, T., Giordano, M., Gras, P., Hegner, B., Heinrich, L.,
Acosta, U. H., Kluth, S., Ling, J., Mato, P., Mikhasenko, M., Briceño, A. M.,
Pivarski, J., Samaras-Tsakiris, K., Schulz, O., Stewart, G. A., Strube, J., &
Vassilev, V. (2023). Potential of the julia programming language for high
energy physics computing. _Computing. Comput Softw Big Sci 7, 10 (2023)_ ,
_7_(1). https://doi.org/10.1007/s41781-023-00104-x

preGlennon, N., Mirasola, A. E., Musoke, N., Neyrinck, M. C., & Prescod-
Weinstein, C. (2023a). Scalar dark matter vortex stabilization with black
holes. _Journal of Cosmology and Astroparticle Physics_ , _2023_(07), 004.
https://doi.org/10.1088/1475-7516/2023/07/004

preGlennon, N., Mirasola, A. E., Musoke, N., Neyrinck, M. C., & Prescod-
Weinstein, C. (2023b). _Supplementary animations for "Scalar dark matter
vortex stabilization with black holes"_. Zenodo.
https://doi.org/10.5281/zenodo.7675830

preGlennon, N., Musoke, N., Nadler, E. O., Prescod-Weinstein, C., & Wechsler,
R. H. (2023). _Supplementary animations for "Dynamical friction in self-
interacting ultralight dark matter"_. Zenodo.
https://doi.org/10.5281/zenodo.7927474

preGlennon, N., Musoke, N., Nadler, E. O., Prescod-Weinstein, C., & Wechsler,
R. H. (2024). Dynamical friction in self-interacting ultralight dark matter.
_Phys. Rev. D_ , _109_(6), 063501. https://doi.org/10.1103/PhysRevD.109.063501

preGlennon, N., Musoke, N., & Prescod-Weinstein, C. (2023a). _Supplementary
animations for "Simulations of multi-field ultralight axion-like dark
matter"_. Zenodo. https://doi.org/10.5281/zenodo.7675775

preGlennon, N., Musoke, N., & Prescod-Weinstein, C. (2023b). Simulations of
multi-field ultralight axion-like dark matter. _Physical Review D_ , _107_(6),
063520. https://doi.org/10.1103/physrevd.107.063520

preGlennon, N., Nadler, E. O., Musoke, N., Banerjee, A., Prescod-Weinstein,
C., & Wechsler, R. H. (2022). Tidal disruption of solitons in self-interacting
ultralight axion dark matter. _Physical Review D_ , _105_(12), 123540.
https://doi.org/10.1103/physrevd.105.123540

preGuth, A. H. (1981). The Inflationary Universe: A Possible Solution to the
Horizon and Flatness Problems. _Phys. Rev. D_ , _23_ , 347–356.
https://doi.org/10.1103/PhysRevD.23.347

preHu, W., Barkana, R., & Gruzinov, A. (2000). Cold and fuzzy dark matter.
_Phys. Rev. Lett._ , _85_ , 1158–1161.
https://doi.org/10.1103/PhysRevLett.85.1158

preJain, M., & Amin, M. A. (2023). i-SPin: An integrator for multicomponent
schrödinger-poisson systems with self-interactions. _JCAP_ , _04_ , 053.
https://doi.org/10.1088/1475-7516/2023/04/053

preJain, M., Amin, M. A., & Pu, H. (2023). i-SPin 2: An integrator for general
spin-s gross-pitaevskii systems. _Phys. Rev. E 108, 055305, 15 November 2023_
, _108_(5), 055305. https://doi.org/10.1103/physreve.108.055305

preLinde, A. D. (1982). A New Inflationary Universe Scenario: A Possible
Solution of the Horizon, Flatness, Homogeneity, Isotropy and Primordial
Monopole Problems. _Phys. Lett. B_ , _108_ , 389–393.
https://doi.org/10.1016/0370-2693(82)91219-9

preMatos, T., Guzman, F. S., & Urena-Lopez, L. A. (2000). Scalar field as dark
matter in the universe. _Class. Quant. Grav._ , _17_ , 1707–1712.
https://doi.org/10.1088/0264-9381/17/7/309

preMina, M., Mota, D. F., & Winther, H. A. (2020). SCALAR: An AMR code to
simulate axion-like dark matter models. _Astron. Astrophys._ , _641_ , A107.
https://doi.org/10.1051/0004-6361/201936272

preMusoke, N., Hotchkiss, S., & Easther, R. (2020). Lighting the dark:
Evolution of the postinflationary universe. _Phys. Rev. Lett._ , _124_(6),
061301. https://doi.org/10.1103/PhysRevLett.124.061301

prePadmanabhan, N., Ronaghan, E., Zagorac, J. L., & Easther, R. (2020).
Simulating ultralight dark matter in chapel. _2020 IEEE International Parallel
and Distributed Processing Symposium Workshops (IPDPSW)_ , 678–678.
https://doi.org/10.1109/ipdpsw50202.2020.00120

prePolanco, J. I. (2021a). _PencilArrays.jl: Distributed Julia arrays using
the MPI protocol_ (Version 0.9.9). https://doi.org/10.5281/zenodo.5148035

prePolanco, J. I. (2021b). _PencilFFTs.jl: FFTs of MPI-distributed Julia
arrays_ (Version 0.12.4). https://doi.org/10.5281/zenodo.3618781

preRoesch, E., Greener, J. G., MacLean, A. L., Nassar, H., Rackauckas, C.,
Holy, T. E., & Stumpf, M. P. H. (2021). _Julia for biologists_.
https://doi.org/10.48550/ARXIV.2109.09973

preSchwabe, B., Gosenca, M., Behrens, C., Niemeyer, J. C., & Easther, R.
(2020). Simulating mixed fuzzy and cold dark matter. _Phys. Rev. D_ ,
_102_(8), 083518. https://doi.org/10.1103/PhysRevD.102.083518

preZhang, J., Liu, H., & Chu, M.-C. (2018). Cosmological simulation for fuzzy
dark matter model. _Frontiers in Astronomy and Space Sciences_ , _5_.
https://doi.org/10.3389/fspas.2018.00048

p