Mean Square Error Estimation of Small Area Predictors by Use of Parametric and Nonparametric Bootstrap

Abstract

In this article, we propose and compare some old and new parametric and nonparametric bootstrap methods for MSE estimation in small area estimation, restricting to the case of the widely used Fay-Herriot model. The parametric method consists of generating parametrically a large number of area bootstrap samples from the model fitted to the original data, re-estimating the model parameters for each bootstrap sample and then estimating the separate components of the MSE. The use of double-bootstrap is also considered. The nonparametric method generates the samples by bootstrapping standardized residuals, estimated from the original sample data. The bootstrap procedures are compared to other methods proposed in the literature in a simulation study, which also examines the robustness of the various methods to non-normality of the model error terms. A design-based MSE estimator for the Fay-Herriot model-dependent predictor is also described and its performance is investigated in a separate simulation study.

AMS subject classification: 62F10, 62F40

Keywords

Design-based MSE EBLUP Fay-Herriot Jackknife Order of bias

1. Introduction

Over the last four decades, there is growing demand all over the world for reliable estimates of small area parameters such as means, counts, proportions or quantiles. The estimates are used for fund allocations, new social and health programs, and more generally, for short and long term planning. Small area estimates (SAE) are also used for testing, correcting and supplementing administrative records. Although commonly known as small area estimation’, the domains of study may consist of socio-demographic subgroups as defined, for example, by gender, age and race, or the intersection of such domains with geographical locations.

The problem of SAE is that the sample sizes in at least some of the domains of study are very small, and often there are no samples available for many or even most of these domains. In such cases, the direct estimates obtained from a survey are unreliable with unacceptable large variances, and no direct survey estimates can be computed for areas with no samples. SAE methodology addresses therefore the following two major problems:

1. How to obtain reliable estimates for each of the areas,

2. How to assess the error of the estimators (MSE, confidence intervals, etc.).

In the present article, we restrict to the popular Fay-Herriot^[7] model and consider the estimation of the MSE of the Empirical Best Linear Unbiased Predictor (EBLUP). Due to time limitation, we only consider areas with samples. The computation of reliable MSE estimators in SAE problems is complicated because the models in use and the small sample sizes within the areas require accounting for the contribution to the error resulting from estimating the model parameters. Several procedures have been proposed in the literature, some of which we consider and compare in a simulation study in the present article. Our main goal is to propose new parametric and nonparametric bootstrap procedures for MSE estimation with correct order of bias, which to the best of our knowledge have not been proposed in the literature in the context of SAE. The parametric method consists of generating parametrically a large number of area bootstrap samples from the model fitted to the original data, re-estimating the model parameters for each bootstrap sample and then estimating separately or jointly the components of the MSE. The double bootstrap procedure is also considered. The nonparametric method generates the samples by bootstrapping standardized residuals computed from the original sample data.

In Section 2, we define the model, the resulting predictors and their theoretical MSEs. In Section 3 we describe two, commonly used estimators of the variance of the random effects, which is a major component of the model. Section 4 contains our proposed parametric and nonparametric bootstrap MSE estimators. Other procedures for MSE estimation proposed in the literature, including estimation of the randomization MSE over all possible sample selections are described in Section 5. In Section 6, we report the results of a simulation study, which compares the MSE estimators considered in the article. We conclude with some brief comments in Section 7.

2. The Fay-Herriot Area Level Model, Estimators and MSE’s

This model is in broad use when the sample information is only available at the area level. It was used originally by Fay and Herriot^[7] (hereafter FH) for predicting the per-capita income in geographical areas of less than 500 residents.

Denote by $y_{i}$ , the direct sample estimator of the mean in area $i$ (based only on the sample from that area), and by $θ_{i}$ the corresponding true area mean. Let $m$ denote the number of areas with observations. The model assumes,

y_{i} = θ_{i} + e_{i}; θ_{i} = x_{i}^{'} β + u_{i}, i = 1, \dots, m

(2.1)

where $x_{i}$ is a $p \times 1$ column vector of known area level characteristics (covariates), $β = (β_{1}, ..., β_{p})^{'}$ is a fixed vector of regression coefficients, $e_{i}$ represents the sampling error, assumed to have zero mean and known design variance $V a r_{D} (e_{i}) = σ_{D i}^{2}$ , and $u_{i}$ is a random effect, assumed to have zero mean and variance $σ_{u}^{2}$ . It is assumed that $E (e_{i} u_{j}) = 0 \forall i, j$ . For known model parameters $(β, σ_{u}^{2})$ and under normality of the error terms ( $e_{i}, u_{i}$ ), the best predictor (minimum MSE) of $θ_{i}$ is,

{\hat{θ}}_{i}^{B P} = γ_{i} y_{i} + (1 - γ_{i}) x_{i}^{'} β = x_{i}^{'} β + γ_{i} (y_{i} - x_{i}^{'} β)

(2.2)

The predictor (2.2) is a composite estimator’ with weight $γ_{i} = σ_{u}^{2} / (σ_{D i}^{2} + σ_{u}^{2})$ , which determines how much weight is assigned to the direct estimator and how much to the synthetic part, $x_{i}^{'} β$ , depending on the corresponding error variances $σ_{D i}^{2}$ and $σ_{u}^{2}$ . The MSE of ${\hat{θ}}_{i}^{B P}$ under the model is,

M S E ({\hat{θ}}_{i}^{B P}) = E {({\hat{θ}}_{i}^{B P} - θ_{i})}^{2} = g_{1 i} (σ_{u}^{2}) = γ_{i} σ_{D i}^{2}

(2.3)

REMARK 1. Equation (2.3) holds also without the normality assumptions. The normality assumptions guarantee that the predictor (2.2) attains the minimum MSE.

When $σ_{u}^{2}$ is known but $β$ is unknown, the best linear unbiased predictor (BLUP) of $θ_{i}$ is obtained by replacing $β$ in (2.2) by the generalized least square estimator (GLS),

{\hat{β}}_{G L S} = {[\sum_{i = 1}^{m} \frac{1}{σ_{u}^{2} + σ_{D_{i}}^{2}} x_{i} x_{i}^{'}]}^{- 1} \sum_{i = 1}^{m} \frac{1}{σ_{u}^{2} + σ_{D_{i}}^{2}} x_{i} y_{i}

(2.4)

The BLUP is thus,

{\hat{θ}}_{i}^{B L U P} = γ_{i} y_{i} + (1 - γ_{i}) x_{i}^{'} {\hat{β}}_{G L S}

(2.5)

The MSE is,

M S E ({\hat{θ}}_{i}^{B L U P}) = E {({\hat{θ}}_{i}^{B L U P} - θ_{i})}^{2} = g_{1 i} (σ_{u}^{2}) + g_{2 i} (σ_{u}^{2})

(2.6)

where $g_{2 i} (σ_{u}^{2})$ represents the additional error resulting from estimating $β$ ;

g_{2 i} (σ_{u}^{2}) = {(1 - γ_{i})}^{2} x_{i}^{'} {[\sum_{i = 1}^{m} \frac{1}{σ_{u}^{2} + σ_{D_{i}}^{2}} x_{i} x_{i}^{'}]}^{- 1} x_{i} = {(1 - γ_{i})}^{2} x_{i}^{'} V a r ({\hat{β}}_{G L S}) x_{i}

(2.7)

REMARK 2. The BLUP property and the MSE expression (2.6) are valid without the normality assumptions of the error terms.

In practice, both $β$ and $σ_{u}^{2}$ are unknown and need to be estimated from the observed data. An empirical BLUP (EBLUP) is obtained by replacing $σ_{u}^{2}$ by an estimator ${\hat{σ}}_{u}^{2} (y)$ in the expression (2.5) of the BLUP, where $y = (y_{1}, ..., y_{m})^{'}$ . The predictor is,

{\hat{θ}}_{i}^{E B L U P} = {\hat{γ}}_{i} y_{i} + (1 - {\hat{γ}}_{i}) x_{i}^{'} {\hat{β}}_{i} ({\hat{σ}}_{u}^{2})

(2.8)

where ${\hat{γ}}_{i}$ and ${\hat{β}}_{i}$ are obtained from $γ_{i}$ and ${\hat{β}}_{G L S}$ by replacing $σ_{u}^{2}$ by ${\hat{σ}}_{u}^{2} (y)$ .

The question arising is how to estimate $σ_{u}^{2}$ and how to estimate the MSE of the resulting EBLUP defined by (2.8), to a correct order of bias.

3. Procedures Proposed in the Literature for Estimating

σ_{u}^{2}

In what follows we describe two procedures that we use in our simulation study:

a) Prasad and Rao^[20] estimator;

{\tilde{σ}}_{P R}^{2} = \frac{1}{(m - p)} [\sum_{i = 1}^{m} {(y_{i} - x_{i}^{'} {\hat{β}}_{O L S})}^{2} - \sum_{i = 1}^{m} σ_{D i}^{2} h_{i}]; h_{i} = (1 - x_{i}^{'} {[\sum_{i = 1}^{m} x_{i} x_{i}^{'}]}^{- 1} x_{i})

(3.1)

where ${\hat{β}}_{OLS} = {[\sum_{i = 1}^{m} x_{i} x_{i}^{'}]}^{- 1} \sum_{i = 1}^{m} x_{i} y_{i}$ . The estimator ${\tilde{σ}}_{P R}^{2}$ can be negative, and so ${\hat{σ}}_{P R}^{2} = \max ({\tilde{σ}}_{P R}^{2}, 0)$ is usually the estimator used in practice. The estimator ${\tilde{σ}}_{P R}^{2}$ (but not ${\hat{σ}}_{P R}^{2}$ ) is unbiased.

b) Fay and Herriot^[7] estimator;

Solve iteratively,

\frac{1}{(m - p)} \sum_{i = 1}^{m} \frac{1}{({\hat{σ}}_{u}^{2} + σ_{D i}^{2})} {[y_{i} - x_{i}^{'} \hat{β} ({\hat{σ}}_{n}^{2})]}^{2} = 1

(3.2)

where $\hat{β} ({\hat{σ}}_{u}^{2})$ is the GLS estimator (2.4), with $σ_{u}^{2}$ replaced by ${\hat{σ}}_{u}^{2}$ . Define the solution by ${\tilde{σ}}_{F H}^{2}$ and set ${\hat{σ}}_{F H}^{2} = \max ({\tilde{σ}}_{F H}^{2}, 0)$ . The rationale of (3.2) is that for ${\hat{σ}}_{u}^{2} = σ_{u}^{2}$ , the expectation of the left-hand side of (3.2) equals 1.

REMARK 3. The estimators defined by (3.1) and (3.2) satisfy (i) $({\hat{σ}}_{u}^{2} - σ_{u}^{2}) = O_{p} (m^{- 0.5}),$ (ii) they are even functions of $y$ such that ${\hat{σ}}_{u}^{2} (y) = {\hat{σ}}_{u}^{2} (- y)$ ; and (iii) they are translation invariant; ${\hat{σ}}_{u}^{2} (y) = {\hat{σ}}_{u}^{2} (y + X d)$ for any vector $d \in R^{p}$ and all $y$ , where $X = [x_{1}, \dots, x_{m}]$ . Under these conditions, the resulting EBLUP predictors remain unbiased.

REMARK 4. Pfeffermann and Nathan^[16] proposed a similar estimator to ${\hat{σ}}_{F H}^{2}$ in the context of regression analysis from a cluster sample with random cluster slopes, and showed some other desirable properties of this estimator.

REMARK 5. Several other procedures have been proposed in the literature for estimating the variance $σ_{u}^{2}$ . Datta and Lahiri^[4] derive maximum likelihood (MLE)- and residual maximum likelihood (REML) estimators for a general mixed linear model under normality of the error terms, which satisfy the regularity conditions in Remark 3 above. These estimators can likewise be negative, particularly with small m. To deal with this problem, Li and Lahiri^[14] propose adjustments to the MLE and REML estimators that produce strictly positive estimates of $σ_{u}^{2}$ . These adjusted estimators also satisfy the regularity conditions in Remark 3. We do not consider further the estimators mentioned in this remark in the present paper. The use of them requires normality of the error terms, but they remain consistent under general conditions, even without the normality assumption (Jiang^[10]). Some of the Jackknife procedures for MSE estimation described in Section 5.2, use a Jackknife estimator of $σ_{u}^{2}$ .

4. MSE Estimation of the EBLUP by Bootstrap

4.1 MSE Decomposition

The EBLUP is defined by (2.8). The prediction error can be decomposed as,

({\hat{θ}}_{i}^{E B L U P} - θ_{i}) = ({\hat{θ}}_{i}^{B L U P} - θ_{i}) + ({\hat{θ}}_{i}^{E B L U P} - {\hat{θ}}_{i}^{B L U P})

(4.1)

where ${\hat{θ}}_{i}^{B L U P}$ is defined by (2.5). Hence, by (2.6),

M S E ({\hat{θ}}_{i}^{E B L U P}) = E {({\hat{θ}}_{i}^{E B L U P} - θ_{i})}^{2} = g_{1 i} (σ_{u}^{2}) + g_{2 i} (σ_{u}^{2}) + E {({\hat{θ}}_{i}^{E B L U P} - {\hat{θ}}_{i}^{B L U P})}^{2} + 2 E ({\hat{θ}}_{i}^{B L U P} - θ_{i}) ({\hat{θ}}_{i}^{E B L U P} - {\hat{θ}}_{i}^{B L U P})

(4.2)

Under normality of the model error terms (sampling errors and random effects), and for estimators ${\hat{σ}}_{u}^{2}$ satisfying the conditions (ii) and (iii) in Remark 3, the cross-product expectation in 4.2 vanishes. $C P E = E ({\hat{θ}}_{i}^{B L U P} - θ_{i}) ({\hat{θ}}_{i}^{E B L U P} - {\hat{θ}}_{i}^{B L U P}) = 0$ (Harville^[8]). However, for other distributions of the model error terms, the cross-product expectation may not vanish and it is of similar magnitude as the second and third terms in the right hand side of (4.2). Lahiri and Rao^[13] developed a second order approximation for the cross-product expectation for the case where ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$ , which only requires that the sampling errors are normally distributed. The approximation involves the fourth moment of the distribution of the random effects.

4.2 Parametric Bootstrap Method for MSE Estimation of the EBLUP

The method consists of the following steps:

P1. For $b = 1, ..., B$ (B large), generate independently normal random effects $u^{b} = (u_{1}^{b}, ... u_{m}^{b})^{'}$ and normal sampling errors $e^{b} = (e_{1}^{b}, ... e_{m}^{b})^{'}$ , and hence bootstrap direct estimators $y^{b} = (y_{1}^{b}, ... y_{m}^{b})^{'}$ from the F-H model (2.1), with hyper-parameters equal to $σ_{D_{i}}^{2}$ , ${\hat{σ}}_{u}^{2} (y)$ and $\hat{β} [y; {\hat{σ}}_{u}^{2} (y)]$ , where $y$ defines the original (parent) sample.

P2. Re-estimate $σ_{u}^{2}$ and $β$ for each of the bootstrap samples, using the same method as used for the original sample, yielding the estimators ${\hat{σ}}_{u}^{2} (y^{b})$ , $\hat{β} [y^{b}; {\hat{σ}}_{u}^{2} (y^{b})]$ , and also $\hat{β} [y^{b}; {\hat{σ}}_{u}^{2} (y)]$ .

P3. Estimate the MSE of the EBLUP as,

M \hat{S} E_{1}^{P B} ({\hat{θ}}_{i}^{E B L U P}) = 2 [g_{1 i} ({\hat{σ}}_{u}^{2} (y)) + g_{2 i} ({\hat{σ}}_{u}^{2} (y))] - {\bar{g}}_{1 i}^{P B} - {\bar{g}}_{2 i}^{P B} + p u c^{P B}

(4.3)

where $p u c^{P B} = B^{- 1} \sum_{b = 1}^{B} {{\hat{θ}}_{i} [y^{b}; {\hat{σ}}_{u}^{2} (y^{b}), \hat{β} (y^{b}; {\hat{σ}}_{u}^{2} (y^{b}))] - {\hat{θ}}_{i} [y^{b}; {\hat{σ}}_{u}^{2} (y), \hat{β} (y^{b}; {\hat{σ}}_{u}^{2} (y))]}^{2}$ and ${\bar{g}}_{t i}^{P B} = B^{- 1} \sum_{b = 1}^{B} g_{t i} ({\hat{σ}}_{u}^{2} (y^{b})); t = 1, 2$ .

The term $p u c^{P B}$ estimates the contribution to the MSE from the parameter uncertainty, as defined by the third term on the right side of (4.2).

Using similar arguments to Pfeffermann and Tiller^[17], it follows that under mild regularity conditions, the MSE estimator (4.3) has bias of order $O (m^{- 2})$ .

REMARK 6. Pfeffermann and Tiller^[17] consider MSE estimation of EBLUP state predictors in the context of state-space models, which contain the FH model as a simple special case.

REMARK 7. Butar and Lahiri^[1] likewise developed the MSE estimator (4.3) although in a different way, and showed that it has bias of order $o (m^{- 1})$ .

The MSE estimator defined by (4.3) assumes that the model error terms are normally distributed and hence that the cross product expectation in (4.2) is zero. When this is not the case, the cross-product expectation may not vanish. Assuming that the true distributions of the random errors are known, one needs to generate the bootstrap samples in Step P1 above by sampling the error terms from their respective distributions, and adding twice the following expression to the estimator (4.3):

\begin{array}{l} C \hat{P} E^{P B} = B^{- 1} \sum_{b = 1}^{B} {{\hat{θ}}_{i} [y^{b}; {\hat{σ}}_{u}^{2} (y^{b}), \hat{β} (y^{b}; {\hat{σ}}_{u}^{2} (y^{b}))] - {\hat{θ}}_{i} [y^{b}; {\hat{σ}}_{u}^{2} (y), \hat{β} (y^{b}; {\hat{σ}}_{u}^{2} (y))]} \\ \times {{\hat{θ}}_{i} [y^{b}; {\hat{σ}}_{u}^{2} (y), \hat{β} (y^{b}; {\hat{σ}}_{u}^{2} (y))] - θ_{i}^{b}}, \end{array}

(4.4)

where $θ_{i}^{b} = x_{i}^{'} \hat{β} (y; {\hat{σ}}_{u}^{2} (y)) + u_{i}^{b}$ is the true’ area mean generated for area $i$ in bootstrap sample $b$ .

An alternative parametric bootstrap estimator, also resulting from Pfeffermann and Tiller^[17], is obtained by replacing (4.3) by

M \hat{S} E_{2}^{P B} ({\hat{θ}}_{i}^{E B L U P}) = [g_{1 i} ({\hat{σ}}_{u}^{2} (y)) + g_{2 i} ({\hat{σ}}_{u}^{2} (y))] - {\bar{g}}_{1 i}^{P B} - {\bar{g}}_{2 i}^{P B} + m s e^{P B}

(4.5)

where $m s e^{P B} = B^{- 1} \sum_{b = 1}^{B} {{\hat{θ}}_{i} [y^{b}; {\hat{σ}}_{u}^{2} (y^{b}), \hat{β} (y^{b}; {\hat{σ}}_{u}^{2} (y^{b}))] - θ_{i}^{b}}^{2}$ is the MSE of the EBLUP under the bootstrap model. It is a naive’ MSE estimator because it ignores the bias resulting from generating the bootstrap samples with a sample estimator ${\hat{σ}}_{u}^{2}$ , rather than with the true value $σ_{u}^{2}$ . For distributions such that the cross-product expectation in (4.2) is of order $O (m^{- 1})$ , the MSE estimator (4.5) has bias of order $O (m^{- 2})$ .

The estimator (4.5) is equivalent asymptotically to the estimator (4.3), but it has the potential advantage of robustness against non-normal distributions of the model error terms. To see this, denote by $E_{*}$ the expectation with respect to the bootstrap model, that is, when generating the area direct estimators with parameters ${\hat{σ}}_{u}^{2} (y), \hat{β} (y; {\hat{σ}}_{u}^{2} (y))$ . Then, in analogy to Eq. (4.2),

E_{*} {({\hat{θ}}_{i}^{b, E B L U P} - θ_{i}^{b})}^{2} = g_{1 i} ({\hat{σ}}_{u}^{2} (y)) + g_{2 i} ({\hat{σ}}_{u}^{2} (y)) + E_{*} {({\hat{θ}}_{i}^{b, E B L U P} - {\hat{θ}}_{i}^{b, B L U P})}^{2} + 2 E_{*} ({\hat{θ}}_{i}^{b, B L U P} - θ_{i}^{b}) ({\hat{θ}}_{i}^{b, E B L U P} - {\hat{θ}}_{i}^{b, B L U P})

(4.6)

Thus, the expression $E {({\hat{θ}}_{i}^{E B L U P} - {\hat{θ}}_{i}^{B L U P})}^{2} + 2 E ({\hat{θ}}_{i}^{B L U P} - θ_{i}) ({\hat{θ}}_{i}^{E B L U P} - {\hat{θ}}_{i}^{B L U P})$ in (4.2) can be estimated by $E_{*} {({\hat{θ}}_{i}^{E B L U P} - θ_{i})}^{2} - g_{1} ({\hat{σ}}_{u}^{2}) - {\hat{g}}_{2} (σ_{u}^{2}) = m s e^{P B} - g_{1} ({\hat{σ}}_{u}^{2}) - {\hat{g}}_{2} (σ_{u}^{2})$ . The estimator (4.5) is obtained by adding the last expression to the bias reduced estimator $2 [g_{1} ({\hat{σ}}_{u}^{2}) + g_{2} ({\hat{σ}}_{u}^{2})] - {\overset{࿽}{g}}_{1 i}^{P B} - {\overset{࿽}{g}}_{2 i}^{P B}$ of $[g_{1} (σ_{u}^{2}) + g_{2} (σ_{u}^{2})]$ . For distributions such that the cross-product term in (4.2) is of order $O (1 / m)$ , the MSE estimator (4.5) has bias of order $O (1 / m^{2})$ .

4.3 Nonparametric Bootstrap Method for MSE Estimation of the EBLUP

For nonparametric bootstrap, we propose using the original estimates of $σ_{u}^{2}$ and $β$ in order to generate bootstrap replications of estimated standardized combined error terms of $(u_{t} + e_{t})$ . The method consists of the following steps:

NP1. Calculate the $m$ estimated standardized residuals

{\hat{r}}_{i} = [y_{i} - x_{i}^{'} \hat{β} (y; {\hat{σ}}_{u}^{2} (y))] / f_{i}^{1 / 2}; f_{i} = ({\hat{σ}}_{u}^{2} (y) + σ_{D_{i}}^{2}) - x_{i}^{'} {[\sum_{i = 1}^{m} \frac{1}{{\hat{σ}}_{u}^{2} (y) + σ_{D_{i}}^{2}} x_{i} x_{i}^{'}]}^{- 1} x_{i}

(4.7)

Note: $f_{i} = V a r [y_{i} - x_{i}^{'} \hat{β} (y; {\hat{σ}}_{u}^{2} (y))]$ under the FH model, with ${\hat{σ}}_{u}^{2} (y), σ_{D_{i}}^{2}$ and $\hat{β} (y; {\hat{σ}}_{u}^{2} (y))$ as the true’ model parameters.

NP2. Sample a large number B of sets of standardized residuals $r^{b} = (r_{1}^{b}, ..., r_{m}^{b}), b = 1, ..., B$ , where each set is a simple random sample with replacement of size $m$ from the standardized residuals ${\hat{r}}_{i}, i = 1, ..., m$ defined by (4.7).

NP3. Calculate the bootstrap direct estimators,

y_{i}^{b} = r_{i}^{b} {(f_{i})}^{1 / 2} + x_{i}^{'} \hat{β} (y; {\hat{σ}}_{u}^{2} (y)); i = 1, ..., m, b = 1, ..., B

(4.8)

Note: $V a r [r_{i}^{b} {(f_{i})}^{1 / 2}] = f_{i}$ , the true variance of the estimated residual term in area $i$ , under the setup above.

NP4. Re-estimate the parameters $σ_{u}^{2}$ and $β$ for each of the bootstrap samples using the same method as used for the original sample, yielding the estimators, ${\hat{σ}}_{u}^{2} (y^{b}), \hat{β} (y^{b}; {\hat{σ}}_{u}^{2} (y^{b}))$ , and $\hat{β} (y^{b}; {\hat{σ}}_{u}^{2} (y))$ . Predict, ${\hat{θ}}_{i} [y^{b}; {\hat{σ}}_{u}^{2} (y^{b}), \hat{β} (y^{b}; {\hat{σ}}_{u}^{2} (y^{b}))]$ and ${\hat{θ}}_{i} [y^{b}; {\hat{σ}}_{u}^{2} (y), \hat{β} (y^{b}; {\hat{σ}}_{u}^{2} (y))]$ .

NP5. Estimate the MSE of the EBLUP as,

M \hat{S} E^{N P B} ({\hat{θ}}_{i}^{E B L U P}) = 2 [g_{1 i} ({\hat{σ}}_{u}^{2} (y)) + g_{2 i} ({\hat{σ}}_{u}^{2} (y))] - {\bar{g}}_{1 i}^{N P B} - {\bar{g}}_{2 i}^{N P B} + p u c

(4.9)

where $p u c^{N P B}$ , ${\bar{g}}_{1 i}^{N P B}$ and ${\bar{g}}_{2 i}^{N P B}$ are defined similarly to in (4.3).

REMARK 8. The estimator (4.9) is essentially the same as (4.3), but based on nonparametric bootstrap. Notice, however, that by bootstrapping the estimated standardized residuals, it is no longer possible to generate true’ bootstrap area means $θ_{i}^{b}$ and hence to compute a nonparametric MSE estimator that is equivalent to (4.5). Also, notice that like the parametric estimator (4.3), the estimator (4.9) assumes that the cross-product expectation in (4.2) is zero, which does not necessarily hold if the true model error terms are not normally distributed.

REMARK 9. In a discussion to an article by Jiang and Lahiri^[12], the late Professor Peter Hall reckons that “small area methods are motivated when data are scarce, and it is exactly in such cases that informed parametric techniques can enjoy statistical advantages over their more adaptive nonparametric cousins. Parametric bootstrap methods therefore have an important role to play’’. Our simulation results in Section 6 do not backup this statement, at least for small m $m$ .

5. EBLUP MSE Estimators Proposed in the Literature

5.1 Estimators Based on Taylor Approximations

Prasad and Rao^[20] show that under normality of the model error terms, the MSE of the EBLUP computed with an estimator ${\hat{σ}}_{u}^{2}$ , can be approximated up to terms of order $o (1 / m)$ as,

M S E [{\hat{θ}}_{i}^{E B L U P} (σ_{u}^{2})] = g_{1 i} (σ_{u}^{2}) + g_{2 i} (σ_{u}^{2}) + g_{3 i} (σ_{u}^{2}) V a r ({\hat{σ}}_{u}^{2})

(5.1)

where $g_{3 i} (σ_{u}^{2}) = σ_{D i}^{4} {(σ_{u}^{2} + σ_{D i}^{2})}^{- 3}$ . For the case where ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$ (Equation 3.1), the authors develop the following estimator for the MSE approximation (5.1), with bias of order $o (1 / m)$ .

M \hat{S} E [{\hat{θ}}_{i}^{E B L U P} ({\hat{σ}}_{P R}^{2})] = g_{1 i} ({\hat{σ}}_{P R}^{2}) + g_{2 i} ({\hat{σ}}_{P R}^{2}) + 2 g_{3 i} ({\hat{σ}}_{P R}^{2}) {\hat{V}}_{P R}

(5.2)

where ${\hat{V}}_{P R} = V \hat{a} r ({\hat{σ}}_{P R}^{2}) = 2 m^{- 2} \sum_{i = 1}^{m} {({\hat{σ}}_{P R}^{2} + σ_{D i}^{2})}^{2}$ .

REMARK 10. Lahiri and Rao^[13] show that the estimator (5.2) is robust to non-normality of the distribution of the model random effects.

Datta et al.^[5] consider the case where $σ_{u}^{2}$ is estimated by the FH method ( ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$ , Equation 3.2). The authors derive the following MSE estimator, with bias correct to the order $o (1 / m),$ ,

M \hat{S} E [{\hat{θ}}_{i}^{E B L U P} ({\hat{σ}}_{F H}^{2})] = g_{1 i} ({\hat{σ}}_{F H}^{2}) + g_{2 i} ({\hat{σ}}_{F H}^{2}) + 2 g_{3 i} ({\hat{σ}}_{F H}^{2}) {\hat{V}}_{F H} - g_{4 i} ({\hat{σ}}_{F H}^{2})

(5.3)

where ${\hat{V}}_{F H} = V \hat{a} r ({\hat{σ}}_{F H}^{2}) = 2 m {[\sum_{i = 1}^{m} {({\hat{σ}}_{F H}^{2} + σ_{D i}^{2})}^{- 1}]}^{- 2}$ and

g_{4 i} ({\hat{σ}}_{F H}^{2}) = 2 {[1 - γ_{i} ({\hat{σ}}_{F H}^{2})]}^{2} [m \sum_{i = 1}^{m} δ_{i}^{- 2} - {(\sum_{i = 1}^{m} δ_{i}^{- 1})}^{2}] {(\sum_{i = 1}^{m} δ_{i}^{- 1})}^{- 3}; δ_{i} = σ_{D i}^{2} + {\hat{σ}}_{F H}^{2}

(5.4)

5.2 Estimators Based on Jackknife Resampling

An alternative approach for estimating the MSE of the EBLUP is the use of Jackknife procedures. Jiang et al.^[11] (hereafter JLW), develop a unified theory for estimation of the MSE of empirical best predictors (EBP) under a general class of mixed models, which includes the FH model as a simple special case. Recall that for the FH model with normal error terms, ${\hat{θ}}_{i}^{E B L U P} = {\hat{θ}}_{i}^{E B P}$ , where ${\hat{θ}}_{i}^{E B P}$ is the empirical best predictor of $θ_{i}$ , with ${\hat{θ}}_{i}^{B P}$ defined by (2.2) (assuming known $β$ and $σ_{u}^{2}$ ).

Similarly to (4.1), $({\hat{θ}}_{i}^{E B L U P} - θ_{i}) = ({\hat{θ}}_{i}^{B P} - θ_{i}) + ({\hat{θ}}_{i}^{E B L U P} - {\hat{θ}}_{i}^{B P})$ . Hence, by (2.3),

$M S E ({\hat{θ}}_{i}^{E B L U P}) = g_{1 i} (σ_{u}^{2}) + E {({\hat{θ}}_{i}^{E B L U P} - {\hat{θ}}_{i}^{B P})}^{2}$ . JLW proposed the following Jackknife MSE estimator:

$M S E ({\hat{θ}}_{i, J L W}^{E B L U P}) = {\hat{M}}_{1 i} + {\hat{M}}_{2 i}$

\begin{array}{l} {\hat{M}}_{1 i} = g_{1 i} ({\hat{σ}}_{u}^{2} (y)) - \frac{m - 1}{m} \sum_{j = 1}^{m} [g_{1 i} ({\hat{σ}}_{u}^{2} (y_{- j})) - g_{1 i} ({\hat{σ}}_{u}^{2} (y))] \\ {\hat{M}}_{2 i} = \frac{m - 1}{m} \sum_{j = 1}^{m} {\{{\hat{θ}}_{i} [y_{i}; {\hat{σ}}_{u}^{2} (y_{- j}), \hat{β} (y_{- j}; {\hat{σ}}_{u}^{2} (y_{- j}))] - {\hat{θ}}_{i} [(y_{i}; {\hat{σ}}_{u}^{2} (y)), \hat{β} (y; {\hat{σ}}_{u}^{2} (y))]\}}^{2} \end{array}

(5.5)

In (5.5), $y_{- j}$ is the vector of observations without the j^th area direct estimator, $\hat{β} (y_{- j}; σ_{u}^{2} (y_{- j}))$ is the GLS estimator of $β$ based on $y_{- j}$ and ${\hat{σ}}_{u}^{2} (y_{- j})$ , and ${\hat{θ}}_{i} (y_{i}; {\hat{σ}}_{u}^{2} (y_{- j}), \hat{β} (y_{- j}; {\hat{σ}}_{u}^{2} (y_{- j}))$ is the EBLUP predictor of $θ_{i}$ based on $y_{- j}$ .

REMARK 11. Lohr and Rao^[15] propose a modification of the estimator (5.5), which is simpler computationally and estimates the conditional MSE, $E [{({\hat{θ}}_{i}^{E B L U P} - θ_{i})}^{2} | y_{i}]$ . Denote $ψ = (σ_{u}^{2}, β)$ and $q_{i} (ψ, y_{i}) = V a r (θ_{i} | y_{i}; ψ)$ . The modification consists of replacing ${\hat{M}}_{1 i}$ in (5.5) by ${\hat{M}}_{1 i, c} = q_{i} (\hat{ψ}, y_{i}) - \sum_{j \neq i}^{m} [q_{i} ({\hat{ψ}}_{(- j)}, y_{i}) - q_{i} (\hat{ψ}, y_{i})]$ . When estitmating $ψ$ by MLE, the modified estimator ${\hat{λ}}_{i, c}^{J K} = {\hat{M}}_{1 i, c} + {\hat{M}}_{2 i}$ has bias of order $o_{p} (1 / m)$ in estimating the conditional MSE and bias of order $o (1 / m)$ in estimating the unconditional MSE. Lohr and Rao^[15] note that for the Fay-Herriot model, the estimator ${\hat{λ}}_{i, c}^{J K}$ is approximately the same as the estimator (5.5).

Chen and Lahiri^[2], develop the following Jackknife estimator:

\begin{array}{l} M \hat{S} E_{1} ({\hat{θ}}_{i J K, C L}^{E B L P P}) = g_{1 i} ({\hat{σ}}_{u}^{2} (y)) + g_{2 i} ({\hat{σ}}_{u}^{2} (y)) \\ - \frac{m - 1}{m} \sum_{j = 1}^{m} [g_{1 i} ({\hat{σ}}_{u}^{2} (y_{- j})) + g_{2 i} ({\hat{σ}}_{u}^{2} (y_{- j})) - g_{1 i} ({\hat{σ}}_{u}^{2} (y)) - g_{2 i} ({\hat{σ}}_{u}^{2} (y))] \\ + \frac{m - 1}{m} \sum_{j = 1}^{m} {\{{\hat{θ}}_{i} [y_{i}; {\hat{σ}}_{u}^{2} (y_{- j}), \hat{β} (y_{- j}; {\hat{σ}}_{u}^{2} (y_{- j}))] - {\hat{θ}}_{i} [(y_{i}; {\hat{σ}}_{u}^{2} (y)), \hat{β} (y; {\hat{σ}}_{u}^{2} (y))]\}}^{2} \end{array}

(5.6)

For the case where the estimator $M \hat{S} E_{1} ({\hat{θ}}_{i J K, C L}^{E B L U P})$ is negative (may happen with small $m$ ), the authors propose to replace the expression in the second row of (5,6) by $σ_{D i}^{4} {[σ_{D i}^{2} + {\hat{σ}}_{u}^{2} (y)]}^{- 3} {\hat{V}}^{J K} ({\hat{σ}}_{u}^{2})$ , where ${\hat{V}}^{J K} ({\hat{σ}}_{u}^{2}) = \frac{m - 1}{m} \sum_{j = 1}^{m} [{\hat{σ}}_{u}^{2} (y_{- j}) - {\hat{σ}}_{u}^{2} (y)]^{2}$ is the Jackknife estimator of $V a r ({\hat{σ}}_{u}^{2})$ .

Chen and Lahiri^[3], approximate additionally the last expression of (5.6) by $σ_{D i}^{4} {[σ_{D i}^{2} + {\hat{σ}}_{u}^{2} (y)]}^{- 4} {[y_{i} - x_{i}^{'} \hat{β} (y; {\hat{σ}}_{u}^{2} (y))]}^{2} {\hat{V}}^{J K} ({\hat{σ}}_{u}^{2})$ . Thus, the MSE estimator proposed by Chen and Lahiri^[3] is,

M \hat{S} E_{2} ({\hat{θ}}_{i J K, C L}^{E B L U P}) = g_{1 i} ({\hat{σ}}_{u}^{2} (y)) + g_{2 i} ({\hat{σ}}_{u}^{2} (y)) + \frac{σ_{D i}^{4}}{{[σ_{D i}^{2} + {\hat{σ}}_{u}^{2} (y)]}^{3}} {\hat{V}}^{J K} ({\hat{σ}}_{u}^{2}) + \frac{σ_{D i}^{4}}{{[σ_{D i}^{2} + {\hat{σ}}_{u}^{2} (y)]}^{4}} {[y_{i} - x_{i}^{'} \hat{β} (y; {\hat{σ}}_{u}^{2} (y))]}^{2} {\hat{V}}^{J K} ({\hat{σ}}_{u}^{2}) .

(5.7)

REMARK 12. Under normality of the error terms $(u_{i}, e_{i})$ , the three Jackknife estimators considered above have bias of order $o (1 / m)$ in estimating the unconditional MSE over the joint distribution of the random effects and the sampling errors.

5.3 Estimator Based on Double Parametric Bootstrap

Hall and Maiti^[9], propose estimating the MSE by use of double-bootstrap. For the FH model (2.1), the procedure consists of the following steps, where we denote by $\hat{ψ} = [{\hat{σ}}_{u}^{2}, \hat{β} ({\hat{σ}}_{u}^{2})]$ the estimators obtained from the original sample.

DB1. Generate a new population of area means from the model (2.1), with parameters $\hat{ψ}$ . Generate a sample and compute the EBLUP based on newly estimated parameters. The new population uses the same covariates as in the original population. Repeat the process independently $B_{1}$ times, with $B_{1}$ large. Denote by $θ_{i, b} (\hat{ψ})$ and ${\hat{θ}}_{i, b}^{E B L U P} ({\hat{ψ}}_{b})$ the true’ mean and corresponding EBLUP for population and sample $b$ , $b = 1, ..., B_{1}$ . Compute the 1^st step bootstrap MSE estimator (same as $m s e^{P B}$ in Eq. 4.5),

M \hat{S} E_{1}^{B S} [{\hat{θ}}_{i}^{E B L U P}] = \frac{1}{B_{1}} \sum_{b = 1}^{B_{1}} {[{\hat{θ}}_{i, b}^{(E B L U P)} ({\hat{ψ}}_{b}) - θ_{i, b} (\hat{ψ})]}^{2}

(5.8)

DB2. For each sample $b$ drawn in Step 1, repeat the computations of Step 1 $B_{2}$ times with $B_{2}$ sufficiently large, yielding new true’ means $θ_{i, b, c} ({\hat{ψ}}_{b})$ and EBLUPs ${\hat{θ}}_{i, b, c}^{E B L U P} ({\hat{ψ}}_{b, c}),$ $b = 1, ..., B_{1}; c = 1, ..., B_{2}$ . Compute the second-step bootstrap MSE estimator,

M \hat{S} E_{2}^{B S} [{\hat{θ}}_{i}^{E B L U P}] = \frac{1}{B_{1}} \sum_{b}^{B_{1}} \frac{1}{B_{2}} \sum_{c = 1}^{B_{2}} {[{\hat{θ}}_{i, b, c}^{(E B L U P)} ({\hat{ψ}}_{b, c}) - θ_{i, b, c} ({\hat{ψ}}_{b})]}^{2}

(5.9)

Denote ${\hat{λ}}_{i, 1}^{B S} = M \hat{S} E_{1}^{B S} ({\hat{θ}}_{i}^{E B L U P})$ , ${\hat{λ}}_{i, 2}^{B S} = M \hat{S} E_{2}^{B S} ({\hat{θ}}_{i}^{E B L U P})$ . The double-bootstrap MSE estimator is obtained by computing bias corrected estimators. For example,

{\hat{λ}}_{i}^{D B S} = \{\begin{matrix} {\hat{λ}}_{i, 1}^{B S} + ({\hat{λ}}_{i, 1}^{B S} - {\hat{λ}}_{i, 2}^{B S}), if {\hat{λ}}_{i, 1}^{B S} \geq {\hat{λ}}_{i, 2}^{B S} \\ {\hat{λ}}_{i, 1}^{B S} \exp [({\hat{λ}}_{i, 1}^{B S} - {\hat{λ}}_{i, 2}^{B S}) / {\hat{λ}}_{i, 2}^{B S}], if {\hat{λ}}_{i, 1}^{B S} < {\hat{λ}}_{i, 2}^{B S} \end{matrix}

(5.10)

REMARK 13. The 1^st step bootstrap estimator (5.8) has bias of order $O (1 / m)$ . The double-bootstrap estimator (5.10) has bias of order $o (1 / m)$ under some mild regularity conditions. However, the use of (5.10) may inflate the variance of the estimator and hence the MSE. To deal with this problem, Hall and Maiti^[9] propose using instead an estimator with a lower bias reduction but smaller MSE. See the article for details.

REMARK 14. The computation of the double bootstrap estimator is very computing intensive when applied with large $B_{2}$ . In the simulation study in Section 6, we follow Erciulescu and Fuller^[6] and set $B_{2} = 1$ , which yields similar results to those obtained with large values of $B_{2}$ .

5.4 Estimation of Randomization Mean Square Error

All the MSE estimators considered so far are model dependent, in our case for the FH model, accounting for all sources of variation. This implies that the target area means are viewed as random, which is different from classical survey sampling theory under which the finite population values, and hence the area means, or other parameters of interest are considered as fixed values. However, users of sample survey estimates are used to measures of error such as MSE, which only account for the variability originating from the randomness of the sample selection (known as the randomization distribution), that is, the MSE over all possible sample selections from the target finite population, with the population values of the survey variables held fixed. We refer to this MSE as the design-based MSE, denoted hereafter, DMSE.

Pfeffermann and Ben-Hur^[19] propose a method for estimating the DMSE of model-dependent small area predictors. The proposed method models the DMSE as a function of known area statistics by repeatedly drawing samples from appropriately generated synthetic populations, and then applies the model to the original sample. The procedure follows a method of bias correction developed by Pfeffermann and Correa^[18] for unit-level model dependent MSE estimation, with appropriate modifications.

The DMSE is defined as,

DMSE ({\hat{θ}}_{i}) = E_{D} [{({\hat{θ}}_{i} - θ_{i})}^{2} | θ_{i}]

(5.11)

where $E_{D}$ is the expectation under the randomization distribution over all possible sample selections from the finite population, with $θ_{i}$ held fixed. Simple calculations show that for the FH model with known model parameters $(σ_{u}^{2}, β)$ ,

λ_{i} (γ_{i}, β, σ_{D_{i}}^{2}) = D M S E ({\hat{θ}}_{i}) = γ_{i}^{2} σ_{D_{i}}^{2} + {(1 - γ_{i})}^{2} {(θ_{i} - {x^{'}}_{i} β)}^{2}

(5.12)

For known parameters, an unbiased estimator of $D M S E ({\hat{θ}}_{i})$ is,

{\hat{λ}}_{i}^{U B} = (2 γ_{i} - 1) σ_{D_{i}}^{2} + {(1 - γ_{i})}^{2} {(y_{i} - {x^{'}}_{i} β)}^{2}

(5.13)

Hence, for large $m$ , an approximately unbiased estimator of $D M S E ({\hat{θ}}_{i})$ is obtained by replacing the unknown parameters in (5.13) by their sample estimates, yielding the estimator,

{\hat{λ}}_{i}^{A U B} = (2 {\hat{γ}}_{i} - 1) σ_{D i}^{2} + {(1 - {\hat{γ}}_{i})}^{2} {(y_{i} - {x^{'}}_{i} {\hat{β}}_{G \hat{L} S})}^{2}

(5.14)

where ${\hat{β}}_{G \hat{L} S}$ is the GLS estimator but with $σ_{u}^{2}$ replaced by ${\hat{σ}}_{u}^{2}$ .

The DMSE estimator ${\hat{λ}}_{i}^{A U B}$ is expected to be unstable for small $m$ and large sampling variances $σ_{D_{i}}^{2}$ . Therefore, Pfeffermann and Ben-Hur^[19] propose an alternative estimator, which is constructed by applying the following steps:

Step 1. Estimate ( ${\hat{σ}}_{u}^{2}, \hat{β}$ ) based on the original sample. Generate a large number R of values $σ_{u r}^{2}, β_{r}$ from neighbourhoods around ${\hat{σ}}_{u}^{2}, \hat{β}$ , which are expected to include the true values underlying the hypothetical model generating the population values.

Step 2. Generate pseudo area means, $θ_{r i} = {x^{'}}_{i} β_{r} + u_{r i}$ ; $r = 1, ..., R;$ $i = 1, ..., m$ , using the same covariates as in the original sample.

Step 3. For each pseudo population of area means, generate $J$ parametric bootstrap samples, $y_{r i j} = θ_{r i} + e_{r i j} = {x^{'}}_{i} β_{r} + u_{r i} + e_{r i j}$ ; $j = 1, ..., J,$ $r = 1, ..., R,$ $i = 1, ..., m$ .

Step 4. For each bootstrap sample, estimate ${\hat{β}}_{r j}, {\hat{σ}}_{u r j}^{2}$ and compute the FH predictor, ${\hat{θ}}_{r i j} = {\hat{γ}}_{r i j} y_{r i j} + (1 - {\hat{γ}}_{r i j}) {x^{'}}_{i} {\hat{β}}_{r j}$ ; ${\hat{γ}}_{r i j} = {\hat{σ}}_{u r j}^{2} {({\hat{σ}}_{u r j}^{2} + σ_{D_{i}}^{2})}^{- 1}$ .

Step 5. Approximate the DMSE of the FH predictor ${\hat{θ}}_{r i} = {\hat{γ}}_{r i} y_{r i} + (1 - {\hat{γ}}_{r i}) {x^{'}}_{i} {\hat{β}}_{r}$ ; ${\hat{γ}}_{r i} = {\hat{σ}}_{u r}^{2} {({\hat{σ}}_{u r}^{2} + σ_{D_{i}}^{2})}^{- 1}$ by

D M S E_{r i} ({\hat{θ}}_{r i}) = \frac{1}{J} \sum_{j = 1}^{J} {({\hat{θ}}_{r i j} - θ_{r i})}^{2}

(5.15)

Step 6. Search for a function $q_{l} (\cdot) = D \hat{M} S E_{q l, r i} ({\hat{θ}}_{r i})$ of known predictors, which best predicts $D M S E_{r i} ({\hat{θ}}_{r i})$ (Eq. 5.15), among plausible functions $q_{l} (\cdot)$ . (See Section 6.2.)

Step 7. Apply the chosen function to the original sample to obtain an estimator of the DMSE of the FH predictor defined in (2.8).

Pfeffermann and Ben-Hur^[19] designed several simulation studies comparing their proposed method to other DMSE estimators proposed in the literature.

6. Simulation Studies

6.1 Simulation Setup and Results for MSE Estimation Under the Model

In order to assess and compare the performance of the model dependent MSE estimators considered in Sections 4 and 5, we conducted a Monte Carlo simulation study, designed as follows:

We generated a large number $Q$ of sets of true population means and corresponding direct estimators ${(θ_{i}^{(q)}, y_{i}^{(q)}, i = 1, ..., m, q = 1, ..., Q}$ , by use of the FH model (2.1). Following Datta et al.^[5] (hereafter DRS), we used for convenience a model with no auxiliary variables, such that $x_{i}^{'} β = μ = 0$ (but assumed unknown and hence estimated for each set of samples). We considered a total of $m = 15$ areas, divided into 5 groups of 3 areas in each group, with different sampling error variances $σ_{D_{i}}^{2}$ in different groups. The sampling error variances are $σ_{D}^{2} (g) = {2.0, 0.6, 0.5, 0.4, 0.2}$ , $g = 1, ..., 5$ , and the variance of the random effects is $σ_{u}^{2} = 1$ , same as under Pattern $b$ in DRS. Since the three areas in each group are exchangeable, the results reported in Tables 1 –6 below are averages over the three areas in each group.

We consider 3 combinations of distributions for the random effects, $u_{i}$ , and the sampling errors, $e_{i}$ : (i) both sets of error terms are generated from normal distributions; (ii) the sampling errors are generated from normal distributions but the random effects from a location exponential distribution; (iii) both sets of error terms are sampled from location exponential distributions. The location exponential distributions were set such that the variances are the same as the variances of the corresponding normal distributions; $E (- 1, 1)$ for the random effects and $E [- σ_{D} (g), σ_{D}^{2} (g)]$ for the sampling errors. The second and third combinations of distributions are considered in order to study the robustness of the various MSE estimators to deviations from the normality assumptions underlying the original FH model.

We started by generating $Q = 50, 000$ sets of true area means and direct estimators for each of the three cases (i)–(iii), and computed the EBLUP, using ${\hat{σ}}_{P R}^{2}$ and ${\hat{σ}}_{F H}^{2}$ for estimating $σ_{u}^{2}$ . This enabled us to approximate the true MSE of the corresponding two EBLUP predictors.

Next, for each of the three combinations of model error distributions and for each of the MSE estimators, we generated $\tilde{Q}$ new sets of direct estimators and used them to compute the MSE estimators and their relative bias and percent RMSE. We generated $\tilde{Q} = 2, 000$ sets for the parametric bootstrap MSE estimators (Eqs. 4.3, 4.5 and 5.10), and $\tilde{Q} = 10, 000$ for the other estimators.

Tables 1 –6 summarize the outcomes of the simulation study. The results refer to the following estimators: Prasad-Rao MSE estimator (Eq. 5.2), the DRS estimator (Eq. 5.3), the Jackknife estimator (Eq. 5.7), the Nonparametric Bootstrap estimator (Eq. 4.9) and the three bootstrap estimators–the parametric bootstrap estimator Eq. (4.3) (adding twice Eq. (4.4) when sampling from the exponential distributions), the Double bootstrap estimator Eq. (5.10) and the first stage bootstrap estimator Eq. (5.8). The simulation results of the alternative parametric bootstrap estimator Eq. (4.5) are not presented as they are very similar to the results of the parametric bootstrap estimator Eq. (4.3).

For the parametric bootstrap estimator and the first stage bootstrap estimator we used $B = 500$ replications. For the Double parametric bootstrap estimator, we used $B_{1} = 250$ and $B_{2} = 1$ replications. (See Remark 14). Note that the total number of replications of the double bootstrap estimator is 500, similarly to the other two bootstrap estimators. All the estimators defined above are presented separately when estimating $σ_{u}^{2}$ by ${\hat{σ}}_{P R}^{2}$ and when estimating $σ_{u}^{2}$ by ${\hat{σ}}_{F H}^{2}$ .

We used the correction ${\hat{σ}}_{P R}^{2} = \max ({\tilde{σ}}_{P R}^{2}, 0)$ for negative estimates in Eq. (3.1). We didn’t face negative variance estimators for ${\hat{σ}}_{F H}^{2}$ , because the iterative algorithm we used produced positive estimators for positive starting values. When the error terms are drawn from the exponential distributions, we show the results obtained for the three bootstrap-based estimators both when the error terms are wrongly assumed to be drawn from the normal distribution and when they are drawn from the correct distributions.

Tables 1 –3 show the true MSEs and the percent bias of the various MSE estimators, separately for each of the sampling variances $σ_{D}^{2} (g)$ and the 3 combinations of the model error distributions defined above. Tables 4 –6 show the corresponding percent root MSE (RMSE) of the MSE estimators.

Table 1 shows the results for the case where the distributions of the two error terms are normal. As can be seen, in this case the biases are low, except for the Taylor estimator with ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$ , where the bias increases quite drastically as the sampling variance decreases, similarly to the results in DRS. All the estimators, except the nonparametric bootstrap estimator NPB and the double bootstrap estimator DPB have generally a smaller bias for the case where ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$ than for the case where ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2},$ particularly for the smaller variances $σ_{D}^{2} (g)$ . The estimator DPB has generally the smallest bias, with NPB and JK-ACL coming next. The parametric bootstrap estimator PB performs likewise well when ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$ , but less so when ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$ . As expected, the first stage bootstrap estimator, $m s e^{B P 1}$ has a relatively large negative bias in all the cases. This estimator has bias of order $O (1 / m)$ .

Table 1.

True MSE of EBLUP and Percent Relative Bias of MSE estimators based on Taylor approximations Eqs. (5.2, 5.3) (Taylor), Jackknife Eq. (5.7) (JK-ACL), Nonparametric Bootstrap Eq. (4.9) (NPB), Parametric Bootstrap Eq. (4.3) (PB), Double Parametric Bootstrap Eq. (5.10) (DPB) and $m s e^{P B 1}$ Eq. (5.8). 15 Areas, model errors generated from normal distributions.

	${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$					${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$
$σ_{D}^{2} (g)$	2.0	0.6	0.5	0.4	0.2	2.0	0.6	0.5	0.4	0.2
100 × MSE	78.3	43.7	38.8	33.4	19.9	76.9	42.1	37.0	31.6	18.1
Taylor	–0.7	6.8	9.3	12.5	34.5	–2.4	–0.8	0.3	0.4	2.7
JK-ACL	–3.2	–2.0	–1.2	–0.4	4.7	–1.6	–1.6	–0.7	–0.6	0.3
NPB	–0.7	–1.6	–1.3	–1.8	–0.3	1.4	0.6	1.7	1.4	2.5
PB	–2.9	–3.7	–3.1	–2.9	–1.3	–0.9	–1.1	–0.2	–0.3	0.7
DPB	–0.6	–1.2	–0.2	–0.7	0.6	–0.9	–0.1	–0.4	–0.7	0.5
$m s e^{P B 1}$	–8.5	–11.0	–10.2	–10.1	–7.9	–5.7	–7.0	–6.2	–6.3	–4.6

Table 2 shows the results obtained when the sampling errors are sampled from normal distributions, but the random effects are sampled from the location exponential distribution. The results of the parametric bootstrap estimators under the correct model, that is, when generating the bootstrap samples by drawing from the correct distribution of the random effects are labelled by E. The results when the random effects are wrongly assumed to be generated from the normal distribution are labelled by N.

Table 2.

True MSE of EBLUP and Percent Relative Bias of MSE estimators based on Taylor approximations Eqs. (5.2, 5.3) (Taylor), Jackknife Eq. (5.7) (JK-ACL), Nonparametric Bootstrap Eq. (4.9) (NPB), Parametric Bootstrap Eq. (4.3) (PB), Double Parametric Bootstrap Eq. (5.10) (DPB) and $m s e^{P B 1}$ Eq. (5.8). 15 Areas, Random Effects generated from location exponential distribution, Sampling Errors generated from normal distribution. Parametric Bootstrap samples generated from the true random effect distribution (E) and by assuming a normal distribution (N).

		${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$					${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$
$σ_{D}^{2} (g)$		2.0	0.6	0.5	0.4	0.2	2.0	0.6	0.5	0.4	0.2
100 × MSE		74.2	40.7	36.9	31.6	19.4	74.7	39.4	35.1	29.9	17.3
Taylor		0.4	14.6	17.2	24.9	69.4	–6.0	0.0	0.5	2.0	8.6
JK-ACL		–0.1	1.7	0.8	2.3	6.7	–0.7	2.2	2.2	2.6	3.6
NPB		3.1	2.5	1.2	2.2	2.1	1.5	4.3	4.2	4.6	5.5
PB		–4.2	–4.3	–5.4	–4.6	–4.9	–5.4	–1.9	–1.9	–1.3	0.1
DPB	N	–6.0	–5.4	–6.4	–5.8	–5.3	–6.0	–1.8	–2.8	–1.8	–0.5
$m s e^{P B 1}$		–8.5	–10.0	–11.1	–10.4	–10.2	–9.0	–7.0	–7.1	–6.5	–5.0
PB*		–3.6	–4.7	–5.8	–4.8	–4.8	–6.1	–3.3	–3.3	–2.6	–0.7
DPB	E	0.1	0.3	–2.2	–0.9	–1.9	–5.5	–3.0	–4.4	–3.4	–1.6
$m s e^{P B 1}$		–12.9	–16.0	–17.2	–15.9	–14.6	–14.1	–13.0	–12.0	–13.0	–9.3

* Added twice Eq. (4.4) to Eq. (4.3)

The results in Table 2 reveal that all the estimators, except for the Taylor based estimator with ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$ perform relatively well in this case as well, despite the non-normality of the random effects, although the biases are generally higher than in Table 1, where the random effects are generated from the normal distribution. The large biases observed for the Taylor based estimator with ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$ are somewhat surprising in view of the theoretical results of Lahiri and Rao^[13] (see Remark 10), but notice that in the present experiment we only consider 15 areas. When ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$ , the estimator DPB E performs best, but the estimator DPB N has relative large biases of 5–6 percent in this case. The estimator JK-ACL generally performs well, but it has a relatively large bias of 6.7 percent when $σ_{D}^{2} (g) = 0.2$ (smallest sampling error variance). When ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$ , JK-ACL performs overall the best. The estimator DPB N performs well except for the case of $σ_{D}^{2} (g) = 2$ where the bias is –6.0 percent but surprisingly, the estimator DPB E which draws the random effects from the correct distribution performs less satisfactorily. The NPB estimator performs somewhat better with ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$ than with ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$ . The PB estimator performs well with ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$ (somewhat better than DPB), except when $σ_{D}^{2} (g) = 2$ , in which case it has a bias of –5.4 and –6.1 percent. As in Table 1, the estimator $m s e^{P B 1}$ has large negative biases, much larger in absolute values than all the other estimators.

Table 3 shows the results obtained for the case where both the random effects and the sampling errors are generated from the location exponential distributions. The relative biases in this table are much larger than in Tables 1 and 2, except in the case of PB E and DPB E with ${\hat{σ}}_{F H}^{2}$ , which use the correct distributions for generating the bootstrap samples. For ${\hat{σ}}_{F H}^{2}$ , the estimators JK-ACL, NPB and to a lesser extent also Taylor, also perform relatively well, except in the case $σ_{D}^{2} (g) = 2$ . Interestingly, for ${\hat{σ}}_{F H}^{2}$ , the PB estimator performs better than DPB E and the estimator $m s e^{P B 1}$ E has smaller percent biases than some of the other estimators, about –11.5 percent (–8.4 percent when $σ_{D}^{2} (g) = 0.2$ ). Thus, at least for a small number of areas as in the present experiment, all the other methods are sensitive to the deviation from normality of the sampling error distribution.

Table 3.

True MSE of EBLUP and Percent Relative Bias of MSE estimators based on Taylor approximations Eqs. (5.2, 5.3) (Taylor), Jackknife Eq. (5.7) (JK-ACL), Nonparametric Bootstrap Eq. (4.9) (NPB), Parametric Bootstrap Eq. (4.3) (PB), Double Parametric Bootstrap Eq. (5.10) (DPB) and $m s e^{P B 1}$ Eq. (5.8). 15 Areas, Random Effects and Sampling Errors generated from location exponential distributions. Parametric Bootstrap samples generated from the true error distributions (E) and by assuming a normal distribution (N).

		${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$					${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$
$σ_{D}^{2} (g)$		2.0	0.6	0.5	0.4	0.2	2.0	0.6	0.5	0.4	0.2
100 × MSE		92.7	45.4	40.3	34.7	20.8	89.9	43.1	38.1	32.2	18.1
Taylor		–22.1	2.8	8.9	18.0	80.5	–23.1	–10.4	–9.1	–6.7	5.0
JK-ACL		–21.3	–11.9	–10.7	–9.6	–0.3	–17.5	–8.1	–7.7	–6.4	–1.9
NPB		–20.4	–12.8	–12.0	–11.5	–9.1	–15.8	–5.6	–5.1	–4.0	0.5
PB		–26.6	–19.1	–18.4	–18.0	–15.9	–23.7	–13.6	13.0	–11.7	–6.9
DPB	N	–28.8	–20.8	–19.5	–18.4	–16.6	–25.3	–14.3	13.3	–12.2	–7.2
$m s e^{P B 1}$		–29.1	–22.9	–22.3	–22.1	–19.7	–26.6	–17.7	–17.0	–16.0	–11.2
PB*		–7.9	–8.1	–8.1	–8.3	–6.8	–4.2	–1.6	–1.9	–1.5	0.8
DPB	E	–8.1	–9.1	–8.0	–7.3	–6.3	–8.1	–3.5	–3.7	–4.8	–1.1
$m s e^{P B 1}$		–14.6	–17.8	–17.7	–17.9	–16.0	–11.6	–11.2	11.4	–11.1	–8.4

* Added twice Eq. (4.4) to Eq. (4.3)

Tables 4 –6 show the percent Root MSE (RMSE) of the MSE estimators, under the three combinations of the distributions of the random effects and the sampling errors.

For the case where both distributions are normal (Table 4), the estimator DPB has the lowest RMSE, but for $σ_{D}^{2} (g) = 0.6, 0.5, 0.4$ , the RMSEs of all the estimators, including DPB and $m s e^{P B 1}$ are of similar magnitude of between 20 to 25 percent. When $σ_{D}^{2} (g) = 2$ , the RMSE of the DPB estimator is about 40.3 percent and the percent RMSE of all the other estimators are in the range of 50 to 57 percent. When $σ_{D}^{2} (g) = 0.2$ , all the percent RMSEs including the DPB are in the range of 30–35 percent, (30 to 31 percent when ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$ ). We also note that unlike in the case of the relative biases, the RMSEs of all the estimators are generally similar when estimating ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$ or ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$ .

Table 4.

True MSE of EBLUP and Percent Root MSE of MSE estimators based on Taylor approximations Eqs. (5.2, 5.3), (Taylor), Jackknife Eq. (5.7) (JK-ACL), Nonparametric Bootstrap Eq. (4.9) (NPB), Parametric Bootstrap Eq. (4.3) (PB), Double Parametric Bootstrap Eq. (5.10) (DPB) and $m s e^{P B 1}$ Eq. (5.8). 15 Areas, errors generated from normal distributions.

	${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$					${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$
$σ_{D}^{2} (g)$	2.0	0.6	0.5	0.4	0.2	2.0	0.6	0.5	0.4	0.2
100 × MSE	78.3	43.7	38.8	33.4	19.9	76.9	42.1	37.0	31.6	18.1
Taylor	50.1	20.1	19.7	21.3	34.5	48.4	21.3	21.0	22.1	29.9
JK-ACL	53.0	23.3	22.7	23.8	32.1	50.5	22.2	21.7	22.8	30.4
NPB	54.9	24.4	23.5	24.0	30.5	57.3	25.9	24.8	24.8	30.7
PB	51.6	23.5	23.0	23.8	30.6	49.6	22.0	21.6	22.7	30.3
DPB	40.2	19.5	19.9	21.7	30.0	40.4	19.7	20.1	21.9	30.2
$m s e^{P B 1}$	46.6	23.0	23.2	24.5	31.5	46.0	21.9	22.1	23.6	31.0

Similar relative performance patterns are found in the case where the sampling errors have a normal distribution but the random effects are generated from the location exponential distribution (Table 5), and in the case where the two errors are generated from the location exponential distribution (Table 6). Note first that the percent RMSEs in Table 5 are of similar magnitude to the percent RMSEs in Table 4, but they are larger in Table 6. Thus, generating the sampling errors from the exponential distribution increases the RMSE of the MSE estimators. The estimator DPB has again the lowest RMSEs and $m s e^{P B 1}$ also performs relatively well in both the tables. For $σ_{D}^{2} (g) = 2$ , the NPB estimator has larger percent RMSEs than the other estimators in all the three tables. The RMSEs of DPB N are similar to the RMSEs of DPB E in Table 5, but much smaller in Table 6 when $σ_{D}^{2} (g) = 2.$

Table 5.

True MSE of EBLUP and Percent Root MSE of MSE estimators based on Taylor approximations Eqs. (5.2, 5.3) (Taylor), Jackknife Eq. (5.7) (JK-ACL), Nonparametric Bootstrap Eq. (4.9) (NPB), Parametric Bootstrap Eq. (4.3) (PB), Double Parametric Bootstrap Eq. (5.10) (DPB) and $m s e^{P B 1}$ Eq. (5.8). 15 Areas, Random Effects generated from location exponential distribution, Sampling Errors generated from normal distribution. Parametric Bootstrap samples taken from the true random effect distribution (E) and by assuming a normal distribution (N).

		${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$					${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$
$σ_{D}^{2} (g)$		2.0	0.6	0.5	0.4	0.2	2.0	0.6	0.5	0.4	0.2
100 × MSE		74.2	40.7	36.9	31.6	19.4	74.7	39.4	35.1	29.9	17.3
Taylor		52.4	19.6	18.9	20.6	37.9	51.2	22.2	21.4	21.9	28.1
JK-ACL		59.9	25.1	23.7	23.7	30.6	58.7	24.2	23.0	23.0	29.0
NPB		61.1	26.4	24.7	24.1	28.6	65.3	28.9	26.9	25.9	29.7
PB		53.8	24.7	23.8	24.0	29.3	52.0	23.3	22.6	23.2	29.4
DPB	N	34.7	18.5	19.3	21.3	28.8	37.0	19.4	20.0	21.7	29.2
$m s e^{P B}$		48.6	23.3	23.1	23.9	29.9	48.8	22.9	22.7	23.7	30.0
PB*		55.0	25.0	24.0	24.1	29.3	52.1	23.7	23.0	23.5	29.6
DPB	E	38.6	18.4	18.9	20.5	28.3	37.0	19.3	20.1	21.9	29.4
$m s e^{P B}$		46.7	23.4	23.4	24.4	30.3	45.2	23.1	23.2	24.5	30.6

* Added twice Eq. (4.4) to Eq. (4.3)

Table 6.

True MSE of EBLUP and Percent Root MSE of MSE estimators based on Taylor approximations Eqs. (5.2, 5.3) (Taylor), Jackknife Eq. (5.7) (JK-ACL), Nonparametric Bootstrap Eq. (4.9) (NPB), Parametric Bootstrap Eq. (4.5) (PB), Double Parametric Bootstrap Eq. (5.10) (DPB) and $m s e^{P B 1}$ Eq. (5.8). 15 Areas, Random Effects and Sampling Errors generated from location exponential distributions. Parametric Bootstrap samples taken from the true error distributions (E) and by assuming a normal distribution (N).

		${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$					${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$
$σ_{D}^{2} (g)$		2.0	0.6	0.5	0.4	0.2	2.0	0.6	0.5	0.4	0.2
100 × MSE		92.7	45.4	40.3	34.7	20.8	89.9	43.1	38.1	32.2	18.1
Taylor		52.6	24.7	25.0	27.3	45.4	52.0	27.7	27.8	28.9	35.1
JK-ACL		65.3	31.9	31.4	32.2	39.7	63.3	29.9	29.3	30.1	36.3
NPB		63.0	33.2	32.5	32.9	37.9	68.9	34.1	32.8	32.6	37.0
PB		55.4	32.2	32.3	33.2	38.8	53.0	29.4	29.5	30.7	36.9
DPB	N	31.6	26.8	28.4	30.7	38.3	33.5	25.4	26.8	29.1	36.7
$m s e^{P B}$		50.4	31.2	31.7	33.1	39.1	49.6	29.2	29.7	31.1	37.5
PB*		63.6	30.9	30.6	31.5	37.4	63.8	28.2	27.9	28.9	35.9
DPB	E	45.9	25.1	26.4	28.5	36.6	45.6	24.5	25.5	27.9	35.9
$m s e^{P B}$		54.6	29.8	30.6	32.2	38.6	56.0	27.9	28.4	30.0	37.0

*Added twice Eq. (4.4) to Eq. (4.3).

6.2 Simulation Setup and Results for Design-Based MSE Estimation

In Section 5.4, we considered the estimation of the design-based MSE of the FH model-dependent EBLUP. We conducted a second simulation study in order to assess the performance of the following three estimators: the approximately unbiased estimator ${\hat{λ}}_{i}^{A U B}$ (Eq. 5.14) based on L=10,000 simulations; the average of the estimators (5.15) over L=500 simulations, denoted $A v (D M S E)$ and the estimator proposed by Pfeffermann and Ben-Hur^[19], denoted $D M S E_{P - B}$ . We used L=500, R=100, J=250 for the computation of the third estimator. See Section 5.4 for the definitions of R and J.

We chose the function $q_{l} (\cdot) = D M S E_{q_{l}, r i} ({\hat{θ}}_{r i})$ (Step 6) among linear regression functions by combination of stepwise regression and cross validation techniques, with the following plausible predictors, ${\hat{a}}_{r}, {\hat{σ}}_{u}^{2}, {\hat{a}}_{r}, {\hat{σ}}_{D_{i}}^{2}, {\hat{γ}}_{r i}, {\hat{γ}}_{r i}^{2}, (1 - {{\hat{γ}}_{r i})}^{2}, {\hat{θ}}_{r i}, {({\hat{θ}}_{r i} - μ_{r})}^{2}, {(y_{r i} - μ_{r})}^{2}$ (is the estimator of the constant term in simulation r), and dependent variables, $D_{i}, \log (D_{i}), \arcsin (\sqrt{D_{i} / 100}), \frac{1}{D_{i}}, \sqrt{D_{i}}, \frac{1}{\sqrt{D_{i}}}$ ; denoting $D M S E_{r i} ({\hat{θ}}_{r i})$ by $D_{i}$ .

The model, number of areas and the distributions of the random effects and the sampling errors are the same as in Section 6.1, estimating $σ_{u}^{2}$ by ${\hat{σ}}_{P R}^{2}$ and ${\hat{σ}}_{F H}^{2}$ . The true design-based MSEs have been computed based on L=50,000 simulated values of true area means and corresponding sample estimators. The selection of the function $q_{l} (\cdot)$ has been applied for each distribution of the random effects and the sampling errors.

Tables 7 –9 present the results obtained for the three distributions of the random effects and the sampling errors. we used the neighbourhoods $μ \in (- 2, 2), σ_{u}^{2} \in (0.1, 2)$ for the computation of the means $θ_{r i} = μ + u_{r i}$ and the direct estimators $y_{r i j} = θ_{r i} + e_{r i j} = μ + u_{r i} + e_{r i j}$ . (Steps 1–3, see Section 5.4).

REMARK 15. The results in Tables 7 ߝ 9 are averages over M=10 finite populations. In our simulation study, we only consider 15 areas, and the true area means are $θ_{i} = u_{i},$ such that it suffices that a small number of the random effects take extreme values to destabilise the results.

The first noteworthy outcome emerging from Tables 7 –9 is that the True DMSEs are systematically somewhat lower than the corresponding true MSEs in Tables 1 –3 under the model, except for the case where ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$ $σ_{D}^{2} (g) = 0.2$ . This outcome can be explained by the fact that the MSEs under the model account also for the distribution of the random effects, which are held fixed under the design-based approach. All the three estimators in Tables 7 –9 have a negative bias, with the absolute percent biases being less than 6 percent in Table 7, less than 10 percent in Table 8, and less than 14 percent in Table 9. The three estimators perform quite similarly in Tables 7 and 8, and generally better with than with ${\hat{σ}}_{F H}^{2}$ than with ${\hat{σ}}_{P R}^{2}$ , but the proposed estimator $D M S E_{P - B}$ dominates the other two estimators in Table 9, when both the random effects and the sampling errors are generated from location exponential distributions. Notice the relative good performance of the approximately unbiased estimator $λ_{i}^{A U B}$ in Tables 7 and 8.

Table 7.

True DMSE of EBLUP and Percent Relative Bias of the following DMSE estimators: the approximately unbiased estimator ${\hat{λ}}_{i}^{A U B}$ (Eq. 5.14), the average estimator $A v (D M S E)$ (average of estimators 5.15) and the proposed estimator $D M S E_{P - B}$ . 15 areas, random effects and sampling errors generated from normal distributions.

	${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$					${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$
$σ_{D}^{2} (g)$	2.0	0.6	0.5	0.4	0.2	2.0	0.6	0.5	0.4	0.2
100 × MSE	70.7	38.0	36.9	31.6	19.3	69.3	37.3	35.4	30.3	18.0
${\hat{λ}}_{i}^{A U B}$	–4.6	–3.9	–3.4	–2.8	–2.2	–4.1	–3.2	–2.9	–2.7	–2.4
$A v D M S E$	–6.1	–5.0	–3.9	–3.3	–2.8	–5.2	–3.9	–2.5	–1.9	–1.7
$D M S E_{P - B}$	–5.9	–4.0	–3.7	–3.6	–1.6	–5.3	–2.7	–2.8	–3.3	–1.2

Table 8.

True DMSE of EBLUP and Percent Relative Bias of the following DMSE estimators: the approximately unbiased estimator ${\hat{λ}}_{i}^{A U B}$ (Eq. 5.14), the average estimator $A v (D M S E)$ (average of estimators 5.15) and the proposed estimator $D M S E_{P - B}$ . 15 areas, random effects generated from location exponential distribution, sampling errors generated from normal distributions.

	${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$					${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$
$σ_{D}^{2} (g)$	2.0	0.6	0.5	0.4	0.2	2.0	0.6	0.5	0.4	0.2
100 × MSE	71.3	38.4	33.7	28.4	19.0	71.4	37.1	32.7	27.6	17.4
${\hat{λ}}_{i}^{A U B}$	–9.2	–8.4	–7.3	–6.5	–6.2	–8.8	–6.9	–5.0	–5.7	–4.8
$A v D M S E_{r i}$	–10.3	–8.9	–8.1	–6.9	–6.1	–10.1	–7.7	–7.3	–6.1	–2.6
$D M S E_{P - B}$	–9.8	–8.1	–7.5	–6.3	–5.2	–9.3	–7.5	–7.1	–4.9	–3.1

Table 9.

True DMSE of EBLUP and Percent Relative Bias of the following DMSE estimators: the approximately unbiased estimator $λ_{i}^{A U B}$ (Eq. 5.14), the average estimator $A v (D M S E)$ (average of estimators 5.15) and the proposed estimator $D M S E_{P - B}$ . 15 areas, random effects and sampling errors generated location exponential distributions.

	${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$					${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$
$σ_{D}^{2} (g)$	2.0	0.6	0.5	0.4	0.2	2.0	0.6	0.5	0.4	0.2
100 × MSE	90.8	41.9	35.6	29.7	20.2	87.5	40.6	34.6	29.1	18.1
${\hat{λ}}_{i}^{A U B}$	–13.7	–12.7	–12.5	–12.1	–11.4	–12.9	–11.8	–11.4	–10.2	–10.4
$A v D M S E$	–13.4	–12.6	–12.2	–12.5	–11.6	–7.1	–9.1	–7.5	–8.6	–7.1
$D M S E_{P - B}$	–11.9	–10.4	–9.5	–11.2	–8.8	–6.1	–7.5	–5.8	–6.4	–5.7

We conclude from this simulation study that it is possible to estimate the design-based MSE of model dependent estimators with good quality, even with a small number of areas. See Pfeffermann and Ben-Hur^[19] for the performance of their proposed estimator in the case of the unit-level generalised linear mixed model.

7. Summary Remarks

In this article, we compare a large number of methods for estimating the MSE of the EBLUP under the Fay-Herriot model. The first important result of this study is that the EBLUP that uses the estimator ${\hat{σ}}_{F H}^{2}$ for estimating the variance of the random effects has somewhat lower true MSEs than the EBLUP that uses the estimator ${\hat{σ}}_{P R}^{2}$ , although not by much (compare the true MSEs in the various tables). On the other hand, no single method of MSE estimation dominates all the other methods in terms of bias and RMSE (of the MSE estimators).

In what follows we discuss first the MSE estimators under the model. When generating the error terms from the correct distributions, the double bootstrap estimator DPB has generally the lowest bias and percent RMSE. The Jackknife estimator JK-ACL and the nonparametric estimator NPB have generally small biases when the sampling errors have a normal distribution, but the biases increase in the case where they are generated from the location exponential distribution, as is the case with all the other estimators. The parametric bootstrap estimator PB has larger biases than the previous two estimators when ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$ , but similar biases when ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$ and the sampling errors are generated from normal distributions, except in Table 2 when $σ_{D}^{2} (g) = 2$ . Much smaller bias when $σ_{D}^{2} (g) = 0.2$ . When the two error terms are generated from the location exponential distribution and ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$ , the PB estimator that accounts for these distributions has somewhat lower biases than the corresponding DPB estimator. The Taylor based estimators perform relatively well in terms of bias when ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$ , but less so in the case where ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{P R}^{2}$ . The Taylor based estimator has a large negative bias in Table 3 when $σ_{D}^{2} (g) = 2$ , as have all the estimators that assume normality of the error terms. However, the percent RMSEs of the MSE estimators are in most cases quite similar with the two variance estimators.

We emphasize again that our results are restricted to 15 small areas. Most of the published studies on the estimation of the MSE of the EBLUP in SAE, focus on the bias of the MSE estimators. Clearly, the bias is the dominant contributor to the MSE when the number of areas is large, but not when it is small, as in the present study. We recognize that analytical comparisons of the MSE of MSE estimators to the right order are complicated, but this fundamental quality of MSE estimators should be explored empirically. As our results indicate, a MSE estimator with negligible bias may actually have a larger variance and hence a larger MSE than another estimator with a large bias.

The present article explores the effect of deviations from normality of the distributions of the model error terms on the performance of the MSE estimators. All the methods, except for PB and DPB when based on the correct distribution of the sampling errors, and to a lesser extent also the JK-ACL, NPB and the Taylor method with ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$ and $σ_{D}^{2} (g) \neq 2$ , yield estimators with large bias when the distribution of the sampling errors is not normal. Clearly, the use of bootstrap estimators based on the correct distribution of the sampling errors requires identifying this distribution, which is not always simple, especially with small number of areas. The JK-ACL, NPB and the Taylor estimator, with ${\hat{σ}}_{u}^{2} = {\hat{σ}}_{F H}^{2}$ , seem to be more robust to deviations from normality of the sampling errors.

In the present article we also study the performance of three plausible estimators of the design-based MSE. All the estimators perform well when the sampling errors are generated from normal distributions, but the proposed estimator $D M S E_{P - B}$ performs better than the other two estimators when the sampling errors are generated from the location exponential distribution. Studying the performance of the three estimators and possibly other estimators proposed in the literature for estimating the design-based MSE, including for areas with no samples, when the distributions of the error terms are different from normal but normality is assumed, need to be investigated.

Footnotes

Declaration of Conflicting Interests

The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The authors received no financial support for the research, authorship and/or publication of this article.

Notes

References

Butar

, Lahiri

On measures of uncertainty of empirical Bayes small area estimators. J Stat Plann Infer 2003; 112: 63–76.

Chen

, Lahiri

A weighted Jackknife MSPE in small area estimation. Proc Surv Res Meth Sec ASA 2002; 433–477. Url http://www.asasrms.org/Proceedings/y2002/Files/JSM2002-001127.pdf

Chen

, Lahiri

A Comparison of different MSPE estimators of EBLUP for the Fay-Herriot model. Proc Surv Res Meth Sec ASA 2003; 905–911.

Datta

, Lahiri

A unified measure of uncertainty of estimated best linear unbiased predictors in small area estimation problems. Stat Sinica 2000; 10: 613–627.

Datta

, Rao

JNK

, Smith

DD.

On measuring the variability of small area estimators a basic area level model. Biometrika 2005; 95: 183–196.

Erciulescu AL

Fuller AF.

Small area prediction under alternative model specifications. Stat Transit 2014; 17: 9–24.

Fay

, Herriot

RA.

Estimation of income from small places: an application of James-Stein procedures to census data. J Amer Stat Assoc 1979; 74: 269–277.

Harville

DA.

Decomposition of prediction error. J Amer Stat Assoc 1985; 80: 132–138.

Hall

, Maiti

On parametric bootstrap methods for small area prediction. J Royal Stat Soc, Series B 2016; 68: 221–238.

10.

Jiang

REML estimation: asymptotic behavior and related topics. Annal Stat 1996; 24: 255–286.

11.

Jiang

, Lahiri

, Wan

A unified jackknife method. Annal Stat 2002; 30: 1782–1810.

12.

Jiang

, Lahiri

Mixed model prediction and small area estimation (with discussion). Test 2006; 15: 1–96.

13.

Lahiri

, Rao

JNK.

Robust estimation of mean square error of small area estimators. J Amer Stat Assoc 1995; 90: 758–766.

14.

, Lahiri

An adjusted maximum likelihood method for solving small area estimation problems. J Multivar Anal 2010; 101: 882–892.

15.

Lohr

, Rao

JNK.

Jackknife estimation of mean squared error of small area predictors in nonlinear mixed models. Biometrika 2009; 96: 457–468.

16.

Pfeffermann

, Nathan

Regression analysis of data from a cluster Sample. J Amer Stat Assoc 1981; 76: 681–689.

17.

Pfeffermann

, Tiller

RB.

Bootstrap approximation to prediction MSE for state-space models with estimated parameters. J Time Series Anal 2005; 26: 893–916.

18.

Pfeffermann

, Correa

Empirical bootstrap bias correction and estimation of prediction mean square error in small area estimation. Biometrika 2012; 99: 1–16.

19.

Pfeffermann

, Ben-Hur

Estimation of randomisation mean square error in small area estimation. Inter Stat Rev 2018; 87: 31–49.

20.

Prasad

NGN

, Rao

JNK.

On estimation of the mean squared error of small-are estimators. J Amer Stat Assoc 1990; 85: 163–171.