arXiv:1612.04025v1 [stat.ME] 13 Dec 2016

Second-order unbiased naive estimator of mean
squared error for EBLUP in small-area estimation
Masayo Yoshimori Hirose
The Institute of Statistical Mathematics

Abstract
An empirical best linear unbiased prediction (EBLUP) estimator is
utilized for efficient inference in small-area estimation. To measure its
uncertainty, we need to estimate its mean squared error (MSE) since
the true MSE cannot generally be derived in a closed form. The naive
MSE estimator, one of the estimators available for small-area inference, is unlikely to be chosen, since it does not achieve the desired
asymptotic property, namely second-order unbiasedness, although it
maintains strict positivity and tractability. Therefore, users tend to
choose the second-order unbiased MSE estimator. In this paper, we
seek a new adjusted maximum-likelihood method to obtain a naive
MSE estimator that achieves the required asymptotic property. To
obtain the result, we also reveal the relationship between the general
adjusted maximum-likelihood method for the model variance parameter and the general functional form of the second-order unbiased, and
strictly positive, MSE estimator. We also compare the performance of
the new method with that of the existing naive estimator through a
Monte Carlo simulation study. The results show that the new method
remedies the underestimation associated with the existing naive estimator.

Keywords, Adjusted maximum-likelihood method; Empirical best linear unbiased prediction; Fay–Herriot model; Linear mixed model; Mean
squared error.

1

Introduction

In recent decades, there has been high demand for reliable statistics on
smaller geographic areas and sub-populations where large samples are not
1

available. Considering the limited number of observations, a design-based
direct estimator is not reliable for such “small areas”—as they are called.
An empirical best linear unbiased prediction (EBLUP) estimator is widely
used as an efficient estimator based on a specific linear mixed model. It
would be quite interesting to use the mean squared error (MSE) of EBLUP
as a measure of its uncertainty. For small-area inference, its MSE needs
to be estimated with high accuracy since it is not generally derived in a
closed form. Given a consistent estimator of an unknown model variance
parameter, the MSE of EBLUP is always larger than that of the best linear unbiased prediction (BLUP) estimator which assumes a known model
variance parameter, under certain conditions (Kackar and Harville, 1984).
In most small-area applications, sufficient accuracy cannot be achieved by
ignoring this difference, which is of the order of O(m−1 ) for large m (number of areas). Moreover, the naive MSE estimator, a consistent estimator
substituted for the model variance parameter in the MSE of BLUP, lacks
second-order unbiasedness for sufficient asymptotic accuracy in small-area
estimation with large m.
Therefore, several second-order unbiased MSE estimators, with some
bias correction, are suggested in place of the naive estimator (Prasad and
Rao, 1990; Datta and Lahiri, 2000; Butar and Lahiri, 2003; Das et al.,
2004; Datta et al, 2005; Hall and Maiti, 2006; Li and Lahiri, 2010; Yoshimori and Lahiri, 2014). In particular, Yoshimori and Lahiri (2014) achieved
strictly positive variance estimation while maintaining a functional form of
the second-order unbiased MSE estimator proposed in Datta and Lahiri
(2000), by using the adjusted maximum-likelihood method.
Incidentally, a relevant question that arises is, Can the naive MSE estimator provide second-order unbiasedness through the adjusted maximumlikelihood method? To answer this question, this paper proposes a new
method for naive MSE estimation as (6)-(7) in Section 4, which achieves
the desired asymptotic property while maintaining strict positivity. To obtain the result, Section 3 provides a theorem to choose a suitable adjusted
maximum-likelihood method for a specified functional form of the secondorder unbiased and strictly positive MSE estimator, and vice versa (to choose
a suitable functional form of the second-order unbiased and strictly positive
MSE estimator for a specified adjusted maximum-likelihood method). Section 5 presents a performance comparison among certain MSE estimators,
including ours. The regularity conditions and all technical proofs are deferred to the appendix.

2

2

The uncertainty of EBLUP under the Fay–Herriot
model

The Fay–Herriot model (Fay and Herriot, 1979) is widely used for small-area
inference. For i = 1, . . . , m,
ind

Level 1 :yi |θi ∼ N (θi , Di );
ind

Level 2 :θi ∼ N (x′i β, A).

(1)

The level-1 model takes into account the sampling distribution of the direct
estimator yi for the ith small area. The true small-area mean for the ith
area, denoted by θi , is linked to providing area-specific auxiliary variables
xi = (xi1 , · · · , xip )′ in the level-2 model. In practice, the coefficient vector
β in Rp and the model variance parameter A in this linking model are
unknown. The assumption of a known Di often follows from the asymptotic
variances of the transformed direct estimates (Efron and Morris, 1975) or
from empirical variance modelling (Fay and Herriot, 1979). This model can
be rewritten as a specific linear mixed model:
yi = θi + ei = x′i β + ui + ei , i = 1, . . . , m,
where ui and ei are mutually independent with the normality assumption
iid
ind
ui ∼ N (0, A) and ei ∼ N (0, Di ). It is well known that among all linear unbiased predictors θ̂i of θi , BLUP yields the minimum MSE, which is defined
as E[(θ̂i − θi )2 ], where the expectation is defined with respect to the joint
distribution of y and θ under the Fay–Herriot model (1). We give the form
of BLUP as follows:
θ̂iB = (1 − Bi )yi + Bi x′i β̃,
Di
is called the shrinkage factor toward x′i β̃ from the direct
where Bi = A+D
i
estimator yi with β̃ = β̃(A) = (X ′ V −1 X)−1 X ′ V −1 y, y = (y1 . . . , ym )′ , X =
(x1 , . . . , xm )′ and V = diag{A + D1 , · · · , A + Dm }.
Since A is unknown in practice, the following EBLUP of θi is widely
used for small-area inference, with A replaced with its consistent estimator,
Â, in θ̂iB :
θ̂iEB = (1 − B̂i )yi + B̂i x′i β̂,

where B̂i = Di and β̂ = β̃(Â). Hereafter, the consistent estimator Â
Â+Di
also denotes an even-translation-invariant estimator for all β and y that
achieve an unbiasedness in the EBLUP, as in Kackar and Harville (1981).
To estimate the model variance parameter A, we can use the method of
3

moments estimator (Fay and Herriot, 1979; Prasad and Rao, 1990) and
the standard maximum-likelihood estimators, such as the profile maximumlikelihood (PML) and the residual maximum-likelihood (REML) estimators.
In particular, the REML estimator of A is preferred in terms of its higherorder asymptotic accuracy for large m. Let ÂRE denote the REML estimator
of A, obtained as
ÂRE = arg max LRE (A|y),
0≤A<∞

where the residual likelihood function is
LRE (A|y) = |X ′ V −1 X|−1/2 |V |−1/2 exp{−y ′ P y/2}
and P = V −1 − V −1 X(X ′ V −1 X)−1 X ′ V −1 .
However, the REML estimator of A has serious problem such that it
could be zero when m (number of small areas) is not large enough, even
though A = 0 is not realistic in the context of small-area estimation. In
order to avoid zero estimates, Li and Lahiri (2010) and Yoshimori and Lahiri
(2014) suggested using the specific adjusted maximum-likelihood estimators.
Let ÂLL and ÂY L denote the respective estimators, given as
Âad = arg max L̃(A)LRE (A),
0≤A<∞

where Âad ∈ {ÂLL , ÂY L }, and L̃(A) are adopted from their specific adjustment factors, L̃(A) = A for ÂLL and L̃(A) = arctan[tr(I − B∗ )]1/m for ÂY L
with B∗ = diag(B1 , . . . , Bm ).
The MSE of BLUP under the Fay–Herriot model can be derived in a
closed form as
h
i
M SEi [θ̂iB ] ≡ E (θ̂iB − θi )2 = g1i (A) + g2i (A),
D2

ADi
and g2i (A) = (A+Di i )2 x′i (X ′ V −1 X)−1 xi . Unlike BLUP,
where g1i (A) = A+D
i
EBLUP cannot generally provide a closed-form MSE, so we need to estimate
the MSE of EBLUP from observed data in order to measure the uncertainty
of EBLUP. One simple MSE estimator, called the naive MSE estimator, can
be constructed by plugging ÂRE into A in M SEiBLU P :

M̃iN [θ̂i (ÂRE )] = g1i (ÂRE ) + g2i (ÂRE ),

(2)

where M̃iN denotes the naive MSE estimator of EBLUP under the REML
method.
4

However, Kackar and Harville (1984) showed that the MSE of BLUP is
smaller than that of EBLUP because the term depends on the variability
of the estimator for A, which is of the order of O(m−1 ) for large m, and it
is not accurate enough to be ignored for small-area inference (Prasad and
Rao, 1990). The result also implies that the bias of M̃iN (ÂRE ) is of order of
O(m−1 ) under certain regularity conditions. To gain more efficiency even for
such situations, Prasad and Rao (1990) obtained an approximation of true
MSE, M SEi , up to the order O(m−1 ), and the second-order unbiased MSE
estimator, M̂i , of EBLUP with a method of moments estimator of A, so as
to satisfy E[M̂i − M SEi ] = o(m−1 ) for large m. Datta and Lahiri (2000)
and Das et al. (2004) provided such approximation and an MSE estimator
of EBLUP with REML based on the Taylor linearization method:
M SEiRE ≡M SEi [θ̂iEB (ÂRE )] = g1i (A) + g2i (A) + g3i (A) + o(m−1 );
M̂iDL =g1i (ÂRE ) + g2i (ÂRE ) + 2g3i (ÂRE ),

(3)
(4)

where g3i (A) = 2Di2 /[(A + Di )3 tr(V −2 )] and M̂iDL is second-order unbiased
under certain regularity conditions such that E[M̂i − M SEiRE ] = o(m−1 ).
As mentioned above, MSE estimators generally require some bias correction methods to provide second-order unbiasedness.

3

General functional form of MSE estimation for
achieving second-order unbiasedness and strict
positivity

As in Hirose (2016), we consider the general functional form of an MSE
estimator, denoted as:
M̂ig (Âi ) = g1i (Âi ) + g2i (Âi ) + ci (Âi )g3i (Âi ),
with some function ci (A), where Âi is a general adjusted maximum-likelihood
estimator, defined as
Âi = arg max L̃i (A)LRE (A|y),
0≤A<∞

with the general adjustment factor L̃i (A), satisfying Condition A1 given in
the appendix.
We also present a theorem on how to select an adjustment factor, L̃i (A),
for the specified functional form of a second-order unbiased and strictly
5

positive MSE estimator using the adjusted maximum-likelihood method.
This theorem also comes in handy to choose a suitable functional form of
the second-order unbiased and strictly positive MSE estimator for a specified
adjusted maximum-likelihood method.
Theorem 1. Under the regularity conditions and Condition A1, when we
use ci (A) for the adjustment factor L̃i (A), such that
2 − ci (A)
∂ log L̃i (A)
=
+ o(1),
∂A
(A + Di )

(5)

with ci (A) ≤ 2 being of the order of O(1) for large m, satisfying
∂ci (A)
(A + Di ) − ci (A) + 2 ≥ 0,
∂A
the following results hold:
(i) M SEig ≡ M SE[θ̂iEB (Âi )] = g1i (A) + g2i (A) + g3i (A) + o(m−1 ),
= M SEiRE + o(m−1 );
(ii) E[M̂ig (Âi )] = M SEig + o(m−1 );
(iii) With additional condition ci (A) ≥ 0, we obtain M̂ig (Âi ) > 0;
(iv) There exists at least one estimate ÂSi for A > 0, with the conditions
ci (A) ≥ 0 and m > p + 4. A more progressive existence
condition is re-i
h
R (A)
quired for the number of areas m > p such that exp log A2 − ci A
dA =
o(A(m−p)/2 ) holds for large A,

where M̂ig (Âi ) = g1i (Âi ) + g2i (Âi ) + ci (Âi )g3i (Âi ) and
ÂSi = arg max L̃i (A)L̃add (A)LRE (A|y)
0≤A<∞

with L̃add (A) satisfying Condition A2-A3. Incidentally, even if Âi replaces
to ÂSi , parts (i)-(iii) hold. If L̃i (A)
= 0 holds, we no longer need to
A=0

consider the L̃add (A) term.
The proof of Part (i) can be obtained along the same lines as in Das et
al. (2004). The proof of Part (iii) follows from the definition of M̂ig . Parts
(ii) and (iv) are deferred to the appendix. Thus, ci (A) should move only
between 0 and 2.
Corollary 1. From Theorem 1, when the REML estimator is used, the
value 2 is selected as a suitable ci (A), corresponding to M̂iDL .
6

4

Second-order unbiased naive MSE estimator

Theorem 1 ensures that the naive MSE estimator provides second-order
unbiasedness and strict positivity by setting ci (A) = 0. Thus, we obtain
a suitable adjustment factor, L̃N
i (A), up to the order of O(1) for large m,
after solving

∂ log L̃i (A)
∂A

=

2
(A+Di )

+ o(1):

2
L̃N
i (A) = C(A + Di ) ,

where C is a generic positive constant.
However, the estimates could be zero since L̃N
i (A) |A=0 6= 0, as described
in Yoshimori and Lahiri (2014). To avoid this problem, we add an additional adjustment factor, L̃add (A), satisfying Conditions A2-3. For example,
L̃add (A) can be adopted as the specific adjustment factor as in Yoshimori
and Lahiri (2014). Thus, we finally obtain the specific adjusted maximumlikelihood estimator, denoted by ÂN
i , to construct the second-order unbiased
naive MSE estimator while maintaining strict positivity from Theorem 1:
N
ÂN
i = arg max L̃i (A)L̃add (A)LRE (A|y).
0≤A<∞

(6)

Let M̂ N denote the new naive MSE estimator:
N
M̂iN = g1i (ÂN
i ) + g2i (Âi ).

(7)

RE +
Additionally, we also show the result such that M SE i [θ̂i (ÂN
i )] = M SEi
o(m−1 ) from Theorem 1 (i).
Next, we obtain the following theorem on the properties of ÂN
i .

Theorem 2. Under the regularity conditions and Conditions A2-3, we have,
for large m,
(i) E[ÂN
i − A] =

4
tr[V −2 ](A+Di )

2
(ii) E[(ÂN
i − A) ] =

2
tr[V −2 ]

+ o(m−1 ),

+ o(m−1 ),

(iii) ÂN
i is strictly positive for m > p + 4.
The proofs of parts (i) and (ii) are similar to those shown in Yoshimori
and Lahiri (2014). For Part (iii), the proof follows from Theorem 1 (iv) by
setting ci (A) = 0.

7

5

Simulation study

In this section, we compare performances among different estimators of both
the variance parameter A and the MSE of the EBLUP, mentioned in the previous section. In order to investigate the effect of m and Bi , we assume that
m = 15 in a balanced case such that Bi = B patterns: {0.1,0.3,0.5,0.7,0.9}
with fixed Di = D = 1 for all areas. We generated 104 independent data
sets {yi , i = 1, . . . , m} from the Fay–Herriot model (1) with x′i β = 0 and
p = 1. In this simulation study, we also estimate this zero mean from a
practical perspective. In terms of MSE evaluation, we compared the MSE
of the EBLUP with two different estimators: REML ÂRE and our new estimator ÂN . We denote them as “REML” and “NRE”, respectively. When
the REML yielded zero estimates, we treated them as 0.01.
Table 1 shows each simulated MSE of the EBLUP multiplied by 100,
based on two variance estimation methods. From this result, the new variance estimator provides very similar performance to REML in terms of MSE
of EBLUP for small or moderate B values. In contrast, the new variance
estimator does not achieve better performance than the REML method for
large B values in terms of MSE.
Table 1: MSE of EBLUP based on REML and NRE methods, multiplied by
100
Bi = B
0.1
0.3
0.5
0.7
0.9
RE 92.26 77.17 59.83 40.27 20.66
NRE 92.26 77.19 60.97 44.69 28.35
Note: The table shows values increased 100-fold to allow easy comparison.
We also report the percentage of the relative biases (PRB) of different
MSE estimators for the MSE of EBLUP with REML in Table 2 considering
the good performance of the MSE of EBLUP with REML, shown in Table
1. PRB is defined as
P RB :

M̂i − M SEiRE
× 100,
M SEiRE

where M̂i denotes an MSE estimator and M SEiRE is defined as in (3).
We now consider three MSE estimators for M̂i , M̃iN (ÂRE ), M̂iDL (ÂRE ),
and M̂iN (ÂN ), defined in (2), (4), and (7). Hereafter, we denote them as
“Naive.RE”, “DL.RE”, and “Naive.N”.
8

Table 2: Percent RB (PRB) of MSE estimates for MSE of EBLUP with
REML
Bi = B Naive.RE DL.RE Naive.N
0.1
-3.48
-0.08
-0.08
0.3
-12.65
-0.57
-0.63
-21.29
3.99
2.78
0.5
0.7
-22.56
26.27
19.23
-3.29 107.40
75.57
0.9
From the table, the naive estimator with REML (Naive.RE) tends to be
underestimated, unlike other estimators. It probably occurs from the absence of a positive bias correction term to achieve second-order unbiasedness.
As regards other estimation methods, the performance of the second-order
unbiased naive MSE estimator is similar to that of DL.RE for small and
moderate B values. Moreover, our naive estimator, M̂iN (ÂN ), remedies the
over-estimation issue caused by DL.RE for large B values.

6

Conclusion

In this paper, we established that the new estimator ÂN
i conduces to a
second-order unbiased naive MSE estimator while maintaining strict positivity. Results show that the new method remedies the under-estimation issue
associated with the existing naive estimator. Moreover, we also revealed the
relationship between the general functional form of MSE estimation and the
general adjustment factor L̃i (A). Consequently, we can, on the one hand,
easily construct a second-order unbiased and strictly positive MSE estimator
for EBLUP using the specified adjusted maximum-likelihood method and,
on the other, select an adjustment factor with the above MSE estimator.

Acknowledgement
The author’s research was supported by Grant-in-Aid for Research Activity
start-up, JSPS Grant Number 26880011.

9

Appendix.A

Conditions

Regularity Conditions
We assume the following regularity conditions:
R1 rank(X) = p is fixed for large m;
R2 The elements of X are uniformly bounded such that supi≥1 hii = O(m−1 ),
where hii = x′i (X ′ X)−1 xi ;
R3 0 < inf i≥1 Di ≤ supi≥1 Di < ∞, 0 < A < ∞;
R4 |Âi | < Cad mλ , where Cad is a generic positive constant and λ is a
small positive constant, where Âi is a general adjustment maximumlikelihood estimator of A.
We also consider the class of adjustment factors L̃i (A) and L̃add (A) as

Conditions A
A1 log L̃i (A) is independent of y and four times continuously differentiable
with respect to A and a strictly monotonically increasing and concave
L̃i (A)
is of the order of O(1) for large
function of A > 0. Moreover, log∂A
k
m with k = 0, 1, 2, 3, 4;
A2 log L̃add (A) is independent of y and four times continuously differentiable with respect to A. Moreover,
for large m with k = 0, 1, 2, 3, 4;

∂ k log L̃add (A)
∂Ak

is of the order of o(1)

A3 log L̃add (A) is a strictly monotonically increasing and concave function
of A > 0 with L̃add (A)
= 0 and L̃add (A) < C on A > 0 with a
A=0
generic positive constant C.

Appendix.B
Appendix.B.1

Proof of Theorem 1 (ii) and (iv)
Proof of Theorem 1 (ii)

From Theorem1 on Yoshimori and Lahiri (2014), we have for large m
E[g1i (Âi ) − g1i (A)] =Bi2

∂ log L̃i (A)
2
− g3i (A) + o(m−1 ).
−2
tr[V ]
∂A

10

(8)

Using Theorem 1 (i) and the result (8),
E[M̂ig − M SEig ] =E[g1i (Âi ) + g2i (Âi ) + ci (Â)g3i (Âi )]
− [g1i (A) + g2i (A) + g3i (A)] + o(m−1 );
=Bi2

2
∂ log L̃i (A)
+ (ci (A) − 2)g3i (A) + o(m−1 ).
−2
tr[V ]
∂A
(9)

If the first two terms on the right-hand side of (9) vanish for second-order
unbiasedness, we obtain the following differential equation:
[2 − ci (A)]
∂ log L̃i (A)
=
+ o(1).
∂A
A + Di
Thus, Theorem 1 (ii) follows.

Appendix.B.2

Proof of Theorem 1 (iv)

We shall first prove with regard to the progressive condition for ÂG
i existence
on A > 0.
From Conditions A1-A3, we have for A > 0
L̃i (A)L̃add (A)LRE (A)

A=0

= 0, and L̃i (A)L̃add (A)LRE (A)

A=0

> 0.

Thus, it suffices to show that for large A,
L̃i (A)L̃add (A)LRE (A) = o(1).

(10)

Let C denote a generic positive constant. From the fact that
p

m

1

LRE (A) < C(A + sup Di ) 2 |X ′ X|− 2 (A + inf Di )− 2 ,
i≥1

i≥1

Condition A3 reduces (10) to the following, for large A:
L̃i (A) = o(A(m−p)/2 ).
h
R
The solution of the differential equation (5), L̃i (A) = C exp log(A + Di )2 −
further reduces to


Z
ci (A)
2
dA = o(A(m−p)/2 ).
(11)
exp log A −
A
11

i

ci (A)
(A+Di ) dA

,

Thus, we need the number of areas m > p to satisfy the above condition
(11) with fixed p and ci (A) for Conditions A1-A3. When ci (A) ≥ 0 holds,
(11) reduces to


Z
ci (A)
2
dA ≤ A2 .
exp log A −
A
Hence, m > p+4 can be a conservative existence condition for ÂG
i for A > 0.

References
[1] Butar, F. B. and Lahiri, P(2003). On measures of uncertainty of empirical Bayes small area estimators.Journal of Statistical Planning and
Inference,112:63–76
[2] Das, K., Jiang, J., and Rao, J. N. K. (2004). Mean squared error of
empirical predictor. The Annals of Statistics, 32, 818-840.
[3] G.S.Datta, and P.Lahiri, (2000), A unified measure of uncertainty of estimated best linear unbiased predictions in small area estimation problems, Statist.Sinica, 10, 613-627
[4] G.S.Datta, J.N.K.Rao, D.D.Smith, (2005), On measuring the variability of small area estimators under a basic area level model, Biometrika,
92 183-196.
[5] Efron, B. and Morris, C. N. (1975). Data analysis using Stein’s estimator and its generalizations. Journal of the American Statistical Association 70: 311-319.
[6] R.E.Fay, R.A.Herriot, (1979). Estimates of income for small places: an
application of James-Stein procedures to census data, Journal of the
American Statistical Association 74, 269-277
[7] Hall, P. and Maiti, T. (2006). On parametric bootstrap methods for
small area prediction. Journal of the Royal Statistical Society: Series
B, 68, 221-238
[8] Hirose, M. Y. (2016). Non-area-specific adjustment factor for secondorder efficient empirical Bayes confidence interval. arXiv preprint
[9] Kackar, R. N. and 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.
12

[10] Kackar, R. N. and Harville, D. A. (1984). Approximations for standard
errors of estimators of fixed and random effects in mixed linear models.
Journal of the American Statistical Association, 79:853–862
[11] H. Li and P. Lahiri, (2010). An adjusted maximum likelihood method
for solving small area estimation problems, Journal of Multivariate
Analysis, 101, 882-892
[12] N.G.N. Prasad and J.N.K. Rao, (1990). The estimation of the mean
squared error of small area estimators, Journal of the American Statistical Association, 85, 163-171
[13] Yoshimori, M. and Lahiri, P. (2014). A new adjusted maximum likelihood method for the Fay–Herriot small area model. Journal of Multivariate Analysis 124 281-294.

13

