# Small Area Estimation using EBLUPs under the Nested Error Regression Model

Ziyang Lyu
UNSW Data Science Hub, School of Mathematics and Statisics
University of New South Wales
A.H. Welsh
Research School of Finance, Actuarial Studies and Statistics,
Australian National University Ziyang Lyu is a postdoctoral research fellow in
UNSW Data Science Hub, School of Mathematics and Statisics,University of New
South Wales, NSW, 2033, Australia. Email<EMAIL_ADDRESS>A.H.Welsh is
E.J. Hannan Professor of Statistics in the Research School of Finance,
Actuarial Studies and Statistics, Australian National University, ACT, 2601,
Australia. Email<EMAIL_ADDRESS>
###### Abstract

Estimating characteristics of domains (referred to as small areas) within a
population from sample surveys of the population is an important problem in
survey statistics. In this paper, we consider model-based small area
estimation under the nested error regression model. We discuss the
construction of mixed model estimators (empirical best linear unbiased
predictors, EBLUPs) of small area means and the conditional linear predictors
of small area means. Under the asymptotic framework of increasing numbers of
small areas and increasing numbers of units in each area, we establish
asymptotic linearity results and central limit theorems for these estimators
which allow us to establish asymptotic equivalences between estimators,
approximate their sampling distributions, obtain simple expressions for and
construct simple estimators of their asymptotic mean squared errors, and
justify asymptotic prediction intervals. We present model-based simulations
that show that in quite small, finite samples, our mean squared error
estimator performs as well or better than the widely-used Prasad & Rao (1990)
estimator and is much simpler, so is easier to interpret. We also carry out a
design-based simulation using real data on consumer expenditure on fresh milk
products to explore the design-based properties of the mixed model estimators.
We explain and interpret some surprising simulation results through analysis
of the population and further design-based simulations. The simulations
highlight important differences between the model- and design-based properties
of mixed model estimators in small area estimation.

Keywords: increasing area size asymptotics, indirect estimator, mean squared
error estimation, mixed model estimator, model-based prediction, prediction
intervals

##  1 Introduction

Estimates of area-level characteristics of interest (such as means, totals,
and quantiles) for areas, domains or clusters within a population (all
intended to be included whenever we refer to areas) obtained from sample
survey data are widely used for resource allocation in social, education and
environmental programs, and as the basis for commercial decisions. Direct
estimates which use only data specific to an area, can have large standard
errors because of relatively small area-specific sample sizes. Small area
estimation is concerned with producing more reliable estimates with valid
measures of uncertainty for the characteristics of interest; recent reviews
have been given by for example Rao (2005, 2008), Lehtonen & Veijanen (2009),
Pfeffermann (2013), and Rao & Molina (2015).

A popular method for obtaining reliable estimates (Fay & Herriot, 1979;
Battese et al., 1988) is to introduce a mixed model for the population which
includes fixed effects (to describe either unit-level or area-level effects)
and random effects (to capture additional between area variation), fit the
model using data from multiple areas and then, use the fitted model to
construct the desired estimates. When we have unit-level data, a simple and
widely used model is the nested error regression or random intercept model
(Battese et al., 1988), and a widely used method for estimating means or
totals is to use empirical best linear unbiased predictors (EBLUPs) obtained
by minimising the (prediction) mean squared error and then estimating the
unknown quantities by maximum likelihood or restricted maximum likelihood
(REML) estimation; see for example Saei & Chambers (2003a, b), Jiang & Lahiri
(2006) and Haslett & Welsh (2019). In addition to the issue of the level at
which the data are available, there is also a subtlety about the target of
estimation. The most commonly studied characteristics of interest are small
area means (or equivalently totals) and, when we assume a mixed model, the
conditional expectations of the small area means given the random effects.
These two targets, the small area means and their conditional expectations are
different and have different EBLUPs with potentially different mean squared
errors, but are treated as interchangeable in small area estimation. They are
both random variables under the model-based framework, so technically they
need to be predicted rather than estimated. However, it is common to use both
“prediction” and “estimation” in small area estimation and we refer to their
EBLUPs as mixed model estimators that are distinguished by their different
targets (Tzavidis et al., 2010); they can also be described as composite and
synthetic estimators respectively.

The model-based variability of small area estimators is usually described by
reporting estimates of their (prediction) mean squared errors or by prediction
intervals, often based on these estimates. Estimation of (prediction) mean
squared errors for mixed model estimators is complicated, even for simple
linear mixed models like the nested error regression model: estimates of
(prediction) mean squared error for mixed model estimators based on treating
the variance parameters as known (i.e. for the BLUPs rather than the EBLUPs)
are underestimates when linear mixed models are fitted to real data; and
simple, analytic expressions for the (prediction) mean squared errors of mixed
model estimators are not available, complicating their estimation. Under
normal linear mixed models (including the nested error regression model), when
the number of areas is allowed to increase while the area sizes are held fixed
(or bounded), Kackar & Harville (1981) and Prasad & Rao (1990) used second
order Taylor expansions to obtain approximations to the (prediction) mean
squared error of the EBLUPs of the conditional expectation of the small area
means, and then constructed mean squared error estimators by replacing the
unknown quantities in these approximations by estimators. The Prasad-Rao
approximation and estimator have been extended to more general models and to
allow additional estimators of the model parameters by Datta & Lahiri (2000)
and Das et al. (2004); see also Za̧dło (2009) and Torabi & Rao (2013).
Alternatives to estimators based on analytic approximations include estimators
obtained using resampling methods. Jiang et al. (2002) proposed and
investigated cluster-level jackknife methods (unusually, treating the small
area means as the characteristics of interest), Hall & Maiti (2006a) proposed
a parametric bootstrap approach for constructing bias-corrected estimates of
the prediction mean squared error and prediction regions and Chatterjee et al.
(2008) used a different approach to construct parametric bootstrap prediction
intervals. For the considerably more complicated non-normal case, Hall & Maiti
(2006b) proposed a moment-matching, double-bootstrap procedure to estimate the
prediction mean squared error.

Many extensions of the basic approach have been developed. These include, for
example, introducing outlier robust estimators (Sinha & Rao, 2009; Tzavidis et
al., 2010), spatial models to allow correlation between small areas (Saei &
Chambers, 2005; Torabi & Jiang, 2020), and different response distributions
through using generalized linear mixed models (Saei & Chambers, 2003a) and
models that allow responses with extra zeros (Chandra & Chambers, 2016). Other
work has incorporated design-based considerations (Jiang & Lahiri, 2006; Jiang
et al., 2011).

The standard asymptotic framework for model-based small area estimation under
the nested error regression model follows Kackar & Harville (1981) and Prasad
& Rao (1990) in allowing the number of areas to increase while holding the
area sizes fixed (or bounded). This framework has several disadvantages, most
notably that the small area estimators are not consistent and their asymptotic
distributions are not known. In turn, this hinders the derivation of
approximate distributions (none are available), complicates the construction
of mean squared error estimates and means that prediction intervals based on
the estimated mean squared errors cannot be shown to achieve their nominal
level even asymptotically. To overcome these difficulties, we need both the
number of areas and the sample size in each area to increase. This appears to
contradict the “small” in “small area estimation”, but the approximations
derived within this framework perform well (the ultimate purpose of deriving
asymptotic approximations) even when some areas have quite small sample size.
In addition, many practical applications include a number of large “small
areas” and no tiny ones, so it is often a practically relevant framework.
Examples occur in clinical research (clustered trials) when we study records
on large groups (areas) of patients (units) with each group treated by a
different medical practitioner or at a different hospital, in educational
research when we look at records on college students (units) grouped within
schools (areas), and in sample surveys when we observe people or households
(units) grouped in defined clusters (areas). For example, Arora & Lahiri
(1997) gave an instance with $43$ areas ranging in size from $95$ to $633$
units, and such examples are common in poverty data (Pratesi, 2016).

We use recent increasing number of areas and increasing area size asymptotic
results obtained by Lyu & Welsh (2022) for estimating the parameters of the
nested error regression model by maximum likelihood or REML estimation and by
Lyu & Welsh (2021) for the EBLUPs for the random effects in the nested error
regression model, to study small area estimators in the same framework.
Without having to assume normality in the model, we obtain simple
approximations to the distributions of the mixed model estimators that involve
the distribution of the characteristic of interest and a normal distribution.
We further obtain strikingly simple expressions for the asymptotic mean
squared errors of the estimators which are easy to estimate and can be used in
prediction intervals which are demonstrated to have correct asymptotic
coverage. Such results are are not available under the standard asymptotic
framework. They are achieved by approximating the estimators directly and
taking the (prediction) mean squared error of the approximation, rather than
directly approximating the (prediction) mean squared error. Our results fill
the practical and theoretical gap around the most widely used mixed model
estimators in small area estimation and suggest how other similar gaps can be
treated.

We describe the nested error regression model, discuss the targets of
estimation and the mixed model estimators we consider in Section 2. We present
our increasing number of areas and increasing area size asymptotic results in
Section 3 and use (model-based) simulation to demonstrate the relevance of
these results to finite samples in Section 4. We include a design-based
simulation using real data on consumer expenditure on fresh milk products, and
then use additional design-based simulations to explore some unexpected
findings in Section 5. We conclude with a brief discussion in Section 6.

##  2 Small area estimation

Consider a population $U=\cup_{i=1}^{g}U_{i}$ of $N$ units, partitioned into
$g$ exclusive areas $U_{i}$, each containing $N_{i}$ units so that
$\sum_{i=1}^{g}N_{i}=N$. Let $y_{ij}$ be a scalar survey variable of interest
and $\mathbf{x}_{ij}$ a vector of auxiliary variables for the $j$th unit in
the $i$th area. The problem of interest is to use data from a sample of units
in $U$ to make inference about the area means
$\bar{y}_{i}=N_{i}^{-1}\sum_{j=1}^{N_{i}}y_{ij}$. We assume that the values of
the auxiliary variables are known for every unit in the population and that
the values of the survey variable are observed on a sample of units
$s=\cup_{i=1}^{g}s_{i}$, where $s_{i}\subseteq U_{i}$ of size $n_{i}\leq
N_{i}$ is the set of sample units within $U_{i}$ and $\sum_{i=1}^{g}n_{i}=n$.
We assume that the units are selected using a non-informative sampling method
such that $n_{i}>n_{L}>0$, so that some units from every area are included in
the sample. (This can be weakened provided the areas with no units are a non-
informatively selected sample of the areas.) The available data $\mathcal{D}$
consists of $y_{ij}$ for $j\in s_{i}$, $i=1,\ldots,g$ and $\mathbf{x}_{ij}$
for $j=1,\ldots,N_{i}$, $i=1,\ldots,g$.

We assume the nested error regression model (or random intercept model) for
the survey variable at the population level, so

$y_{ij}=\mu(\mathbf{x}_{ij})+\alpha_{i}+e_{ij},\qquad\text{for
$j=1,\ldots,N_{i},\,i=1,\ldots,g$,}$ (2.1)

where $\mu(\mathbf{x}_{ij})$ is the regression function (the conditional mean
of the response) given $\mathbf{x}_{ij}$, $\alpha_{i}$ is a random effect
representing a random intercept or cluster effect and $e_{ij}$ is a random
error. We assume that the $\\{\alpha_{i}\\}$ and $\\{e_{ij}\\}$ are all
mutually independent with mean zero and variances (called the variance
components) $\sigma_{\alpha}^{2}$ and $\sigma_{e}^{2}$, respectively, and
write $\boldsymbol{\theta}=[\sigma_{\alpha}^{2},\sigma_{e}^{2}]^{T}$. These
random variables do not have to be normally distributed so the responses
$y_{i}$ are not necessarily normally distributed. Lyu & Welsh (2022) showed
that, in asymptotic theory with increasing area size, we need to distinguish
between area variables (variables that vary within areas) which we place in
the $p_{w}$-vector $\mathbf{x}_{ij}^{(w)}$ and the between area variables
(variables that are constant within areas) which we place in the
$p_{b}$-vector $\mathbf{x}_{i}^{(b)}$. We then write the regression function
as

$\mu(\mathbf{x}_{ij})=\beta_{0}+\mathbf{x}_{i}^{(b)T}\boldsymbol{\beta}_{1}+\mathbf{x}_{ij}^{(w)T}\boldsymbol{\beta}_{2}=\mathbf{u}_{i}^{T}\boldsymbol{\xi}+\mathbf{x}_{ij}^{(w)T}\boldsymbol{\beta}_{2},$
(2.2)

where $\beta_{0}$ is the unknown intercept, $\boldsymbol{\beta}_{1}$ is the
unknown between area slope parameter and $\boldsymbol{\beta}_{2}$ is the
unknown within area slope parameter. The intercept is a between area variable
so it is convenient to group the intercept and the between area slope terms in
the $(p_{b}+1)-$vectors $\mathbf{u}_{i}=[1,\,\,\mathbf{x}_{i}^{(b)T}]^{T}$ and
$\boldsymbol{\xi}=[\beta_{0},\boldsymbol{\beta}_{1}^{T}]^{T}$. It is also
sometimes useful to write the regression function in terms of
$\mathbf{z}_{ij}=[1,\,\,\mathbf{x}_{i}^{(b)T},\,\,\mathbf{x}_{ij}^{(w)T}]^{T}=[\mathbf{u}_{i}^{T},\,\,\mathbf{x}_{ij}^{(w)T}]^{T}$
and
$\boldsymbol{\beta}=[\beta_{0},\boldsymbol{\beta}_{1}^{T},\boldsymbol{\beta}_{2}^{T}]^{T}=[\boldsymbol{\xi}^{T},\boldsymbol{\beta}_{2}^{T}]^{T}$.
Finally, grouping the between and within area parameters, we write the full
set of model parameters as
$\boldsymbol{\omega}=[\beta_{0},\boldsymbol{\beta}_{1}^{T},\sigma_{\alpha}^{2},\boldsymbol{\beta}_{2}^{T},\sigma_{e}^{2}]^{T}$.

The form (2.2) for the regression function includes the three cases with
either or both between and within area variables. It also gives us the option
of replacing $\mathbf{x}_{ij}^{(w)}$ by the population-area-mean-centered
within area variables $\mathbf{x}_{ij}^{(w)}-\bar{\mathbf{x}}_{i}^{(w)}$,
where
$\bar{\mathbf{x}}_{i}^{(w)}=N_{i}^{-1}\sum_{k=1}^{N_{i}}\mathbf{x}_{ik}^{(w)}$,
for $j=1,\ldots,N_{i}$, $i=1,\ldots,g$, and including the area means
$\bar{\mathbf{x}}_{i}^{(w)}$ as contextual effects with the between area
variables from $\mathbf{x}_{ij}$ in $\mathbf{x}_{i}^{(b)}$. This centering
makes the between and within area variables orthogonal and allows a very
simple interpretation of the parameters from separate between and within area
models; see Yoon & Welsh (2020) and the references therein. The maximum
likelihood and REML estimators of the parameters are asymptotically orthogonal
under the conditions of Theorem1 of Lyu & Welsh (2022), so the effect is on
the model parameters rather than the estimators; this is the reason we center
about the population area means rather than the sample area means.

The two common targets when we are interested in the small area means are the
actual small area means
$\bar{y}_{i}=\mathbf{u}_{i}^{T}\boldsymbol{\xi}+\bar{\mathbf{x}}_{i}^{(w)^{T}}\boldsymbol{\beta}_{2}+\alpha_{i}+\bar{e}_{i}$
and the conditional linear predictors of the small area means

${\eta}_{i}=\operatorname{E}(\bar{y}_{i}|\mathbf{u}_{i},\mathbf{x}_{i1}^{w},\ldots,\mathbf{x}_{iN_{i}}^{(w)},\alpha_{i})=\mathbf{u}_{i}^{T}{\boldsymbol{\xi}}+\bar{\mathbf{x}}_{i}^{(w)T}{\boldsymbol{\beta}}_{2}+\alpha_{i},$

$i=1,\ldots,g$. The two targets satisfy $\bar{y}_{i}-\eta_{i}=\bar{e}_{i}$ so
are not identical and the differences do not decrease in the standard
asymptotic framework with fixed area size; the targets have the same
expectations $\operatorname{E}(\bar{y}_{i})=\operatorname{E}(\eta_{i})$ but
different variances
$\operatorname{Var}(\bar{y}_{i})=\sigma_{\alpha}^{2}+N_{i}^{-1}\sigma_{e}^{2}$
and $\operatorname{Var}(\eta_{i})=\sigma_{\alpha}^{2}$. This has practical
consequences because unbiased predictors of one will also be unbiased
predictors of the other, but the (prediction) mean squared errors of unbiased
predictors will differ with the target. We prefer to predict $\bar{y}_{i}$
rather than $\eta_{i}$ because $\bar{y}_{i}$ is an easily interpretable finite
population parameter that is is not tied to a specific population model and
can be evaluated without any prediction error when $U_{i}$ is completely
enumerated. In contrast, $\eta_{i}$ is more difficult to interpret, is not a
finite population parameter, is tied to the specific population model we are
using and is subject to prediction error even when $U_{i}$ is completely
enumerated. The use of $\eta_{i}$ as a target for prediction in small area
estimation arguably reflects the fact that early work on prediction for mixed
models like (2.1) and (2.2) (Kackar & Harville, 1981; Prasad & Rao, 1990)
which was brought into the finite population context for small area estimation
was actually done in an infinite population framework where $\bar{y}_{i}$ was
observed (so was a trivial target)) and prediction was to other (infinite)
populations (assumed to follow the same model) rather than within a given
finite population. In finite populations, prediction within the given
population is the more relevant problem. One advantage of working under our
increasing area size framework is that the two targets are asymptotically the
same up to order $N_{i}^{-1/2}$, because

$\bar{y}_{i}-\eta_{i}=\bar{e}_{i}=O_{p}(N_{i}^{-1/2}),\qquad\mbox{ as
}N_{i}\rightarrow\infty.$

This means that the predictors are quite similar in large areas and may
explain why the distinction between the two targets has been largely ignored
in practice.

The prediction mean squared error for predicting a target random variable is
minimised by the conditional expectation of the target given the observed data
$\mathcal{D}$. Typically, the distribution of the target given $\mathcal{D}$
is derived from a model with unknown parameters (such as (2.1) and (2.2)) so
the conditional expectation depends on these unknown parameters and a feasible
predictor requires replacing the unknown parameters by estimators; in our
case, we use the maximum likelihood or REML estimators of the parameters in
the model (2.1) and (2.2).

To simplify notation, we write $j\not\in s_{i}$ to mean $j\in U_{i}\backslash
s_{i}$, $k_{i}=(N_{i}-n_{i})/N_{i}$, and use subscripts $(s)$ and $(r)$ to
denote quantities related to sampled and non-sampled units. Specifically, let
$\bar{y}_{i(s)}=n_{i}^{-1}\sum_{j\in s_{i}}y_{ij}$,
$\bar{y}_{i(r)}=(N_{i}-n_{i})^{-1}\sum_{j\not\in s_{i}}y_{ij}$,
$\bar{\mathbf{x}}_{i(s)}^{(w)}=n_{i}^{-1}\sum_{j\in
s_{i}}\mathbf{x}_{ij}^{(w)}$,
$\bar{\mathbf{x}}_{i(r)}^{(w)}=(N_{i}-n_{i})^{-1}\sum_{j\not\in
s_{i}}\mathbf{x}_{ij}^{(w)}$, $\bar{e}_{i(s)}=n_{i}^{-1}\sum_{j\in
s_{i}}e_{ij}$ and $\bar{e}_{i(r)}=(N_{i}-n_{i})^{-1}\sum_{j\not\in
s_{i}}e_{ij}$. Then, under the model (2.1) and (2.2), we can write the actual
small area means $\bar{y}_{i}$ as

$\bar{y}_{i}=(1-k_{i})\bar{y}_{i(s)}+k_{i}\bar{y}_{i(r)}=(1-k_{i})\bar{y}_{i(s)}+k_{i}(\mathbf{u}_{i}^{T}\boldsymbol{\xi}+\bar{\mathbf{x}}_{i(r)}^{(w)T}\boldsymbol{\beta}_{2}+\alpha_{i}+\bar{e}_{i(r)}).$
(2.3)

Taking the conditional expectation given the data $\mathcal{D}$ of (2.3) and
substituting the maximum likelihood or restricted maximum likelihood (REML)
estimators $\hat{\boldsymbol{\xi}}$ and $\hat{\boldsymbol{\beta}}_{2}$ for
$\boldsymbol{\xi}$ and $\boldsymbol{\beta}_{2}$, respectively, and predicting
$\operatorname{E}(\alpha_{i}|\mathcal{D})$ by the empirical best linear
unbiased predictor (EBLUP)

$\hat{\alpha}_{i}=\hat{\gamma}_{i}\\{\bar{y}_{i(s)}-\mathbf{u}_{i}^{T}\hat{\boldsymbol{\xi}}-\bar{\mathbf{x}}_{i(s)}^{(w)T}\hat{\boldsymbol{\beta}}_{2}\\},\qquad\mbox{
with
}\,\hat{\gamma}_{i}=s_{i}\hat{\sigma}_{\alpha}^{2}/(\hat{\sigma}_{e}^{2}+s_{i}\hat{\sigma}_{\alpha}^{2}),$
(2.4)

we obtain the predictor

$\begin{split}\hat{M}_{i}^{sam}&=(1-k_{i})\bar{y}_{i(s)}+k_{i}\\{\mathbf{u}_{i}^{T}\hat{\boldsymbol{\xi}}+\bar{\mathbf{x}}_{i(r)}^{(w)T}\hat{\boldsymbol{\beta}}_{2}+\hat{\alpha}_{i}\\},\end{split}$
(2.5)

which the mixed model estimator of $\bar{y}_{i}$ (Tzavidis et al., 2010); the
superscript ’sam’ shows that the target is the ‘small area mean’. It can also
be called an EBLUP or a composite estimator (Costa et al., 2003) as it
combines the sample mean $\bar{y}_{i(s)}$ with a synthetic component
$\mathbf{u}_{i}^{T}\hat{\boldsymbol{\xi}}+\bar{\mathbf{x}}_{i(r)}^{(w)T}\hat{\boldsymbol{\beta}}_{2}+\hat{\alpha}_{i}$
(Rao, 2008). The mixed model estimator of $\eta_{i}$ obtained as above is

$\begin{split}\hat{M}_{i}^{\text{clp}}&=\mathbf{u}^{T}_{i}\hat{\boldsymbol{\xi}}+\bar{\mathbf{x}}^{(w)T}_{i}\hat{\boldsymbol{\beta}}_{2}+\hat{\alpha}_{i}=(1-k_{i})\\{\mathbf{u}_{i}^{T}\hat{\boldsymbol{\xi}}+\bar{\mathbf{x}}_{i(s)}^{(w)T}\hat{\boldsymbol{\beta}}_{2}+\hat{\alpha}_{i}\\}+k_{i}\\{\mathbf{u}_{i}^{T}\hat{\boldsymbol{\xi}}+\bar{\mathbf{x}}_{i(r)}^{(w)T}\hat{\boldsymbol{\beta}}_{2}+\hat{\alpha}_{i}\\},\end{split}$
(2.6)

which is an EBLUP and a fully synthetic or indirect estimator (Prasad & Rao,
1990; Lahiri & Rao, 1995; Jiang et al., 2011); the superscript ‘clp’ shows
that the target is the ‘conditional linear predictor’.

The main difference between the mixed model estimator of $\bar{y}_{i}$ (2.5)
and the mixed model estimator of $\eta_{i}$ (2.6) estimators is that the
former uses the observed $\bar{y}_{i(s)}$ whereas the latter uses a model-
based prediction for this quantity. The difference between the estimators of
the two targets can be expressed in terms of the EBLUP for the random effect
(2.4) as

$\begin{split}\hat{M}_{i}^{\text{sam}}-\hat{M}_{i}^{\text{clp}}&=(1-k_{i})\\{\bar{y}_{i(s)}-\mathbf{u}_{i}^{T}\hat{\boldsymbol{\xi}}-\bar{\mathbf{x}}_{i(s)}^{(w)T}\hat{\boldsymbol{\beta}}_{2}-\hat{\alpha}_{i}\\}=\frac{n_{i}}{N_{i}}\left\\{\frac{\hat{\alpha}_{i}}{\hat{\gamma}_{i}}-\hat{\alpha}_{i}\right\\}=\frac{\hat{\sigma}_{e}^{2}}{\hat{\sigma}_{\alpha}^{2}}\frac{\hat{\alpha}_{i}}{N_{i}}.\end{split}$
(2.7)

This difference is often quite small, but it can be large for the areas with
extreme EBLUPs $\hat{\alpha}_{i}$, particularly if these are small areas and
the estimated within area variance is much larger than the estimated between
area variance so $\hat{\sigma}_{e}^{2}>\hat{\sigma}_{\alpha}^{2}$.
Asymptotically, the difference is $O_{p}(N_{i}^{-1})$, so the estimators are
asymptotically equivalent up to this order and hence asymptotically closer
than their respective targets; see Section 3 for details. Again, these
properties only hold when the area sizes are increasing and do not hold for
fixed area-size asymptotics.

We have emphasised that the predictor depends on the target random variable,
but of course it also depends on the data and the model. We note that both
predictors (2.5) and (2.6) can be computed from the sample data
$\mathcal{D}_{(s)}=\\{(y_{ij},\mathbf{x}_{ij}^{(w)T})^{T},j=1,\ldots
n_{i},\mathbf{u}_{i},i=1,\ldots,g\\}$ and the the population means
$\bar{\mathbf{x}}_{i}^{(w)}$, $i=1,\ldots,g$. If, instead of $\mathcal{D}$, we
only observe $\mathcal{D}_{(s)}$, the population mean
$\bar{\mathbf{x}}_{i}^{(w)}$ and hence the non-sample mean of the within area
variables $\bar{\mathbf{x}}_{i(r)}^{(w)}$ is unknown and also needs to be
predicted. If we estimate $\bar{\mathbf{x}}_{i(r)}^{(w)}$ by the simple
nonparametric estimator $\bar{\mathbf{x}}_{i(s)}^{(w)}$, we obtain the mixed
model estimator/predictor of the small area mean $\bar{y}_{i}$

$\begin{split}\hat{M}_{i}^{*sam}&=(1-k_{i})\bar{y}_{i(s)}+k_{i}\\{\mathbf{u}_{i}^{T}\hat{\boldsymbol{\xi}}+\bar{\mathbf{x}}_{i(s)}^{(w)T}\hat{\boldsymbol{\beta}}_{2}+\hat{\alpha}_{i}\\}.\end{split}$

Clearly, precision in the description of predictors (especially optimal
predictors) is needed to avoid confusion.

The asymptotic results (especially Therem 3) in Section 3 below establish
that, as $g,n_{L}\to\infty$, the mixed model estimators are asymptotically
equivalent predictors of the small area means $\bar{y}_{i}$, they are
asymptotically unbiased predictors and their asymptotic prediction mean
squared errors are
$\operatorname{MSE}_{\text{LW},i}=n_{i}^{-1}k_{i}\sigma_{e}^{2}$. We can
estimate $\operatorname{MSE}_{\text{LW},i}$ by substituting the maximum
likelihood or REML estimator for $\sigma_{e}^{2}$ to obtain

$\displaystyle\widehat{\operatorname{MSE}}_{\text{LW},i}=n_{i}^{-1}k_{i}{{\hat{\sigma}_{e}}}^{2}.$
(2.8)

We can then construct simple, asymptotic $100(1-\varepsilon)\%$ prediction
intervals for $\bar{y}_{i}$ which we denote sam-LW and clp-LW, respectively,
as

$\displaystyle[\hat{M}_{i}^{\text{sam}}-\Phi^{-1}(1-\varepsilon/2)\widehat{\operatorname{MSE}}_{\text{LW},i}^{1/2},\,\hat{M}_{i}^{\text{sam}}+\Phi^{-1}(1-\varepsilon/2)\widehat{\operatorname{MSE}}_{\text{LW},i}^{1/2}],$
(2.9)
$\displaystyle[\hat{M}_{i}^{\text{clp}}-\Phi^{-1}(1-\varepsilon/2)\widehat{\operatorname{MSE}}_{\text{LW},i}^{1/2},\,\hat{M}_{i}^{\text{clp}}+\Phi^{-1}(1-\varepsilon/2)\widehat{\operatorname{MSE}}_{\text{LW},i}^{1/2}],$
(2.10)

where $\Phi^{-1}$ is the inverse of the standard normal cumulative
distribution function. The asymptotic coverage of the intervals (2.9) and
(2.10) is guaranteed by Theorem 1, Slutsky’s Theorem and Theorem 3, which do
not require the assumption of normality in the model.

##  3 Increasing area-size asymptotic results

We assume throughout that the true model describing the actual data generating
mechanism is given by (2.1) and (2.2) with true parameter
$\dot{\boldsymbol{\omega}}=[\dot{\beta}_{0},\dot{\boldsymbol{\beta}}_{1}^{T},{\dot{\sigma}_{\alpha}}^{2},\dot{\boldsymbol{\beta}}_{2}^{T},{\dot{\sigma}_{e}}^{2}]^{T}$
and take all expectations under the true model. Let
$\hat{\boldsymbol{\omega}}$ denote the normal maximum likelihood estimator
(MLE) of $\dot{\boldsymbol{\omega}}$ obtained by maximizing the normal
likelihood based on a sample $s$. Similarly, let
$\hat{\boldsymbol{\theta}}_{R}$ be the normal REML estimator of
$\dot{\boldsymbol{\theta}}=[{\dot{\sigma}_{\alpha}}^{2},{\dot{\sigma}_{e}}^{2}]^{T}$
obtained by maximizing the normal REML criterion function and let
$\hat{\boldsymbol{\beta}}(\boldsymbol{\theta})=[\hat{\beta}_{0}(\boldsymbol{\theta}),\hat{\boldsymbol{\beta}}_{1}(\boldsymbol{\theta})^{T},\hat{\boldsymbol{\beta}}_{2}(\boldsymbol{\theta})^{T}]^{T}$
be the profile likelihood estimator of
$\dot{\boldsymbol{\beta}}=[\dot{\beta}_{0},\dot{\boldsymbol{\beta}}_{1}^{T},\dot{\boldsymbol{\beta}}_{2}^{T}]^{T}$
obtained by maximizing the normal likelihood with $\boldsymbol{\theta}$ held
fixed. We call
$\hat{\boldsymbol{\beta}}_{R}=\hat{\boldsymbol{\beta}}(\hat{\boldsymbol{\theta}}_{R})$
the normal REML estimator of $\dot{\boldsymbol{\beta}}$, and
$\hat{\boldsymbol{\omega}}_{R}=(\hat{\beta}_{R0},\hat{\boldsymbol{\beta}}_{R1}^{T},\hat{\sigma}_{R\alpha}^{2},\hat{\boldsymbol{\beta}}_{R2}^{T},\hat{\sigma}_{Re}^{2})^{T}$
the normal REML estimator of $\dot{\boldsymbol{\omega}}$.

Following Lyu & Welsh (2022), we obtain the asymptotic properties of the
maximum likelihood and REML estimators of the model parameters under the
following conditions:

Condition A

  1. 1.

The model (2.1) and (2.2) holds with true parameters
$\dot{\boldsymbol{\omega}}$ inside the parameter space $\Omega$.

  2. 2.

The number of clusters $g\to\infty$ and the minimum area sample size
$n_{L}\to\infty$.

  3. 3.

The random variables $\\{\alpha_{i}\\}$ and $\\{e_{ij}\\}$ are independent and
identically distributed and there is a $\delta>0$ such that
$\operatorname{E}|\alpha_{i}|^{4+\delta}<\infty$ and
$\operatorname{E}|e_{ij}|^{4+\delta}<\infty$ for all $j=1,\ldots,{N_{i}}$,
$i=1,\ldots,g$.

  4. 4.

The limits
$\mathbf{b}_{1}=\lim_{g\rightarrow\infty}g^{-1}\sum_{i=1}^{g}\mathbf{x}_{i}^{(b)}$,
$\mathbf{B}_{2}=\lim_{g\rightarrow\infty}g^{-1}\sum_{i=1}^{g}\mathbf{x}_{i}^{(b)}\mathbf{x}_{i}^{(b)^{T}}$
and
$\mathbf{B}_{3}=\lim_{g\to\infty}\lim_{s_{L}\to\infty}n^{-1}\sum_{i=1}^{g}\sum_{j\in
s_{i}}(\mathbf{x}_{ij}^{(w)}-\bar{\mathbf{x}}_{i}^{(w)})(\mathbf{x}_{ij}^{(w)}-\bar{\mathbf{x}}_{i}^{(w)})^{T}$
exist, and the matrices $\mathbf{B}_{2}$ and $\mathbf{B}_{3}$ are positive
definite. Further,
$\lim_{g\to\infty}g^{-1}\sum_{i=1}^{g}|\bar{\mathbf{x}}_{i}^{(w)}|^{2}<\infty$,
and there exists a $\delta>0$ such that
$\lim_{g\to\infty}g^{-1}\sum_{i=1}^{g}|\mathbf{x}_{i}^{(b)}|^{2+\delta}<\infty$
and $\lim_{g\to\infty}\lim_{s_{L}\to\infty}n^{-1}\sum_{i=1}^{g}\sum_{j\in
s_{i}}|\mathbf{x}_{ij}^{(w)}-\bar{\mathbf{x}}_{i}^{(w)}|^{2+\delta}<\infty$.

As noted in Lyu & Welsh (2022), these are mild conditions. Conditions A3 and
A4 ensure that limits needed to ensure the existence of the asymptotic
variance of the estimating function exist, and that we can establish a
Lyapounov condition and hence a central limit theorem for the estimating
function. Condition A4 ensures the matrix

$\mathbf{B}=\mbox{block
diag}[\mathbf{B}_{u}/{\dot{\sigma}_{\alpha}}^{2},1/(2{\dot{\sigma}_{\alpha}}^{4}),\mathbf{B}_{3}/{\dot{\sigma}_{e}}^{2},1/(2{\dot{\sigma}_{e}}^{4})],\qquad\mbox{
with }\mathbf{B}_{u}=\left[\begin{matrix}1&\mathbf{b}_{1}^{T}\\\
\mathbf{b}_{1}&\mathbf{B}_{2}\end{matrix}\right],$

is positive definite. For later, note that
$\mathbf{B}_{u}=\lim_{g\to\infty}g^{-1}\sum_{i=1}^{g}\mathbf{u}_{i}\mathbf{u}_{i}^{T}$.
Lyu & Welsh (2022) established the following theorem for the normal maximum
likelihood and normal REML estimators.

######  Theorem 1.

(Lyu and Welsh, 2022) Suppose Condition A holds. Then, as $g,n_{L}\to\infty$,
there are solutions $\hat{\boldsymbol{\omega}}$ to the normal maximum
likelihood estimating equations and $\hat{\boldsymbol{\omega}}_{R}$ to the
normal REML estimating equations satisfying
$|\mathbf{K}^{1/2}(\hat{\boldsymbol{\omega}}-\hat{\boldsymbol{\omega}}_{R})|=O_{p}(1)$,
where
$\mathbf{K}=\operatorname{diag}(g\boldsymbol{1}_{p_{b}+2}^{T},n\boldsymbol{1}_{p_{w}+1}^{T})$.
Moreover, for both estimators, we have
$\hat{\boldsymbol{\xi}}-\dot{\boldsymbol{\xi}}=\mathbf{B}_{u}^{-1}g^{-1}\sum_{i=1}^{g}\mathbf{u}_{i}\alpha_{i}+o_{p}(g^{-1/2})$,
$\hat{\sigma}_{\alpha}^{2}-{\dot{\sigma}_{\alpha}}^{2}=g^{-1}\sum_{i=1}^{g}(\alpha_{i}^{2}-{\dot{\sigma}_{\alpha}}^{2})+o_{p}(g^{-1/2})$,
$\hat{\boldsymbol{\beta}}_{2}-\dot{\boldsymbol{\beta}}_{2}=\mathbf{B}_{3}^{-1}n^{-1}\sum_{i=1}^{g}\sum_{j\in
s_{i}}(\mathbf{x}_{ij}^{(w)}-\bar{\mathbf{x}}_{i}^{(w)})e_{ij}+o_{p}(n^{-1/2})$
and $\hat{\sigma}_{e}^{2}-\dot{\sigma}_{e}^{2}=n^{-1}\sum_{i=1}^{g}\sum_{j\in
s_{i}}(e_{ij}^{2}-{\dot{\sigma}_{e}}^{2})+o_{p}(n^{-1/2})$, as
$g,n_{L}\rightarrow\infty$, so

$\mathbf{K}^{1/2}(\hat{\boldsymbol{\omega}}-\dot{\boldsymbol{\omega}})\xrightarrow{D}N(\boldsymbol{0},\mathbf{C}),$

where

$\begin{split}\mathbf{C}&=\left[\begin{matrix}{\dot{\sigma}_{\alpha}}^{2}/(1-\mathbf{b}_{1}^{T}\mathbf{B}_{2}^{-1}\mathbf{b}_{1})&-{\dot{\sigma}_{\alpha}}^{2}\mathbf{b}_{1}^{T}\mathbf{B}_{2}^{-1}/(1-\mathbf{b}_{1}^{T}\mathbf{B}_{2}^{-1}\mathbf{b}_{1})&\operatorname{E}\alpha_{1}^{3}&\boldsymbol{0}_{[1:p_{w}]}&0\\\
-{\dot{\sigma}_{\alpha}}^{2}\mathbf{B}_{2}^{-1}\mathbf{b}_{1}/(1-\mathbf{b}_{1}^{T}\mathbf{B}_{2}^{-1}\mathbf{b}_{1})&{\dot{\sigma}_{\alpha}}^{2}\\{\mathbf{B}_{2}^{-1}+\mathbf{B}_{2}^{-1}\mathbf{b}_{1}\mathbf{b}_{1}^{T}\mathbf{B}_{2}^{-1}/(1-\mathbf{b}_{1}^{T}\mathbf{B}_{2}^{-1}\mathbf{b}_{1})\\}&\boldsymbol{0}_{[p_{b}:1]}&\boldsymbol{0}_{[p_{b}:p_{w}]}&\boldsymbol{0}_{[p_{b}:1]}\\\
\operatorname{E}\alpha_{1}^{3}&\boldsymbol{0}_{[1:p_{b}]}&\operatorname{E}\alpha_{1}^{4}-{\dot{\sigma}_{\alpha}}^{4}&\boldsymbol{0}_{[1:p_{w}]}&0\\\
\boldsymbol{0}_{[p_{w}:1]}&\boldsymbol{0}_{[p_{w}:p_{b}]}&\boldsymbol{0}_{[p_{w}:1]}&{\dot{\sigma}_{e}}^{2}\mathbf{B}_{3}^{-1}&\boldsymbol{0}_{[p_{w}:1]}\\\
0&\boldsymbol{0}_{[1:p_{b}]}&0&0&\operatorname{E}e_{ij}^{4}-{\dot{\sigma}_{e}}^{4}\end{matrix}\right]\end{split}$

and $\boldsymbol{0}_{[p:q]}$ denotes the $p\times q$-matrix of zeros.

We use Theorem 1 to derive the asymptotic distribution of the mixed model
estimators (2.5) and (2.6). We begin by establishing asymptotic linearity
results for the estimators. Recall that $k_{i}=(N_{i}-n_{i})/N_{i}$.

######  Theorem 2.

Suppose Condition A holds. Then, as $g,n_{L}\to\infty$, we have

$\begin{split}&\hat{M}_{i}^{\text{sam}}-\bar{y}_{i}=k_{i}\\{\bar{e}_{i(s)}-\bar{e}_{i(r)}\\}+k_{i}O_{p}(n_{L}^{-1}+g^{-1/2}n_{L}^{-1/2}),\\\
&\hat{M}_{i}^{\text{clp}}-\bar{y}_{i}=k_{i}\\{\bar{e}_{i(s)}-\bar{e}_{i(r)}\\}+O_{p}(n_{L}^{-1}+g^{-1/2}n_{L}^{-1/2})\quad\text{and}\\\
&\hat{M}_{i}^{\text{clp}}-\dot{\eta}_{i}=\bar{e}_{i(s)}+O_{p}(n_{L}^{-1}+g^{-1/2}n_{L}^{-1/2}).\end{split}$

###### Proof.

From (2.5) and (2.3), we can write

$\hat{M}_{i}^{\text{sam}}-\bar{y}_{i}=k_{i}\left\\{\mathbf{u}_{i}^{T}(\hat{\boldsymbol{\xi}}-\dot{\boldsymbol{\xi}})+\bar{\mathbf{x}}_{i(r)}^{(w)T}(\hat{\boldsymbol{\beta}}_{2}-\dot{\boldsymbol{\beta}}_{2})+\hat{\alpha}_{i}-\alpha_{i}+\bar{e}_{i(r)}\right\\}.$

Using the approximation

$\hat{\alpha}_{i}=\alpha_{i}+\bar{e}_{i(s)}-\mathbf{u}_{i}^{T}(\hat{\boldsymbol{\xi}}-\dot{\boldsymbol{\xi}})+O_{p}(n_{L}^{-1}+g^{-1/2}n_{L}^{-1/2})$
(3.1)

obtained by Lyu & Welsh (2021) in the proof of their Theorem 2 (see also their
Supplementary Material page 14), and the result from Theorem 1 that
$\hat{\boldsymbol{\beta}}_{2}-\dot{\boldsymbol{\beta}}_{2}=O_{p}(n^{-1/2})=O_{p}(g^{-1/2}n_{L}^{-1/2})$,
as $g,n_{L}\to\infty$, we obtain the approximation

$\begin{split}\hat{M}_{i}^{\text{sam}}-\bar{y}_{i}&=k_{i}\left\\{\bar{e}_{i(s)}-\bar{e}_{i(r)}+\bar{\mathbf{x}}_{i(r)}^{(w)T}(\hat{\boldsymbol{\beta}}_{2}-\dot{\boldsymbol{\beta}}_{2})+O_{p}(n_{L}^{-1}+g^{-1/2}n_{L}^{-1/2})\right\\}\\\
&=k_{i}\\{\bar{e}_{i(s)}-\bar{e}_{i(r)}\\}+k_{i}O_{p}(n_{L}^{-1}+g^{-1/2}n_{L}^{-1/2}).\end{split}$

Similarly, from (2.6) and (2.3), we have

$\hat{M}_{i}^{\text{clp}}-\bar{y}_{i}=\mathbf{u}_{i}^{T}(\hat{\boldsymbol{\xi}}-\dot{\boldsymbol{\xi}})+\bar{\mathbf{x}}^{(w)T}_{i}(\hat{\boldsymbol{\beta}}_{2}-\dot{\boldsymbol{\beta}}_{2})+(\hat{\alpha}_{i}-\alpha_{i})-\bar{e}_{i}.$

Substituting the approximation (3.1) for $\hat{\alpha}_{i}$, we then obtain

$\begin{split}\hat{M}_{i}^{\text{clp}}-\bar{y}_{i}&=\bar{e}_{i(s)}-\bar{e}_{i}+\bar{\mathbf{x}}_{i(r)}^{(w)T}(\hat{\boldsymbol{\beta}}_{2}-\dot{\boldsymbol{\beta}}_{2})+O_{p}(n_{L}^{-1}+g^{-1/2}n_{L}^{-1/2})\\\
&=\bar{e}_{i(s)}-\bar{e}_{i}+O_{p}(n_{L}^{-1}+g^{-1/2}n_{L}^{-1/2})\\\
&=k_{i}\\{\bar{e}_{i(s)}-\bar{e}_{i(r)}\\}+O_{p}(n_{L}^{-1}+g^{-1/2}n_{L}^{-1/2}),\end{split}$
(3.2)

because
$k_{i}\bar{e}_{i(s)}-k_{i}\bar{e}_{i(r)}=k_{i}\bar{e}_{i(s)}-\\{\bar{e}_{i}-(1-k_{i})\bar{e}_{i(s)}\\}=\bar{e}_{i(s)}-\bar{e}_{i}$.

The final result follows from the fact that
$\bar{y}_{i}=\dot{\eta}_{i}+\bar{e}_{i}$, so we have

$\begin{split}&\hat{M}_{i}^{\text{clp}}-\dot{\eta}_{i}=\hat{M}_{i}^{\text{clp}}-(\bar{y}_{i}-\bar{e}_{i})=\bar{e}_{i(s)}+O_{p}(n_{L}^{-1}+g^{-1/2}n_{L}^{-1/2}),\end{split}$

using (3.2). ∎

We can also consider using the estimator of the small area mean to estimate
$\dot{\eta}_{i}$. We obtain the same leading term as for the estimator of the
conditional linear predictor because
$k_{i}\\{\bar{e}_{i(s)}-\bar{e}_{i(r)}\\}+\bar{e}_{i}=\bar{e}_{i(s)}-(n_{i}/N_{i})\bar{e}_{i(s)}-k_{i}\bar{e}_{i(r)}+\bar{e}_{i}=\bar{e}_{i(s)}$.
However, it would be unusual in practice to use the estimator of the small
area mean for this purpose, so we only state the formal result for the
estimator of the conditional linear predictor. The estimators of the small
area mean and the conditional linear predictor have the same leading terms (so
are asymptotically equivalent to first order. However the remainders for the
two estimators are different, showing that there can be higher order
differences between them. In particular, If we sample the entire $i$th cluster
(so $n_{i}=N_{i}$ and $k_{i}=0$), we have
$\hat{M}_{i}^{\text{sam}}-\bar{y}_{i}=0$ but
$\hat{M}_{i}^{\text{clp}}-\bar{y}_{i}=O_{p}(n_{L}^{-1}+g^{-1/2}n_{L}^{-1/2})$,
so $\hat{M}_{i}^{\text{clp}}-\bar{y}_{i}$ is only asymptotically zero.

The asymptotic distribution of the estimators is given by the following
theorem.

######  Theorem 3.

Suppose Condition A holds, $n_{i}^{1/2}/n_{L}\rightarrow 0$ and
$n_{i}/N_{i}\to f_{i}$ for $0\leq f_{i}<1$. Then as $g,n_{L}\to\infty$, we
have

$\displaystyle
n_{i}^{1/2}(\hat{M}_{i}^{\text{sam}}-\bar{y}_{i})\xrightarrow{D}N(0,(1-f_{i}){\dot{\sigma}_{e}}^{2}),\quad
n_{i}^{1/2}(\hat{M}_{i}^{\text{clp}}-\bar{y}_{i})\xrightarrow{D}N(0,(1-f_{i}){\dot{\sigma}_{e}}^{2}),$
$\displaystyle\mbox{ and }\qquad
n_{i}^{1/2}(\hat{M}_{i}^{\text{clp}}-\dot{\eta}_{i})\xrightarrow{D}N(0,{\dot{\sigma}_{e}}^{2}).$

###### Proof.

Since the $e_{ij}$ are independent, $\bar{e}_{i(s)}$ and $\bar{e}_{i(r)}$ are
independent, and the central theorem ensures that
$n_{i}^{1/2}\bar{e}_{i(s)}\xrightarrow{D}N(0,{\dot{\sigma}_{e}}^{2})$ and
$(N_{i}-n_{i})^{1/2}\bar{e}_{i(r)}\xrightarrow{D}N(0,{\dot{\sigma}_{e}}^{2})$.
Provided $f_{i}<1$, it follows that

$\displaystyle
n_{i}^{1/2}\bar{e}_{i(r)}=k_{i}^{-1/2}(1-k_{i})^{1/2}(N_{i}-n_{i})^{1/2}\bar{e}_{i(r)}\xrightarrow{D}N(0,(1-f_{i})^{-1}f_{i}{\dot{\sigma}_{e}}^{2}).$

Consequently,

$\displaystyle
k_{i}n_{i}^{1/2}(\bar{e}_{i(s)}-\bar{e}_{i(r)})\xrightarrow{D}N(0,(1-f_{i})^{2}\\{(1-f_{i})^{-1}f_{i}+1\\}{\dot{\sigma}_{e}}^{2})=N(0,(1-f_{i}){\dot{\sigma}_{e}}^{2}),$

and the theorem follows from Theorem 2. ∎

Therem 3 allows us to reach several interesting conclusions.

  * i)

The asymptotic distribution of both estimators is the distribution of the
target characteristic of interest $F_{\bar{y}_{i}}$ or $F_{\dot{\eta}_{i}}$;
Therem 3 suggests that a better approximation is the distribution of
$k_{i}(\bar{e}_{i(s)}-\bar{e}_{i(r)})$ plus $\bar{y}_{i}$ when the target is
$\bar{y}_{i}$, or the distribution of $\bar{e}_{i(s)}$ plus $\dot{\eta}_{i}$
when the target is $\dot{\eta}_{i}$. The first pair of random variables are
uncorrelated while the second are independent (so the distribution is the
convolution of the $N(0,n_{i}^{-1}\dot{\sigma}_{e}^{2})$ distribution with
$F_{\dot{\eta}_{i}}$). These distributions are not the same in general, but
when the random effects and error in the model are all normally distributed,
these approximations are the same, being
$N(\mathbf{u}_{i}^{T}\dot{\boldsymbol{\xi}}+\bar{\mathbf{x}}_{i}^{(w)T}\dot{\boldsymbol{\beta}}_{2},\,\dot{\sigma}_{\alpha}^{2}+N_{i}^{-1}{\dot{\sigma}_{e}}^{2}+(1-f_{i})n_{i}^{-1}{\dot{\sigma}_{e}}^{2})=N(\mathbf{u}_{i}^{T}\dot{\boldsymbol{\xi}}+\bar{\mathbf{x}}_{i}^{(w)T}\dot{\boldsymbol{\beta}}_{2},\,\dot{\sigma}_{\alpha}^{2}+n_{i}^{-1}{\dot{\sigma}_{e}}^{2})$
and
$N(\mathbf{u}_{i}^{T}\dot{\boldsymbol{\xi}}+\bar{\mathbf{x}}_{i}^{(w)T}\dot{\boldsymbol{\beta}}_{2},\,\dot{\sigma}_{\alpha}^{2}+n_{i}^{-1}{\dot{\sigma}_{e}}^{2})$,
respectively.

  * ii)

The asymptotic mean squared error of $\hat{M}_{i}^{\text{clp}}$ for estimating
$\dot{\eta}_{i}$ is $n_{i}^{-1}\dot{\sigma}_{e}^{2}$ which is greater than or
equal to $n_{i}^{-1}(1-f_{i})\dot{\sigma}_{e}^{2}$, the asymptotic mean
squared error of $\hat{M}_{i}^{\text{clp}}$ (or $\hat{M}_{i}^{\text{sam}}$)
for estimating $\bar{y}_{i}$. Thus, using the asymptotic mean squared error of
$\hat{M}_{i}^{\text{clp}}$ for estimating $\dot{\eta}_{i}$ when we are
estimating $\bar{y}_{i}$ is conservative. This gives an intuitive explanation
for why we might expect the Prasad-Rao estimator of the mean squared error for
estimating $\dot{\eta}_{i}$ (Prasad & Rao, 1990, see the Appendix) to tend to
be greater than our mean squared error estimator.

  * iii)

We can also use Theorem 1 to describe the rate at which we can estimate the
asymptotic mean squared errors. For example, we have
$\widehat{\operatorname{MSE}}_{\text{LW},i}-\operatorname{MSE}_{\text{LW},i}=n_{i}^{-1}k_{i}({{\hat{\sigma}_{e}}}^{2}-{\dot{\sigma}_{e}}^{2}),$
so it follows from Theorem 1 that
$k_{i}^{-1/2}n_{i}^{1/2}n^{1/2}\left(\widehat{\operatorname{MSE}}_{\text{LW},i}-\operatorname{MSE}_{\text{LW},i}\right)=n^{1/2}({{\hat{\sigma}_{e}}}^{2}-{\dot{\sigma}_{e}}^{2})\xrightarrow{D}N(0,\operatorname{E}e_{11}^{4}-{\dot{\sigma}_{e}}^{4}).$

  * iv)

Theorem 3 establishes that the prediction intervals (2.9) and (2.10) have the
correct asymptotic level.

  * v)

The result for estimating $\dot{\eta}_{i}$ also holds when $f_{i}=1$. In this
case, for estimating $\bar{y}_{i}$, provided $N_{i}-n_{i}\rightarrow\infty$,
we have
$(N_{i}-n_{i})^{1/2}(\hat{M}_{i}^{\text{sam}}-\bar{y}_{i})\xrightarrow{P}0$
and
$(N_{i}-n_{i})^{1/2}(\hat{M}_{i}^{\text{clp}}-\bar{y}_{i})\xrightarrow{P}0.$

##  4 Simulation study

In this section we present results from a model-based simulation study to
compare the performance of the prediction intervals (2.9) and (2.10) with the
Prasad-Rao interval (see Supplementary) based on the mixed model estimator of
the small area mean (Sam), the mixed model estimator of the conditional linear
predictor (Clp) and the proposed (LW) and Prasad-Rao (PR) estimators of their
(prediction) root mean squared errors (RMSEs) discussed in the preceding
sections.

We generated population data with $g\in\\{15,30,50\\}$ small areas and $N_{i}$
units in each small area by making area 1 the smallest area ($N_{1}=N_{L}=40$)
and then setting the remaining $N_{i}$ equal to the integer parts of $g-1$
independent uniform $[40,400]$ random variables. The $N_{i}$ were generated
once for each simulation setting so that each setting involved populations of
fixed $N_{i}$ and hence fixed size $N=\sum_{i=1}^{g}N_{i}$. For each setting,
we generated $N$ population values of an auxiliary variable $x_{ij}$ with a
cluster structure by setting $x_{ij}=3+2u_{i}+4v_{ij}$, where $u_{i}$ and
$v_{ij}$ are independent standard normal random variables. We centered the
$x_{ij}$ about their small area means $\bar{x}_{i}$ to obtain the within small
area variable $x_{ij}-\bar{x}_{i}$ and also included $\bar{x}_{i}$ as a
between small area variable. The $N$ population values for $y$ were generated
from the model

$y_{ij}=\beta_{0}+\beta_{1}\bar{x}_{i}+\beta_{2}(x_{ij}-\bar{x}_{i})+\alpha_{i}+e_{ij},\quad
j=1,\ldots,N_{i},\,i=1,\ldots,g,$ (4.1)

where $\\{\alpha_{i}\\}$ were generated independently from $F_{\alpha}$ with
$\operatorname{E}(\alpha_{i})=0$ and
$\operatorname{Var}(\alpha_{i})=\sigma_{\alpha}^{2}$, and independently,
$\\{e_{ij}\\}$ were generated independently from $F_{e}$ with
$\operatorname{E}(e_{ij})=0$ and $\operatorname{Var}(e_{ij})=\sigma_{e}^{2}$.
We set the true regression parameters $\dot{\boldsymbol{\beta}}=[5,7,3]^{T}$,
the variance components ${\dot{\sigma}_{\alpha}}^{2}\in\\{4,64\\}$ and
${\dot{\sigma}_{e}}^{2}\in\\{25,100\\}$, and the distributions
$F_{\alpha}=N(0,\dot{\sigma}_{\alpha}^{2})$ or
$F_{\alpha}=0.3N(0.5,1)+0.7N\big{(}\mu,\\{\dot{\sigma}_{\alpha}^{2}-0.375-0.7\mu^{2}\\}/0.7\big{)}$,
and $F_{e}=N(0,\dot{\sigma}_{e}^{2})$ or
$F_{e}=0.3N(0.5,1)+0.7N(\dot{\mu},\\{\dot{\sigma}_{e}^{2}-0.375-0.7\dot{\mu}^{2}\big{\\}}/0.7)$
with $\dot{\mu}=-0.3\times 0.5/0.7$. The $3$ values of $g$ and $2$ for each of
${\dot{\sigma}_{\alpha}}^{2}$, ${\dot{\sigma}_{e}}^{2}$, $F_{\alpha}$ and
$F_{e}$ produced $48$ different simulation settings.

For each of the $48$ simulation settings, we generated $1000$ populations and
then selected one sample via simple random sampling without replacement from
each population. We set the sample size in each small area to $n_{i}=25$ if
$N_{i}<50$, $n_{i}=\lfloor 0.5*N_{i}\rfloor$ if $50\leq N_{i}\leq 100$, and
$n_{i}=\lfloor 0.25*N_{i}\rfloor$ if $N_{i}>100$, where $\lfloor\quad\rfloor$
is the integer part function, and selected the units in each area
independently by simple random sampling without replacement. For each sample,
we fitted the model (4.1) using REML in lmer and computed the mixed model
estimators of the small area means (Sam) and conditional linear predictors
(Clp) and the root mean squared error estimates (LW and PR). We also computed
the $95\%$ prediction intervals based on these estimates described in (2.9)
and (2.10). We compared the performance of the intervals by examining the
empirical coverage (Cvge) and the relative expected length (RLen) of the
prediction intervals.

For each simulation setting, we report the size $N_{i}$ and the sample size
$n_{i}$ for each small area. For the three intervals, for every small area we
computed the empirical coverage probabilities (Cvge) and the average relative
lengths (Rlen) defined as

$\mbox{Rlen}=|\mbox{Alen}-\mbox{RMSE}_{S}|/\mbox{RMSE}_{S},$

where Alen is the average of the estimates of the root mean squared errors of
the estimators and RMSES is the square root of the average of the squared
differences between the estimator and the true value over the $1000$
simulations. (We suppress the subscript $i$ for simplicity.) The measure Alen
is proportional to the average length of the intervals and RLen measures how
close the lengths of the intervals are to what they should be. Note that ALen
is the same for Sam-LW and Clp-LW because they have the same root mean squared
error estimators, but RLen is different for these two intervals because the
different point estimators have different RMSES. On the other hand, Rlen is
proportional to Alen for Clp-LW and Clp-PR because these both have the same
RMSES. This means that we can evaluate the difference between LW and PR by
comparing Clp-LW and Clp-PR; we can also gain insight into the difference
between using Sam and Clp by comparing Sam-LW with Clp-LW.

The full set of results is given in the Supplementary Material; we present and
discuss illustrative cases below. Table 1 shows the empirical coverage and the
relative length of the three intervals for the setting with variances
${\dot{\sigma}_{\alpha}}^{2}=4$ and ${\dot{\sigma}_{e}}^{2}=100$ (weak within
cluster correlation $0.04$) when $\alpha_{i}$ and $e_{ij}$ have normal
distributions; Table 2 shows the results for the setting with variances
${\dot{\sigma}_{\alpha}}^{2}=64$ and ${\dot{\sigma}_{e}}^{2}=100$ (moderate
within cluster correlation $0.39$) when $\alpha_{i}$ has a mixture
distribution and $e_{ij}$ has a normal distribution; results for the remaining
48 simulation settings are similar. The areas are presented and labeled in
order of increasing size. Simulation standard errors for the coverage
probabilities can be obtained as
$\\{\text{Cvge}(1-\text{Cvge})/1000\\}^{1/2}$; they are approximately 0.008 or
smaller.

Table 1: Simulated coverage and length of prediction intervals when $\alpha_{i}$ and $e_{ij}$ have normal distributions with variances ${\dot{\sigma}_{\alpha}}^{2}=4$ and ${\dot{\sigma}_{e}}^{2}=100$, respectively. |  | Sam-LW |  | Clp-LW |  | Clp-PR
---|---|---|---|---|---|---
Area | N_i | n_i |  | Cvge | Rlen |  | Cvge | Rlen |  | Cvge | Rlen
1 | 40 | 20 |  | 0.965 | 0.10 |  | 0.925 | 0.09 |  | 0.941 | 0.18
2 | 51 | 26 |  | 0.965 | 0.08 |  | 0.966 | 0.05 |  | 0.976 | 0.30
3 | 82 | 41 |  | 0.962 | 0.08 |  | 0.934 | 0.05 |  | 0.987 | 0.46
4 | 86 | 43 |  | 0.980 | 0.19 |  | 0.975 | 0.13 |  | 0.976 | 0.33
5 | 100 | 25 |  | 0.971 | 0.17 |  | 0.965 | 0.10 |  | 0.973 | 0.33
6 | 110 | 28 |  | 0.976 | 0.18 |  | 0.967 | 0.13 |  | 0.961 | 0.27
7 | 113 | 28 |  | 0.972 | 0.17 |  | 0.97 | 0.13 |  | 0.978 | 0.30
8 | 120 | 30 |  | 0.967 | 0.15 |  | 0.966 | 0.11 |  | 0.967 | 0.29
9 | 122 | 30 |  | 0.971 | 0.09 |  | 0.933 | 0.06 |  | 0.969 | 0.32
10 | 135 | 34 |  | 0.962 | 0.07 |  | 0.931 | 0.07 |  | 0.993 | 0.45
11 | 141 | 35 |  | 0.973 | 0.15 |  | 0.966 | 0.11 |  | 0.975 | 0.36
12 | 147 | 37 |  | 0.963 | 0.10 |  | 0.962 | 0.06 |  | 0.988 | 0.35
13 | 150 | 38 |  | 0.967 | 0.11 |  | 0.962 | 0.06 |  | 0.968 | 0.24
14 | 152 | 38 |  | 0.976 | 0.11 |  | 0.968 | 0.06 |  | 0.979 | 0.29
15 | 175 | 44 |  | 0.967 | 0.13 |  | 0.96 | 0.10 |  | 0.975 | 0.36
Table 2: Simulated coverage and length of prediction intervals when $\alpha_{i}$ has a mixture distribution and $e_{ij}$ has a normal distribution with variances ${\dot{\sigma}_{\alpha}}^{2}=64$ and ${\dot{\sigma}_{e}}^{2}=100$, respectively. Method |  | Sam-LW |  | Clp-LW |  | Clp-PR
---|---|---|---|---|---|---
Area | N_i | n_i |  | Cvge | Rlen |  | Cvge | Rlen |  | Cvge | Rlen
1 | 40 | 20 |  | 0.957 | 0.05 |  | 0.934 | 0.03 |  | 0.985 | 0.50
2 | 45 | 20 |  | 0.950 | 0.03 |  | 0.937 | 0.03 |  | 0.987 | 0.42
3 | 58 | 29 |  | 0.957 | 0.01 |  | 0.945 | 0.06 |  | 0.985 | 0.47
4 | 65 | 32 |  | 0.965 | 0.02 |  | 0.945 | 0.03 |  | 0.993 | 0.51
5 | 81 | 40 |  | 0.950 | 0.01 |  | 0.936 | 0.06 |  | 0.989 | 0.46
6 | 85 | 42 |  | 0.957 | 0.01 |  | 0.947 | 0.02 |  | 0.995 | 0.52
7 | 103 | 26 |  | 0.960 | 0.05 |  | 0.949 | 0.03 |  | 0.984 | 0.31
8 | 109 | 27 |  | 0.965 | 0.06 |  | 0.960 | 0.04 |  | 0.986 | 0.31
9 | 115 | 29 |  | 0.957 | 0.06 |  | 0.952 | 0.04 |  | 0.991 | 0.33
10 | 150 | 38 |  | 0.963 | 0.08 |  | 0.959 | 0.07 |  | 0.984 | 0.36
11 | 151 | 38 |  | 0.946 | 0.01 |  | 0.945 | 0.01 |  | 0.977 | 0.28
12 | 162 | 40 |  | 0.952 | 0.06 |  | 0.951 | 0.05 |  | 0.986 | 0.33
13 | 163 | 41 |  | 0.946 | 0.02 |  | 0.947 | 0.01 |  | 0.988 | 0.28
14 | 180 | 45 |  | 0.947 | 0.00 |  | 0.947 | 0.01 |  | 0.985 | 0.25
15 | 193 | 48 |  | 0.961 | 0.06 |  | 0.958 | 0.06 |  | 0.989 | 0.34

The simulation results show that our asymptotic results based on both $g$ and
$n_{L}$ going to infinity provide useful approximations that work well even
when $g=15$ and $n_{L}=20$. The empirical coverages of all three intervals are
close to the nominal level and tend to the nominal level as $g$ and $n_{L}$
increase, confirming our large $g$ and $n_{L}$ asymptotic results. Comparing
the empirical coverage and relative length for Clp-LW with Clp-PR, we see that
the Prasad & Rao (1990) mean squared error estimator PR is typically larger
than our proposed estimator LW, so Clp-PR is more conservative (with wider
intervals) than Clp-LW. The differences in performance between Sam-LW and Clp-
LW are smaller with both performing well.

##  5 Consumer expenditure on fresh milk products

We obtained data from the Dairy Survey component of the 2002 Consumer
Expenditure Survey conducted by the U.S. Bureau of the Census for the U.S.
Bureau of Labor Statistics; the data are available from
https://www.bls.gov/cex/pumd_data.htm. We treated the consumer expenditure on
fresh milk products (MILKPROD) as the survey variable of interest and
considered the problem of estimating the average consumer expenditure on fresh
milk products in different states (small areas). We used the total expenditure
on food (FOODTOT), the number of persons under age 18 in the family (PERSLT18)
and the total family income before taxes in the last 12 months (FINCBEFX) as
the auxiliary variables. This data set is similar to that used in Arora &
Lahiri (1997); they used the earlier 1989 survey, focussed on the expenditure
on fresh whole milk, and potentially used different auxiliary variables. If we
knew the means of the auxiliary variables for each state, we could use our
methods to estimate the average expenditure on fresh milk products in 2002 in
each state. As we do not have this information, we instead treated the data
set as a pseudo-population and sampled from it. We repeated this sampling
$1000$ times, effectively implementing a design-based simulation from our
fixed population to evaluate the design-based properties of our proposed
model-based methods.

In creating the population, we discarded 6 states with fewer than 10
observations, leaving us with $N=4022$ observations from $g=34$ states with
between $N_{L}=36$ and $N_{U}=397$ observations from each state. We centered
the auxiliary variables about their cluster means (adding $cent$ to their
variable name) and then included the cluster means (adding $avg$ to the
variable name) as between state variables so that we have $p_{b}=3$ plus
$p_{w}=3$ auxiliary variables. The cluster means give the average per family
of each variable for each state. We selected the $1000$ samples independently
by simple random sampling without replacement from each state with $n_{L}=20$
by setting $n_{i}=20$ if $N_{i}<50$, $n_{i}=\lfloor 0.5*N_{i}\rfloor$ if
$50\leq N_{i}\leq 100$, and $n_{i}=\lfloor 0.25*N_{i}\rfloor$ if $N_{i}>100$,
where $\lfloor\quad\rfloor$ is the integer part function. In each sample, we
fitted the model (5.1)

$\displaystyle MILKPROD_{ij}=$
$\displaystyle\beta_{0}+\beta_{1}FOODTOTavg_{i}+\beta_{2}PERSLT18avg_{i}+\beta_{3}FINCBEFXavg_{i}$
(5.1)
$\displaystyle+\beta_{4}FOODTOTcent_{ij}+\beta_{5}PERSLT18cent_{ij}+\beta_{6}FINCBEFXcent_{ij}+\alpha_{i}+e_{ij},$

using lmer and then used the parameter estimates to compute the mixed model
estimators of the small area means (Sam) and the conditional linear predictors
(Clp) and their prediction mean squared errors. We then computed the $95\%$
model-based prediction intervals (2.9) and (2.10). Similar to the model-based
simulation, we examined the empirical design-coverage (Cvge) and the relative
design-expected length (Rlen) of the prediction intervals; the results over
$1000$ samples together with the standardised population EBLUPs
$\hat{\alpha}_{i}/\hat{\sigma}_{\alpha}$ for each state are shown in Table 3.
The relative design-bias and the design RMSEs of Sam and Clp, together with
the design-averages of the LW and PR estimators of the RMSEs are available in
the Supplementary Information.

Table 3: Simulated design-coverage and design-expected length of nominal $95\%$ confidence intervals for the average consumer expenditure on fresh milk products in each state in 2002. ∗ identifies states in Group 3 and † identifies states in Group 2.  $\operatorname{RMSE}$ |  | Sam-LW |  | Clp-LW |  | Clp-PR
---|---|---|---|---|---|---
STATE | N_i | n_i | $\hat{\alpha}_{i}/\hat{\sigma}_{\alpha}$ |  | Cvge | Rlen |  | Cvge | Rlen |  | Cvge | Rlen
16 | 36 | 20 | 0.14 |  | 0.997 | 0.57 |  | 1.000 | 0.94 |  | 0.996 | 1.13
50† | 37 | 20 | 0.52 |  | 0.912 | 0.12 |  | 0.926 | 0.25 |  | 0.855 | 0.19
31 | 42 | 20 | -0.33 |  | 1.000 | 0.58 |  | 1.000 | 0.28 |  | 0.981 | 0.29
22† | 43 | 20 | -0.48 |  | 1.000 | 0.56 |  | 1.000 | 0.08 |  | 0.917 | 0.08
21 | 44 | 20 | -0.17 |  | 1.000 | 0.92 |  | 1.000 | 0.85 |  | 0.992 | 0.84
15 | 45 | 20 | 0.10 |  | 1.000 | 0.85 |  | 1.000 | 1.22 |  | 0.995 | 1.18
32† | 46 | 20 | 0.71 |  | 0.921 | 0.14 |  | 0.832 | 0.35 |  | 0.698 | 0.37
37† | 47 | 20 | -0.48 |  | 1.000 | 0.36 |  | 1.000 | 0.05 |  | 0.915 | 0.09
1 | 52 | 26 | -0.01 |  | 1.000 | 1.54 |  | 1.000 | 2.06 |  | 1.000 | 2.53
45 | 53 | 26 | 0.04 |  | 0.973 | 0.14 |  | 0.997 | 0.43 |  | 0.991 | 0.63
2∗ | 58 | 29 | 1.22 |  | 0.795 | 0.32 |  | 0.414 | 0.56 |  | 0.571 | 0.47
9∗ | 59 | 30 | -1.08 |  | 0.745 | 0.37 |  | 0.360 | 0.59 |  | 0.554 | 0.49
41 | 65 | 32 | 0.12 |  | 0.984 | 0.25 |  | 0.997 | 0.70 |  | 0.998 | 1.10
49 | 67 | 34 | 0.00 |  | 0.995 | 0.48 |  | 1.000 | 0.75 |  | 0.999 | 1.24
18 | 71 | 36 | -0.27 |  | 0.971 | 0.02 |  | 0.996 | 0.25 |  | 0.996 | 0.63
27∗ | 76 | 38 | 1.31 |  | 0.740 | 0.39 |  | 0.449 | 0.56 |  | 0.658 | 0.43
8† | 82 | 41 | 0.66 |  | 0.867 | 0.23 |  | 0.840 | 0.28 |  | 0.931 | 0.03
13 | 93 | 46 | -0.09 |  | 0.991 | 0.30 |  | 0.998 | 0.76 |  | 1.000 | 1.43
24† | 94 | 47 | -0.91 |  | 0.933 | 0.06 |  | 0.790 | 0.35 |  | 0.908 | 0.09
29 | 98 | 49 | -0.26 |  | 1.000 | 0.78 |  | 1.000 | 0.75 |  | 1.000 | 1.47
53 | 99 | 50 | 0.24 |  | 0.998 | 0.63 |  | 0.999 | 0.79 |  | 1.000 | 1.56
55 | 119 | 30 | -0.49 |  | 1.000 | 0.47 |  | 1.000 | 0.31 |  | 0.976 | 0.30
51 | 122 | 30 | 0.01 |  | 0.998 | 0.80 |  | 0.998 | 1.21 |  | 1.000 | 1.19
25 | 126 | 32 | 0.06 |  | 1.000 | 1.06 |  | 1.000 | 1.55 |  | 1.000 | 1.59
4 | 133 | 33 | 0.13 |  | 0.999 | 0.55 |  | 1.000 | 0.79 |  | 0.997 | 0.83
26 | 139 | 35 | -0.29 |  | 0.997 | 0.51 |  | 1.000 | 0.60 |  | 0.996 | 0.68
34 | 160 | 40 | -0.16 |  | 0.998 | 0.61 |  | 1.000 | 0.84 |  | 1.000 | 1.00
17 | 161 | 40 | -0.28 |  | 1.000 | 0.81 |  | 1.000 | 0.86 |  | 1.000 | 1.02
39 | 229 | 57 | -0.16 |  | 0.999 | 0.64 |  | 0.999 | 0.87 |  | 1.000 | 1.23
42 | 261 | 65 | -0.57 |  | 0.981 | 0.15 |  | 0.988 | 0.11 |  | 0.993 | 0.36
36 | 280 | 70 | 0.02 |  | 0.975 | 0.16 |  | 0.984 | 0.30 |  | 0.997 | 0.63
12 | 283 | 71 | 0.58 |  | 0.982 | 0.12 |  | 0.979 | 0.09 |  | 0.991 | 0.36
48 | 305 | 76 | -0.44 |  | 0.986 | 0.27 |  | 0.991 | 0.30 |  | 1.000 | 0.65
6 | 397 | 99 | 0.59 |  | 0.964 | 0.05 |  | 0.968 | 0.04 |  | 0.998 | 0.37

To interpret the results, we partition the states into three groups: Group 3
with three states $\\{2,9,27\\}$ for which the design-coverage of all three
intervals is well below the nominal level; Group 2 with six states
$\\{8,22,24,32,37,50\\}$ for which at least one interval has design-coverage
below the nominal level, but not all intervals perform poorly; and Group 3
with the remaining twenty-five states for which the design-coverage of all
three intervals is above the nominal level.

For the states in Group 1, the design-coverages tend to be conservative and
similar across all three intervals. The design-biases for Sam are smaller than
for Clp except in state 49 where the design-biases are very similar. The
design-average LW and PR estimators of the RMSE are greater than the
simulation RMSEs of the Sam and Clp estimators, respectively. The design-
average LW estimator is smaller than the design-average PR estimator in all
except 4 states (15, 21, 51 and 55) where PR has a slightly smaller design-
average than LW.

Group 3 comprises the three states for which all three intervals have design-
coverage well below the nominal level. Sam-LW has the best design-coverage,
followed by Clp-PR and then Clp-LW. Both point estimates Sam and Clp are
design-biased with Clp having greater design-bias than Sam. The design-average
RMSE estimator PR is on design-average greater than LW, explaining why Clp-PR
has better design-coverage than Clp-LW. However, Sam-LW has better design-
coverage than Clp-PR because PR is not large enough to overcome the greater
design-bias in Clp than in Sam.

Group 2 comprises the states for which at least one but not all three
intervals have design-coverage below the nominal level. Usually Clp-LW but
sometimes Clp-PR performs most poorly; Sam-LW usually performs better, but in
state 8 Sam-LW also has low design-coverage. Clp has larger design-bias than
Sam for Group 2. The design-average PR estimator of the RMSE is smaller than
the simulation RMSE of the Clp estimator, while the design-average LW
estimators of the RMSE is close to the simulation RMSEs of the Sam estimators.
The RMSE estimator LW is on design-average smaller than PR; the larger PR
reduces the effect of the design-bias of Clp on the design-coverage in state
8.

These results show that the Sam-LW intervals are generally preferable to the
Clp intervals for making design-based inference in our consumer expenditure on
fresh milk products population. Nonetheless, we did not expect to see results
like those in Group 3 so we felt that further exploration of why these results
occur would be interesting and useful.

As we have access to the whole population (which is not usual in practice), we
are able to explore the population. We used lmer from the R package lme4 to
fit the nested error regression model (5.1) to the population data. Figure 1
shows a normal QQ-plot of the EBLUPs (with approximate $95\%$ prediction
intervals computed as in Lyu & Welsh (2021)) and a normal QQ-plot of the
errors from the fitted model. Figure 1 shows that it is plausible to treat the
random effects as approximately normally distributed; the errors have an
asymmetric long-tailed distribution and it is not plausible to treat the
errors in the model as normally distributed. In a model-based analysis, rather
than simply relying on the asymptotic theory, we could consider using
transformations to improve the fit, make predictions on the transformed scale
and then invert the prediction intervals (possibly adjusting for back-
transformation bias). However, for our design-based analysis, as is arguably
usual in practice, we keep the data on the raw scale. Table 4 shows the
parameter estimates and standard errors for the fitted model. Although some
coefficients are not significant, we retain them in the model. We see that
$\hat{\sigma}_{e}^{2}/\hat{\sigma}_{a}^{2}=48.72$ is large so that the within
cluster correlation is very small (approximately $0.02$).

Table 4: Parameter estimates (REML) for modelling the consumer expenditure on fresh milk products population. Effect | Estimate | Std Error
---|---|---
(Intercept) | 2.6219 | 0.7622
FOODTOTavg | 1.1617 | 5.0871
PERSLT18avg | 0.9235 | 0.5526
FINCBEFXavg | 0.0155 | 0.0083
FOODTOTcent | 4.7930 | 0.3656
PERSLT18cent | 0.7098 | 0.0399
FINCBEFXcent | 0.0015 | 0.0010
$\sigma_{\alpha}^{2}$ | 0.1757 |
$\sigma_{e}^{2}$ | 8.5600 |
NOTE: lmer does not compute standard
errors for the variance components.

Figure 1: Normal QQ-plots for the EBLUPs $\hat{\alpha}_{i}$ and estimated
errors $\hat{e}_{ij}$ from the model (5.1) fitted to the consumer expenditure
on fresh milk products in 2002 data.

The normal QQ-plot of the EBLUPs in Figure 1 shows that the Group 3 states
have extreme EBLUPs and the Group 2 states have EBLUPs in the tails of the
distribution but these are mixed in with some EBLUPs from Group 1 states. This
mixing suggests that the EBLUPS alone do not identify the group to which a
state belongs. Another potentially important value is the sample size $n_{i}$
in each state. Figure 2 shows the population standardised EBLUPS
$\hat{\alpha}_{i}/\hat{\sigma}_{\alpha}$ plotted against the sample size
$n_{i}$ for each state; these variables are also included in Table 3. States
plotted in the top and bottom left of the plot (high standardised EBLUPs and
small to moderate sample size) are in Groups 2 and 3; states at the top or
bottom right (small standardised EBLUPs or large standardised EBLUPS with
large sample sizes) are in Group 1 with all the others states. Specifically,
states 6 ($n_{6}=99$), 12 ($n_{12}=71$) and 42 ($n_{42}=65$)) have relatively
extreme EBLUPs but larger sample sizes so these states are in Group 1. This
suggests that both the magnitude of the standardised EBLUPs and the sample
size determine the difficulty of estimating a particular state.

Figure 2: Plot of the population standardised EBLUPS
$\hat{\alpha}_{i}/\hat{\sigma}_{\alpha}$ against sample size in each state for
the consumer expenditure on fresh milk products. Group 1 states are plotted as
circles, Group 2 as triangles and Group 3 as crosses. Selected states are
labelled by state number.

Another possibility we considered is that the within state variances may not
be constant, either because of lack of fit in the mean in some states or
because the within state variance simply differs in different states.
Incorrectly assuming constant variance as we have done in model (5.1) would
mean that the LW estimator is too small in states with large within state
variance, so the intervals should perform poorly in these states. We computed
the sample variance of the errors in each state to produce rough estimates of
the within state variances. These variances range from $1.48$ to $20.46$. In
Group 3, states 2 and 27 have large variances, but state 9 does not. In Group
2, states 8, 32 and 50 have larger error variances (group 8 has the larges
variance of all the states), but states 22, 24 and 37 actually have very small
error variances. States 18 and 36 in Group 1 have variances comparable to the
larger variances in Group 2 while the majority of states have small variances.
The patterns are less consistent than those based on the standardised EBLUPs
and sample size so we focussed more on these variables.

We cannot rule out the possibility that the results of the design-based
simulation using the consumer expenditure on fresh milk products data are at
least in part due to failures of the model (5.1). However, we can explore
whether similar results occur when the model is correct by carrying out a set
of additional design-based simulations. We used the same 48 settings as in our
model-based simulation, but we generated a single population for each setting
and then selected $1000$ samples from it by simple random sampling without
replacement. We fitted the model and computed the same $95\%$ prediction
intervals as before and evaluated their design-based properties. The full set
of results are available in the Supplementary material; we present and discuss
illustrative cases below.

Table 5: Simulated design-coverage and design-average length of confidence intervals when $\alpha_{i}$ and $e_{ij}$ have normal distributions with variances ${\dot{\sigma}_{\alpha}}^{2}=4$ and ${\dot{\sigma}_{e}}^{2}=100$, respectively. ∗ identifies states in Group 3 and † identifies states in Group 2. Method |  | Sam-LW |  | Clp-LW |  | Clp-PR
---|---|---|---|---|---|---
Area | N_i | n_i | $\hat{\alpha}_{i}/\hat{\sigma}_{\alpha}$ |  | Cvge | Rlen |  | Cvge | Rlen |  | Cvge | Rlen
1 | 40 | 20 | 0.37 |  | 0.993 | 0.30 |  | 0.999 | 0.27 |  | 1.000 | 0.65
2 | 40 | 20 | 0.11 |  | 0.997 | 0.59 |  | 1.000 | 1.25 |  | 1.000 | 1.93
3 | 43 | 20 | 0.12 |  | 0.971 | 0.15 |  | 1.000 | 0.69 |  | 1.000 | 1.13
4 | 50 | 25 | 0.48 |  | 0.980 | 0.21 |  | 0.991 | 0.10 |  | 0.999 | 0.53
5 | 61 | 30 | 0.24 |  | 0.994 | 0.36 |  | 1.000 | 0.61 |  | 1.000 | 1.32
6 | 72 | 36 | 0.37 |  | 0.981 | 0.25 |  | 0.992 | 0.34 |  | 1.000 | 1.02
7 | 74 | 37 | 0.51 |  | 0.979 | 0.09 |  | 0.987 | 0.09 |  | 0.999 | 0.65
8∗ | 79 | 40 | -1.45 |  | 0.860 | 0.25 |  | 0.554 | 0.50 |  | 0.901 | 0.22
9 | 81 | 40 | -0.60 |  | 0.978 | 0.14 |  | 0.981 | 0.06 |  | 1.000 | 0.62
10† | 88 | 44 | 1.03 |  | 0.925 | 0.07 |  | 0.847 | 0.30 |  | 0.989 | 0.10
11 | 89 | 44 | 0.71 |  | 0.982 | 0.16 |  | 0.978 | 0.03 |  | 0.999 | 0.50
12 | 92 | 46 | -0.43 |  | 0.992 | 0.35 |  | 0.996 | 0.32 |  | 1.000 | 1.08
13 | 96 | 48 | -0.61 |  | 0.962 | 0.04 |  | 0.971 | 0.01 |  | 1.000 | 0.59
14 | 101 | 25 | 0.14 |  | 0.997 | 0.49 |  | 0.999 | 0.79 |  | 1.000 | 1.03
15∗ | 105 | 26 | 1.13 |  | 0.909 | 0.19 |  | 0.831 | 0.32 |  | 0.880 | 0.22
16 | 113 | 28 | 0.26 |  | 0.992 | 0.32 |  | 0.995 | 0.41 |  | 0.999 | 0.65
17 | 127 | 32 | 0.07 |  | 0.995 | 0.44 |  | 0.999 | 0.71 |  | 1.000 | 1.05
18 | 129 | 32 | 0.54 |  | 0.986 | 0.23 |  | 0.989 | 0.18 |  | 0.999 | 0.42
19 | 131 | 33 | -0.31 |  | 0.992 | 0.30 |  | 0.996 | 0.41 |  | 1.000 | 0.71
20† | 138 | 34 | -0.91 |  | 0.934 | 0.08 |  | 0.914 | 0.18 |  | 0.962 | 0.01
21 | 152 | 38 | 0.15 |  | 0.996 | 0.47 |  | 0.998 | 0.66 |  | 1.000 | 1.06
22† | 162 | 40 | -1.25 |  | 0.927 | 0.16 |  | 0.854 | 0.30 |  | 0.950 | 0.12
23 | 169 | 42 | -0.74 |  | 0.969 | 0.07 |  | 0.971 | 0.01 |  | 0.992 | 0.28
24 | 170 | 42 | -0.01 |  | 0.987 | 0.36 |  | 0.996 | 0.55 |  | 1.000 | 0.96
25 | 173 | 43 | -0.11 |  | 0.995 | 0.47 |  | 0.997 | 0.67 |  | 1.000 | 1.12
26 | 177 | 44 | 0.78 |  | 0.965 | 0.05 |  | 0.966 | 0.01 |  | 0.996 | 0.28
27 | 186 | 46 | -0.08 |  | 0.991 | 0.34 |  | 0.997 | 0.51 |  | 1.000 | 0.93
28 | 190 | 48 | 0.09 |  | 0.997 | 0.42 |  | 0.999 | 0.59 |  | 1.000 | 1.05
29 | 197 | 49 | -0.84 |  | 0.968 | 0.04 |  | 0.962 | 0.03 |  | 0.992 | 0.25
30 | 199 | 50 | 0.23 |  | 0.989 | 0.25 |  | 0.991 | 0.34 |  | 1.000 | 0.74
Table 6: Simulated design-coverage and design-average length of confidence intervals when $\alpha_{i}$ has a mixture distribution and $e_{ij}$ has a normal distribution with variances ${\dot{\sigma}_{\alpha}}^{2}=4$ and ${\dot{\sigma}_{e}}^{2}=100$, respectively. ∗ identifies states in Group 3 and † identifies states in Group 2. Method |  | Sam-LW |  | Clp-LW |  | Clp-PR
---|---|---|---|---|---|---
Area | N_i | n_i | $\hat{\alpha}_{i}/\hat{\sigma}_{\alpha}$ |  | Cvge | Rlen |  | Cvge | Rlen |  | Cvge | Rlen
1∗ | 40 | 20 | -1.27 |  | 0.883 | 0.27 |  | 0.205 | 0.59 |  | 0.509 | 0.5
2 | 44 | 20 | 0.26 |  | 0.997 | 0.49 |  | 1.000 | 0.65 |  | 0.999 | 0.91
3 | 54 | 27 | 0.03 |  | 0.996 | 0.50 |  | 1.000 | 1.13 |  | 1.000 | 1.86
4 | 58 | 29 | -0.24 |  | 0.976 | 0.16 |  | 0.998 | 0.49 |  | 1.000 | 1.04
5 | 68 | 34 | 0.35 |  | 0.978 | 0.13 |  | 0.994 | 0.24 |  | 0.999 | 0.78
6 | 80 | 40 | 0.15 |  | 0.987 | 0.22 |  | 0.998 | 0.56 |  | 1.000 | 1.31
7 | 80 | 40 | -0.68 |  | 0.958 | 0.03 |  | 0.948 | 0.10 |  | 0.996 | 0.34
8† | 91 | 46 | -0.99 |  | 0.971 | 0.04 |  | 0.872 | 0.28 |  | 0.987 | 0.10
9† | 95 | 48 | 0.91 |  | 0.963 | 0.01 |  | 0.887 | 0.24 |  | 0.989 | 0.17
10† | 96 | 48 | 1.04 |  | 0.920 | 0.14 |  | 0.794 | 0.34 |  | 0.967 | 0.01
11 | 101 | 25 | 0.34 |  | 0.999 | 0.45 |  | 1.000 | 0.48 |  | 0.999 | 0.59
12 | 104 | 26 | 0.60 |  | 0.997 | 0.23 |  | 0.999 | 0.11 |  | 0.995 | 0.20
13 | 109 | 27 | -0.15 |  | 0.991 | 0.33 |  | 0.997 | 0.56 |  | 0.999 | 0.70
14 | 109 | 27 | -0.08 |  | 0.994 | 0.48 |  | 0.997 | 0.83 |  | 1.000 | 1.00
15 | 111 | 28 | -0.16 |  | 0.995 | 0.38 |  | 1.000 | 0.56 |  | 1.000 | 0.73
16 | 112 | 28 | -0.38 |  | 1.000 | 0.54 |  | 1.000 | 0.55 |  | 1.000 | 0.71
17 | 117 | 29 | 0.30 |  | 0.994 | 0.34 |  | 1.000 | 0.43 |  | 0.999 | 0.60
18 | 119 | 30 | 0.28 |  | 0.995 | 0.42 |  | 0.997 | 0.55 |  | 0.999 | 0.75
19 | 123 | 31 | -0.14 |  | 0.995 | 0.59 |  | 0.999 | 0.89 |  | 1.000 | 1.16
20 | 136 | 34 | 0.51 |  | 0.998 | 0.36 |  | 0.999 | 0.31 |  | 0.999 | 0.52
21 | 137 | 34 | -0.39 |  | 0.994 | 0.46 |  | 0.997 | 0.51 |  | 0.999 | 0.76
22 | 147 | 37 | 0.55 |  | 0.991 | 0.27 |  | 0.996 | 0.24 |  | 0.998 | 0.48
23 | 159 | 40 | 0.19 |  | 0.995 | 0.43 |  | 0.998 | 0.63 |  | 1.000 | 0.98
24∗ | 189 | 47 | -1.15 |  | 0.905 | 0.18 |  | 0.852 | 0.3 |  | 0.938 | 0.12
25 | 190 | 48 | -0.02 |  | 0.991 | 0.39 |  | 0.998 | 0.55 |  | 1.000 | 0.95
26 | 194 | 48 | 0.15 |  | 0.988 | 0.33 |  | 0.996 | 0.51 |  | 1.000 | 0.89
27 | 196 | 49 | -0.98 |  | 0.959 | 0.03 |  | 0.934 | 0.15 |  | 0.980 | 0.08
28 | 199 | 50 | 0.62 |  | 0.986 | 0.19 |  | 0.991 | 0.16 |  | 0.999 | 0.46
29 | 199 | 50 | 0.39 |  | 0.989 | 0.31 |  | 0.993 | 0.36 |  | 0.999 | 0.73
30 | 200 | 50 | -0.08 |  | 0.997 | 0.37 |  | 0.998 | 0.58 |  | 1.000 | 0.99

Tables 5 and 6 show the empirical design-coverage and the design-average
relative length of the intervals for the settings with variances
${\dot{\sigma}_{\alpha}}^{2}=4$ and ${\dot{\sigma}_{e}}^{2}=100$ when $e_{ij}$
has a normal distribution and $\alpha_{i}$ has either a normal distribution or
a mixture distribution; settings with non-normal $e_{ij}$ are included in the
supplementary material. The areas are again presented and labeled in order of
increasing size. The Monte Carlo standard errors for the design-coverage
probabilities are approximately less than 0.01.

In all settings, the model holds and the within cluster variances are
constant. Nonetheless, when $\sigma_{e}^{2}/\sigma_{a}^{2}$ is large (as in
Tables 5 and 6), the kind of results we saw in the previous simulation occur.
In Table 5, there are 2 Group 3 areas, numbers 8
($\hat{\alpha}_{8}/\hat{\sigma}_{\alpha}=-1.45$, $n_{8}=40$) and 15
($\hat{\alpha}_{15}/\hat{\sigma}_{\alpha}=1.13$, $n_{15}=26$), and 3 Group 2
areas, numbers 10 ($\hat{\alpha}_{10}/\hat{\sigma}_{\alpha}=1.03$,
$n_{10}=44$), 20 ($\hat{\alpha}_{20}/\hat{\sigma}_{\alpha}=-0.91$, and
$n_{20}=34$), 22 ($\hat{\alpha}_{22}/\hat{\sigma}_{\alpha}=-1.25$,
$n_{22}=40$). In Table 6, there are 2 Group 3 areas, numbers 1
($\hat{\alpha}_{1}/\hat{\sigma}_{\alpha}=-1.27$, $n_{1}=20$), 24
($\hat{\alpha}_{24}/\hat{\sigma}_{\alpha}=-1.15$, $n_{24}=47$), and 3 Group 2
areas, numbers 8 ($\hat{\alpha}_{8}/\hat{\sigma}_{\alpha}=-0.99$, $n_{8}=46$),
9 ($\hat{\alpha}_{9}/\hat{\sigma}_{\alpha}=0.91$, $n_{9}=48$) and 10
($\hat{\alpha}_{10}/\hat{\sigma}_{\alpha}=1.04$, $n_{10}=48$). On the other
hand, when ${\dot{\sigma}_{e}}^{2}/{\dot{\sigma}_{\alpha}}^{2}$ is not large,
there are no Group 2 or 3 areas. Overall, we see that when
${\dot{\sigma}_{e}}^{2}/{\dot{\sigma}_{\alpha}}^{2}$ is large, areas with
extreme EBLUPs and small to moderate sample sizes are more difficult than
other areas to estimate well inthe design-based framework.

Why does the size of the random effect for an area matter in the design-based
framework but not in the model-based framework? In the model-based framework,
the population and hence the random effect for an area is generated anew for
each replication. This means that in each sample we are estimating a
realisation of an independent and identically distributed random variable so,
over model-based replications, we are estimating the expected value of the
random variable which is zero under the model. In the design-based framework,
the population is fixed and the replication is over independent samples from
this fixed population. Once generated, the random effects are fixed. This
means that in areas with extreme random effects, the EBLUPs are estimating the
expected values of extreme order statistics which are not zero and difficult
to estimate. This then flows through into estimating the area mean of the
survey variable for areas with large random effects. In our design-based
simulation, we used the population EBLUPs to assess the difficulty in
estimation, but in practice, as shown in Lyu & Welsh (2021), we would use the
sample EBLUPs to estimate the population random effects.

The above discussion suggests that when we want to achieve good (model-
assisted) design-based rather than model-based performance, we should treat
the random effects in the model as fixed. To check this intuition (and
indirectly confirm the argument above), we repeated our design-based
simulations treating $\alpha_{i}$ as fixed and examined the empirical design
coverage and the relative design-expected length of the model-based prediction
intervals constructed under the fixed area effects model. Explicitly, the
general form of the approach is to rewrite the model (2.1) and (2.2) with
fixed area effects as a regression model

$y_{ij}=\mathbf{z}_{ij}^{*T}\boldsymbol{\chi}+e_{ij},\qquad\text{for
$j=1,\ldots,N_{i},\,i=1,\ldots,g$,}$ (5.2)

where
$\mathbf{z}_{ij}^{*}=[\mathbf{u}_{i}^{T},\mathbf{x}_{ij}^{(w)T},\mathbf{v}_{i}^{T}]^{T}$
with $\mathbf{v}_{i}$ a $g$-vector of zeros with a one in position $i$, and
$\boldsymbol{\chi}=[\boldsymbol{\xi}^{T},\boldsymbol{\beta}_{2}^{T},\alpha_{1},\ldots,\alpha_{g}]^{T}$
with $\sum_{i=1}^{g}\alpha_{i}=0$. The normal maximum likelihood estimator
$\hat{\boldsymbol{\chi}}$ of $\boldsymbol{\chi}$ is the (constrained) ordinary
least squares estimator which can be computed using lm in R with the sum to
zero constraint on the $\alpha_{i}$. The optimal model-based predictor of
$\bar{y}_{i}$ under (5.2) is the composite estimator

$\bar{y}_{i}^{\text{com-
fixed}}=(1-k_{i})\bar{y}_{i(s)}+k_{i}\bar{\mathbf{z}}_{i(r)}^{*T}\hat{\boldsymbol{\chi}},$

where
$\bar{\mathbf{z}}_{i(r)}^{*}=[\mathbf{u}_{i}^{T},\bar{\mathbf{x}}_{i(r)}^{(w)T},\boldsymbol{\nu}_{i}^{T}]^{T}$
and an approximate $100(1-\epsilon)\%$ prediction interval for $\bar{y}_{i}$
is

$\Big{[}\bar{y}_{i}^{\text{com-
fixed}}-\Phi^{-1}(1-\epsilon/2)k_{i}\left\\{\frac{{{\hat{\sigma}_{e}}}^{2}}{N_{i}-n_{i}}+\bar{\mathbf{z}}_{i(r)}^{*T}\hat{\mathbf{V}}_{\boldsymbol{\chi}}\bar{\mathbf{z}}_{i(r)}^{*}\right\\}^{1/2},\,\bar{y}_{i}^{\text{com-
fixed}}+\Phi^{-1}(1-\epsilon/2)k_{i}\left\\{\frac{{{\hat{\sigma}_{e}}}^{2}}{N_{i}-n_{i}}+\bar{\mathbf{z}}_{i(r)}^{*T}\hat{\mathbf{V}}_{\boldsymbol{\chi}}\bar{\mathbf{z}}_{i(r)}^{*}\right\\}^{1/2}\Big{]},$

where ${{\hat{\sigma}_{e}}}^{2}$ estimates $\sigma_{e}^{2}$, and
$\hat{\mathbf{V}}_{\boldsymbol{\chi}}$ estimates the variance of
$\hat{\boldsymbol{\chi}}$. Similarly, the synthetic estimator is
$\bar{y}^{\text{syn-fixed}}_{i}=\bar{\mathbf{z}}_{i}^{*T}\boldsymbol{\chi}$,
where
$\bar{\mathbf{z}}_{i}^{*}=[\mathbf{u}_{i}^{T},\bar{\mathbf{x}}_{i}^{(w)T},\boldsymbol{\nu}_{i}^{T}]^{T}$,
and an approximate $100(1-\epsilon)\%$ prediction interval for $\bar{y}_{i}$
is

$\Big{[}\bar{y}_{i}^{\text{com-
fixed}}-\Phi^{-1}(1-\epsilon/2)(\bar{\mathbf{z}}_{i}^{*T}\hat{\mathbf{V}}_{\boldsymbol{\chi}}\bar{\mathbf{z}}_{i})^{1/2},\,\bar{y}_{i}^{\text{com-
fixed}}+\Phi^{-1}(1-\epsilon/2)(\bar{\mathbf{z}}_{i}^{*T}\hat{\mathbf{V}}_{\boldsymbol{\chi}}\bar{\mathbf{z}}_{i})^{1/2}\Big{]}.$

The full set of results is included in the Supplementary material. The design-
coverage results for both methods are generally good; the relative design-
expected lengths of the intervals for the composite approach are mostly
smaller than those for the synthetic approach.

##  6 Discussion

In this paper, we considered model-based small area estimation under the
nested error regression model. We discussed two target characteristics of
interest, the small area means and the conditional linear predictors of the
small area means, and the construction of mixed model estimators (EBLUPs) of
these two targets. We established asymptotic linearity results and central
limit theorems for these estimators which allow us to establish asymptotic
equivalences between estimators, to approximate their sampling distributions,
obtain simple expressions for and construct simple estimators of their
asymptotic mean squared errors, and justify asymptotic prediction intervals.
Our new results are established under the asymptotic framework of increasing
numbers of small areas and increasing numbers of units in each area, a
framework that has not previously been applied in small area estimation. We
report model-based simulations that show that these results are applicable in
quite small, finite samples, establishing that they fill important theoretical
gaps and are useful in practice. In particular, our mean squared error
estimator performs as well or better than the widely-used Prasad & Rao (1990)
estimator and is much simpler, so it is easier to interpret and consequently
provides more insight. We also carried out a design-based simulation using
real data on consumer expenditure on fresh milk products to explore the
design-based properties of the mixed model estimators. This simulation
produced some surprising results which we managed to explain and interpret
through analysis of the population and further design-based simulations. The
simulations together highlight under-appreciated differences between the
model- and design-based properties of mixed model estimators in small area
estimation.

Given the extensive literature on small area estimation, it is important to
acknowledge that in this paper we have considered only one of many interesting
and important scenarios. Future work should include applying the asymptotic
approach to area level models, outlier robust estimators and the extensions to
the basic nested error regression model discussed in the Introduction.

## References

  * (1)
  * Arora & Lahiri (1997) Arora, V. & Lahiri, P. (1997), ‘On the superiority of the Bayesian method over the BLUP in small area estimation problems’, Statistica Sinica pp. 1053–1063.
  * 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, 28–36.
  * Chandra & Chambers (2016) Chandra, H. & Chambers, R. (2016), ‘Small area estimation for semicontinuous data’, Biometrical Journal 58(2), 303–319.
  * Chatterjee et al. (2008) Chatterjee, S., Lahiri, P. & Li, H. (2008), ‘Parametric bootstrap approximation to the distribution of eblup and related prediction intervals in linear mixed models’, Annals of Statistics 36, 1221–1245.
  * Costa et al. (2003) Costa, A., Satorra, A. & Ventura, E. (2003), ‘Using composite estimators to improve both domain and total area estimation’, Technical report .
  * Das et al. (2004) Das, K., Jiang, J. & Rao, J. (2004), ‘Mean squared error of empirical predictor’, The Annals of Statistics 32(2), 818–840.
  * Datta & Lahiri (2000) Datta, G. S. & Lahiri, P. (2000), ‘A unified measure of uncertainty of estimated best linear unbiased predictors in small area estimation problems’, Statistica Sinica pp. 613–627.
  * Fay & Herriot (1979) Fay, R. & Herriot, R. (1979), ‘Estimates of income from small places: an application of James-Stein procedures to census data’, Journal of the American Statistical Association 74, 269–277.
  * Hall & Maiti (2006a) Hall, P. & Maiti, T. (2006a), ‘On parametric bootstrap methods for small area prediction’, Journal of the Royal Statistical Society Series B 68, 221–238.
  * Hall & Maiti (2006b) Hall, P. & Maiti, T. (2006b), ‘Nonparametric estimation of mean-squared prediction error in nested-error regression models’, Annals of Statistics 34, 1733–1750.
  * Haslett & Welsh (2019) Haslett, S. & Welsh, A. (2019), ‘EBLUPs: Empirical best linear unbiased predictors’, Wiley StatsRef: Statistics Reference Online .
  * Jiang & Lahiri (2006) Jiang, J. & Lahiri, P. (2006), ‘Estimation of finite population domain means: A model-assisted empirical best prediction approach’, Journal of the American Statistical Association 101(473), 301–311.
  * Jiang et al. (2002) Jiang, J., Lahiri, P. & Wan, S.-M. (2002), ‘A unified jackknife theory for empirical best prediction with M-estimation’, The Annals of Statistics 30, 1782–1810.
  * Jiang et al. (2011) Jiang, J., Nguyen, T. & Rao, J. S. (2011), ‘Best predictive small area estimation’, Journal of the American Statistical Association 106(494), 732–745.
  * Kackar & Harville (1981) Kackar, R. N. & Harville, D. A. (1981), ‘Unbiasedness of two-stage estimation and prediction procedures for mixed linear models’, Communications in Statistics-theory and Methods 10, 1249–1261.
  * Lahiri & Rao (1995) Lahiri, P. & Rao, J. (1995), ‘Robust estimation of mean squared error of small area estimators’, Journal of the American Statistical Association 90(430), 758–766.
  * Lehtonen & Veijanen (2009) Lehtonen, R. & Veijanen, A. (2009), Design-based methods of estimation for domains and small areas, in ‘Handbook of statistics’, Vol. 29, Elsevier, pp. 219–249.
  * Lyu & Welsh (2021) Lyu, Z. & Welsh, A. (2021), ‘Asymptotics for EBLUPs: Nested error regression models’, Journal of the American Statistical Association pp. 1–15.
  * Lyu & Welsh (2022) Lyu, Z. & Welsh, A. (2022), ‘Increasing cluster size asymptotics for nested error regression models’, Journal of Statistical Planning and Inference 217, 52–68.
  * Pfeffermann (2013) Pfeffermann, D. (2013), ‘New important developments in small area estimation’, Statistical Science 28(1), 40–68.
  * Prasad & Rao (1990) Prasad, N. & Rao, J. (1990), ‘The estimation of the mean squared error of small-area estimators’, Journal of the American Statistical Association 85, 163–171.
  * Pratesi (2016) Pratesi, M. (2016), ‘Analysis of poverty data by small area estimation’.
  * Rao (2005) Rao, J. (2005), ‘Inferential issues in small area estimation: some new developments’, Statistics in Transition 7(3), 513–526.
  * Rao (2008) Rao, J. (2008), ‘Some methods for small area estimation’, Rivista Internazionale di Scienze Sociali pp. 387–406.
  * Rao & Molina (2015) Rao, J. & Molina, I. (2015), Small Area Estimation, New Jersey: John Wiley & Sons.
  * Saei & Chambers (2003a) Saei, A. & Chambers, R. (2003a), ‘Small area estimation under linear and generalized linear mixed models with time and area effects’, Southampton Statistical Sciences Research Institute Methodology Working Paper M03/15, University of Southampton .
  * Saei & Chambers (2003b) Saei, A. & Chambers, R. (2003b), ‘Small area estimation: A review of methods based on the application of mixed models’, Southampton Statistical Sciences Research Institute Methodology Working Paper M03/16, University of Southampton .
  * Saei & Chambers (2005) Saei, A. & Chambers, R. (2005), ‘Empirical best linear unbiased prediction for out of sample areas’, Southampton Statistical Sciences Research Institute Methodology Working Paper M05/03, University of Southampton .
  * Sinha & Rao (2009) Sinha, S. K. & Rao, J. (2009), ‘Robust small area estimation’, Canadian Journal of Statistics 37(3), 381–399.
  * Torabi & Jiang (2020) Torabi, M. & Jiang, J. (2020), ‘Estimation of mean squared prediction error of empirically spatial predictor of small area means under a linear mixed model’, Journal of Statistical Planning and Inference 208, 82–93.
  * Torabi & Rao (2013) Torabi, M. & Rao, J. (2013), ‘Estimation of mean squared error of model-based estimators of small area means under a nested error linear regression model’, Journal of Multivariate Analysis 117, 76–87.
  * Tzavidis et al. (2010) Tzavidis, N., Marchetti, S. & Chambers, R. (2010), ‘Robust estimation of small-area means and quantiles’, Australian & New Zealand Journal of Statistics 52(2), 167–186.
  * Yoon & Welsh (2020) Yoon, H.-J. & Welsh, A. H. (2020), ‘On the effect of ignoring correlation in the covariates when fitting linear mixed models’, Journal of Statistical Planning and Inference 204, 18–34.
  * Za̧dło (2009) Za̧dło, T. (2009), ‘On MSE of EBLUP’, Statistical Papers 50, 101–118.