Mean-Squared errors of small area estimators
under a multivariate linear model for repeated
measures data
Innocent Ngaruye, Dietrich von Rosen and Martin Singull
The self-archived postprint version of this journal article is available at Linköping
University Institutional Repository (DiVA):
http://urn.kb.se/resolve?urn=urn:nbn:se:liu:diva-159007
N.B.: When citing this work, cite the original publication.
This is an electronic version of an article published in:
Ngaruye, I., von Rosen, D., Singull, M., (2019), Mean-Squared errors of small area estimators under a multivariate linear model for repeated measures data, Communications in Statistics - Theory and Methods, 48(8), 2060-2073. https://doi.org/10.1080/03610926.2018.1444178
Original publication available at:
https://doi.org/10.1080/03610926.2018.1444178
Copyright: Taylor & Francis (STM, Behavioural Science and Public Health Titles)
http://www.tandf.co.uk/journals/default.asp
Mean-Squared errors of small area estimators under
a multivariate linear model for repeated
measures data
INNOCENT NGARUYE* 1,3, DIETRICH VON ROSEN1,2, and MARTIN SINGULL1 1Department of Mathematics, Linköping University
2Department of Energy and Technology, Swedish University of Agricultural Sciences
3Department of Mathematics, College of Science and Technology, University of Rwanda
Abstract
In this paper, we discuss the derivation of the first and second moments for the proposed small area estimators under a multivariate linear model for repeated measures data. The aim is to use these moments to estimate the mean-squared errors (MSE) for the predicted small area means as a measure of precision. At the first stage, we derive the MSE when the covariance matrices are known. At the second stage, a method based on parametric boot-strap is proposed for bias correction and for prediction error that reflects the uncertainty when the unknown covariance is replaced by its suitable estimator.
Keywords:Mean-squared errors, Multivariate linear model, Parametric bootstrap, Repeated measures data, Small area estimation
Mathematics Subject Classifications: 62F12, 62F40, 62D05
*Address correspondence to Innocent Ngaruye, Department of Mathematics, Linköping University, SE-581 83
1 Introduction
Reliable information about various population characteristics of interest are needed by policy and decision makers for planning. Therefore, there is a great need to estimate these character-istics of interest via survey sampling, not only for the total target population, but also for local sub-population units (domains). However, most sampling surveys are designed to target much larger populations. Then, the derived direct survey estimators obtained using data only from the target small domain of interest have been found to be with lack of precision due to small sample size connected to this domain. The development of estimation techniques that provide reliable estimates for such a small domain or small area and standard errors of estimates have been a big concern in recent years. These techniques are commonly known as Small Area Esti-mation (SAE) methods. For comprehensive reviews of SAE, one can refer to Rao (2003); Rao and Isabel (2015); Pfeffermann (2002, 2013).
Longitudinal surveys with repeated measures data over time are developed to study pattern of changes and trends over time. The demand for SAE statistics is not only for cross-sectional data, but also for repeated measures data. Ngaruye et al. (2017) have proposed a multivariate linear model for repeated measures data within small area estimation settings which accounts for grouped response units and random effects variations. One of the methods to ensure the precision of model-based estimators is the assessment of its mean-squared error (MSE). There is an extensive literature about the approaches used for the estimation of MSE of model-based small area estimators. The second-order approximation of asymptotically unbiased MSE based on Taylor series expansion has been considered by various authors such as Kackar and Harville (1984); Prasad and Rao (1990); Datta et al. (1999); Das et al. (2004); Baillo and Molina (2009), among others. However, as pointed out by Kubokawa and Nagashima (2012), the Taylor se-ries expansion is sometimes complicated to implement for complicated models with many un-known parameters since it requires the computation of asymptotic bias and asymptotic vari-ance and covarivari-ance for estimators of unknown parameters.
unknown parameters in a special case of the model considered by Ngaruye et al. (2017) and we use these moments to derive the MSE for the predicted small area means. Further, following Butar and Lahiri (2003); Kubokawa and Nagashima (2012) we propose an unbiased estimator of MSE based on the parametric bootstrap method.
The article is organized as follows. In Section 2, the description of the considered model is reviewed. In Section 3, the approach used for estimation and prediction is presented. In Section 4, some preliminary basic results that are referred to in the next sections are provided. The first and second moments of the proposed estimators are presented in Section 5 and in Section 6 an unbiased estimator of MSE of predicted small area means under a multivariate linear model for repeated measures data is derived.
2 Description of the model
Consider the multivariate linear regression model for repeated measurements as defined in Ngaruye et al. (2017). Assume that a p-vector of measurements over time for a finite popu-lation of size N divided into m non-overlapping small areas of size Ni, i = 1,...,m together with
r -vector of auxiliary variables are available for all units in the population. Suppose also that the
target population is composed of k group units and denote by Ni g the population size of the
g -th group units, g = 1,...,k such thatPk
g =1Ni g = Ni and that the mean growth of the j th unit,
j = 1,..., Ni, in area i for each group to be a polynomial in time of degree q − 1. Then, the unit
level regression model for the j -th unit coming from the small area i at time t which applies for each one of all k group units can be expressed by
yi j t= β0+ β1t + ··· + βqtq−1+ γ0xi j+ ui t+ ei j t,
i = 1,...,m; j = 1,..., Ni; t = t1, . . . , tp,
where the random errors ei j tand random effects ui tare independent and assumed to be i.i.d.
coefficients representing the effects of auxiliary variables. Theβ0, . . . ,βqare unknown
parame-ters. For all time points, the model can be written in matrix form as
yi j= Aβ + 1pγ0xi j+ ui+ ei j, i = 1,...,m; j = 1,..., Ni;
where 1p is a p-vector of ones and ui is assumed to be multivariate normally distributed with
zero mean and unknown positive definite covariance matrixΣu. In this article we assume A =
Ip, meaning that we do not consider trends over time.
Hence, the associated multivariate linear regression model for all units coming from the i -th small area belonging to the g -th group units can be expressed by
Yi g= βg10Ni g+ 1pγ 0X
i g+ uiz0i g+ Ei g, i = 1,...,m, g = 1,...,k, (1)
and the model at small area level for all k group units together, belonging to the i -th small area can be expressed as
Yi=BCi+ 1pγ0Xi+ uiz0i+ Ei,
ui ∼ Np(0,Σu), Ei ∼ Np,Ni(0,Σe, INi), (2)
where Yi = (Yi 1, ··· ,Yi k); B = (β1, ··· ,βk); Xi = (Xi 1, ··· , Xi k); zi = p1N
i1Ni; Ei = (ei 1, ··· ,ei k); and Ci = blkdiag(10Ni 1, . . . , 10Ni k), where the notation blkdiag(A1, . . . , Ak) is a block diagonal
ma-trix with the given matrices Aion the diagonal. The corresponding model combining all disjoint
m small areas and all N units divided into k non-overlapping group units is given by
Y =BC + 1pγ0X +U Z + E, (3)
where Y = (Y1, ··· ,Ym); C = (10m⊗ Ik)CD= (C1, ··· ,Cm); CD= blkdiag(C1, . . . ,Cm);
X = (X1, ··· , Xm); U = (u1, ··· ,um); E = (E1, ··· ,Em); Z = blkdiag(z01, . . . , z0m) and
E ∼ Np,N(0,Σe, IN), U ∼ Np,m(0,Σu, Im),
with p ≤ m and Σuis an arbitrary positive definite matrix. The symbol ⊗ denotes the Kronecker
an orthogonal transformation and partitioning of the model for estimation purpose.
The matrices Z ,C ,CD are of full row rank,C (Z0) ⊆ C (C0D) and Z Z0= Im, where whereC (A)
denotes the column vector space generated by the columns of an arbitrary matrix A.
In model (3), Y : p ×N is the data matrix, B : q ×k is unknown parameter matrix, CD: mk ×N
with rank(CD)+p ≤ N and p ≤ m is the between individual design matrix accounting for group
effects, the matrix U : p × m is a matrix of random effects, Z : m × N is the design matrix for random effects and E is the error matrix. The matrix C is the between individual design matrix that captures all k group units. More details about model formulation can be found in Ngaruye et al. (2017).
3 Estimation and prediction
Model (3) is considered as a random effects model with covariates. For a comprehensive review of different considerations of the random effects model, see for e.g., Yokoyama and Fujikoshi (1992); Yokoyama (1995); Nummi (1997); Pan et al. (2002). The estimation and prediction are performed with a likelihood based approach. In what follows, for an arbitrary matrix A, Ao stands for any matrix of full rank spanning C (A)⊥, i.e., C (Ao) = C (A)⊥, where C (A)⊥ is an orthogonal complement to C (A). Moreover, A− denotes an arbitrary generalized inverse of the matrix A such that A A−A = A. We also denote by PA = A(A0A)−A0 and QA = I − PA the
orthogonal projection matrices onto the column spaceC (A) and onto its orthogonal comple-mentC (A)⊥, respectively. Derivation of estimators and predictors of model (3) are developed in Ngaruye et al. (2017).
We make an orthogonal transformation of model (3) and partition it into three indepen-dent models. This partition is based on orthogonal diagonalization of the idempotent matrix (CDC0D)−1/2CDZ0Z CD0 (CDC0D)−1/2byΓ = ¡Γ1 Γ2¢, the orthogonal matrix of eigenvectors for m and N − m elements. We use the following notations to shorten matrix expressions:
With this orthogonal transformation, we obtain the three independently distributed models V1=Y R1∼ Np,m ³ B K1+ 1pγ0X R1,Σu+ Σe, Im ´ , V2=Y R2∼ Np,mk−m ³ B K2+ 1pγ0X R2,Σe, Imk−m ´ , V3=Y C0D o ∼ Np,N −mk ³ 1pγ0X C0D o ,Σe, IN −mk ´ .
In the following the estimation of mean and covariance is performed using a likelihood based approach while the prediction of random effect is performed using Henderson’s approach consisting of the maximization of the joint density f (Y ,U ) with respect to U under the assump-tion of knownΣe andΣu (Henderson, 1973). We now have the following theorem for the
esti-mators and predictor.
Theorem 3.1. (Ngaruye et al. (2017)) Consider the model (3). Assume that the matrices X
and Ko20K1are of full rank. Then the estimators forγ, Σu, the linear combination BC and the
predictor of U are given by
b γ =1 p(X P1X 0)−1X P 1Y01p, b BC =³Y −1 p1p1 0 pY P1X0(X P1X0)−1X ´ R2K02(K2K20)−C + W K01P2C , where P1=C0D o (C0Do)0+ R2QK0 2R 0 2, P2=Ko2(Ko20K1K01Ko2)−1K o 20, W =³Y −1 p1p1 0 pY P1X0(X P1X0)−1X ´ GR1, G =IN− R2K02(K2K02)−C . and b Σu= 1 m − 1W QK01Ko2W 0− Σ
e, assumed to be positive definite,
b U =(Σe+ bΣu)−1ΣbuW QK0 1K o 2R 0 1Z0.
For the details of the proof we refer to Ngaruye et al. (2017). Here we have corrected the estimatorΣbuto be an unbiased estimator forΣu. We note that the estimators given in Theorem 3.1 are unique. The following two lemmas discuss the uniqueness of estimatorsBC andb U , theb proof of the uniqueness of other estimators is straightforward.
Lemma 3.1. The estimatorBC given in Theorem 3.1 is invariant with respect to the choice ofb
generalized inverse.
Proof. By replacing W with its expression in the Theorem 3.1, we can rewriteBC asb b BC =³Y − 1 p1p1 0 pY P1X0(X P1X0)−1X ´³ R2K02(K2K20)−+GR1K01P2 ´ C =³Y − 1 p1p1 0 pY P1X0(X P1X0)−1X ´³ R2K02(K2K02)−¡Ik− K1K01P2 ¢ + R1K01P2 ´ C . We can put K1K01P2=K1K01K o 2(K o 20K1K01K o 2)−1K o 20 =Ik− K2 ³ K02(K1K01)−1K2 ´− K02(K1K01)−1. (5) Therefore, b BC =³Y −1 p1p1 0 pY P1X0(X P1X0)−1X ´³ R2K02(K2K02)−K2 ³ K02(K1K01)−1K2 ´− × K02(K1K01)−1+ R1K01P2 ´ C
is unique, which completes the proof of the lemma.
Lemma 3.2. The predictorU given in Theorem 3.1 is invariant with respect to the choice ofb
generalized inverse.
Proof. From the expression ofU in Theorem 3.1, inserting the value ofb Σbuyields b U =(Σe+ bΣu)−1ΣbuW QK0 1Ko2R 0 1Z0= ³ Ip− (m − 1)σ2e¡W QK01Ko2W0 ¢−1´ W QK0 1Ko2R 0 1Z0.
First observe that W QK0 1K o 2 = ³ Y −1 p1p1 0 pY P1X0(X P1X0)−1X ´ R1QK0 1K o 2 −³Y − 1 p1p1 0 pY P1X0(X P1X0)−1X ´ R2K02(K2K02)−K1QK01Ko2
and also note that K01P2K1= PK0
1Ko2. From relation (5), it follows that
K1QK0 1K o 2 =K1− K1PK01K o 2=¡Ik− K1K 0 1P2 ´ K1 =K2 ³ K02(K1K01)−1K2 ´− K02(K1K01)−1K1. Thus, W QK0 1K o
2 does not depend on a choice of a generalized inverse and the uniqueness of
W QK0 1K
o
2 implies the uniqueness ofU .b
3.1 Prediction of small area means
The prediction of small area means is based on the prediction approach to finite population under model-based theory. By this approach, the target population under study is considered as a random sample from a larger population characterized by a suitable model and the pre-dictive distribution of the values for non sampled units is obtained given the realized values of sampled units (Bolfarine and Zacks, 1992). The model (3) considered in this paper belongs to the extensions of unit level model, often known as nested linear regression model which was originally proposed by Battese et al. (1988) for prediction of mean per-capital income in small geographical areas within counties in the United States. Following these authors, the estima-tion of a populaestima-tion mean from the sample returns to the predicestima-tion of a mean of non-sampled values.
For k group units in all small areas, we consider the partition of Ni units into Ni g, g = 1,...,k,
and ni sampled units into ni g such that Ni =Pkg =1Ni g and ni =Pkg =1ni g and similarly for Yi
point for each group unit are given by b µi g = 1 Ni g ³ Y(s)i g1ni g+ bY (r ) i g1Ni g−ni g ´ , (6)
where Y(s)i = (yi 1, ··· , yi ni) : p ×ni, standing for the sampled ni observations from the i -th small area and Yb
(r )
i = (yi ni +1, ··· , yi Ni) : p × (Ni− ni), corresponds to the predicted values for non-sampled (Ni−ni) units from the i -th small area. The first term of the expression (6) on the right
side is known from the sampled observations and the second term is the prediction of non-sampled observations obtained using the considered model (Henderson, 1975) and is given by
b Y(r )i g1Ni g−ni g = (1 − ni g Ni g )bβg+ 1 Ni g 1p10Ni g−ni gX (r )0 i g γ +b pNi g− ni g Ni g b ui, (7)
where X(r )i g stands for the matrix of auxiliary information of non-sampled units in the i -th area belonging to the group units g .
Note that bβg is the estimator ofβg which is the g -th column of the estimated matrixB andb b
ui is the i -th column of the predicted matrixU . Throughout this article, we are interested inb the estimation of mean-squared errors (MSE) for the predicted small area means given in above relation (7). Therefore, we need to calculate the moments of the proposed estimators in order to derive an estimate of the MSE.
4 Preparations for moments calculation
In this section a lemma is given for some technical results that are referred to later on for the moment derivations and other calculations. The following standard definition of the covariance between two random matrices, will be used:
cov(X , Y ) = E[vecX vec0Y ] − E[vecX ]E[vec0Y ],
where vec is the usual columnwise vectorization. Note that the dispersion matrixD[X ] is de-fined byD[X ] = cov(X , X ).
Lemma 4.1. From the definition of the matrices K1, K2, R1and R2given in (4) together with
the fact thatC (Z0) ⊆ C (C0
D) and Z Z0= Im, the following identities hold:
(i) R01R1Z = Z ,
(ii) Z0Z R1= R1, (iii) K1K01= C Z0Z C0,
(iv) R02R1= R02Z0= 0 and hence P1Z0= 0, where P1is defined in Theorem 3.1, and (v) C P1= 0.
Proof. We prove the first three identities, the others are obtained by straightforward
calcula-tions. From the orthogonal diagonalization
(CDC0D)−1/2CDZ0Z C0D(CDC0D)−1/2= ΓDΓ0= µ Γ1 Γ2 ¶ Im 0 0 0 Γ0 1 Γ0 2 . It follows that Γ1(CDC0D)−1/2CDZ0Z C0D(CDC0D)−1/2Γ01= R01Z0Z R1= Im.
Furthermore, since Z Z0 = R01R1 = Im, it follows that Z Z0Z = R01R1Z = Z . Similarly, from R01Z0Z R1 = R01R1, we deduce that Z0Z R1 = R1. Moreover, since K1 is a full row rank, then K1K01= K1R10Z0Z R1K01= C R1R01Z0Z R1R01C0= C Z0Z C0.
5 Moments of proposed estimators
In this section, we present the first and second moments of the proposed estimators. All proofs will be given in the end of this article in an Appendix.
5.1 Theoretical moments with known
Σu
For the purpose of moment calculations, from now on, we suppose that we have complete knowledge on the covariance matrixΣuor it has been estimated previously so that it is taken to
be as known. In that regard, we use the following predictor of U
e
U = (Σe+ Σu)−1ΣuW QK10Ko2R01Z0. (8)
Theorem 5.1. Given the estimators in Theorem 3.1. Then,Y = bb BC + 1pγb0X + eU Z is an
unbi-ased predictor, i.e.,E[Y ] = E[Y ].b
The following two theorems give the main results of the paper about moments of the pro-posed estimators.
Theorem 5.2. Given the estimators in Theorem 3.1. The dispersion matricesD[bγ], D[BC ] andb D[U ] are given bye D[bγ] =σ2e p (X P1X 0)−1, D[BC ] =Cb 0P2C ⊗ ¡ Σu+ Σe¢ + ³ C0¡Ik− P2K1K01¢(K2K02)− סIk− K1K01P2¢C ´ ⊗ Σe+ σ2 e p C 0³(K 2K02)−K2R20 + P2K1R01G0 ´ × X0(X P1X0)−1X ³ R2K02(K2K20)−+GR1K01P2 ´ C ⊗ 1p10p, D[U ] =Z Re 1QK0 1K o 2¡Im+ K1(K2K 0 2)−K1¢QK0 1K o 2R 0 1Z0⊗ Σu(Σe+ Σu)−1Σu +σ 2 e p Z R1QK01K o 2R 0 1G0X0(X P1X0)−1X GR1QK0 1K o 2R 0 1Z0 ⊗ Σu(Σe+ Σu)−11p10p(Σe+ Σu)−1Σu.
In Theorem 5.2 the dispersion matrices for the estimators and predictor are given. However, for prediction purposes, it is also of interest to derive the covariances between them.
cov[BC ,b U ] and cov[e BC ,b γ] are given byb cov[U ,e γ] = −b σ2 e p Z R1QK01K o 2R 0 1G0X0(X P1X0)−1⊗ Σu(Σe+ Σu)−11p, cov[BC ,b U ] =e σ2 e p C 0¡(K 2K02)−K2R02+ P2K1R01G0¢X0(X P1X0)−1X GR1QK0 1K o 2R 0 1Z0 ⊗ 1p10p(Σe+ Σu)−1Σu, cov[BC ,b γ] = −b σ2 e p C 0¡(K 2K02)−K2R02+ P2K1R01G0¢X0(X P1X0)−1⊗ 1p.
In the next section we will use Theorem 5.2 and Theorem 5.3 to derive the mean-squared errors of predicted small area means.
5.2 Simulation study with the empirical moments
To put some light on the moment expressions in Section 5.1 we provide a simulation study comparing the first and second moments ofB andb γ. Assume we have 8 small areas, i.e, m = 8b with k = 3 groups and given sample sizes
n11= 2, n12= 3, n13= 6, n21= 7, n22= 9, n23= 13, n31= 10, n32= 4, n33= 5, n41= 3, n42= 9, n43= 7, n51= 4, n52= 6, n53= 2, n61= 10, n62= 11, n63= 18, n71= 7, n72= 4, n73= 4, n81= 6, n82= 9, n83= 5. Furthermore, let p = q = 2 with
B = 8 10 12 9 11 13 , γ = µ 1 2 3 ¶0 , Σu= 2 1 1 2 ,
Σe= Ip, i.e.,σ2e= 1, and the design matrix X is chosen randomly as xi j∼ U [0, 1]. Given model
sam-ple sizes above with the factors 1, 2, 4, 8, 12, 16, 20, i.e., we have n11= 2, 4, 8, 16, 24, 32, 40, and so on. For 10 000 replicates, for each set up of sample sizes, we calculate the estimatesB andb b
γ, respectively. Based on these simulated estimates we derive the empirical moments and the
Frobenius norm of the difference between the empirical (bE[·] andD[·]) and true moments (E[·]b andD[·]) given in Section 5.1. As we can see in Table 1, the estimates are close to the true values and getting better when the sample sizes increases (i.e., estimators seems to be consistent).
Factor ||B − bE[B ]||b F ||D[ bB ] − bD[B ]||b F ||γ − bE[bγ]||F ||D[γ] −b D[bb γ]||F 1 0.0106 0.0176 0.0072 0.0039 2 0.0202 0.0037 0.0074 0.0025 4 0.0052 0.0026 0.0045 0.0008 8 0.0084 0.0012 0.0038 0.0005 12 0.0035 0.0005 0.0034 0.0002 16 0.0036 0.0005 0.0020 0.0002 20 0.0031 0.0005 0.0016 0.0001
Table 1: Frobenius norm of the difference between the empirical (bE[·] andD[·]) and true mo-b ments (E[·] and D[·]).
6 Mean-squared errors of predicted small area means
6.1 Derivation of MSE(
τ
b
i g)
In this section, following Butar and Lahiri (2003) we estimate the mean squared error for pre-dicted small area means in two steps. At the first step, under the assumption of known covari-ance matricesΣeandΣu, the derivation of MSE is presented. At the second step, a parametric
bootstrap approach is proposed for bias correction and approximation of the uncertainty due to the estimation ofΣu. Put
K =³1 −ni g Ni g ´ I , L = 1 Ni g 1p10Ni g−ni gX (r )0 i g , M = pNi g− ni g Ni g Ip, (9) i = 1,...,m, g = 1,...,k.
Then, the linear prediction and empirical linear prediction quantities from small area means given in (7) can be written by
e
τi g=K bβg+ Lγ + Mb uei, i = 1,...,m, g = 1,...,k
b
τi g=K bβg+ Lγ + Mb ubi, i = 1,...,m, g = 1,...,k. (10)
Let eg : k × 1 and fi : m × 1 be the unit basis vectors, i.e., k and m vectors with 1 in the g th and
i th position, respectively, and 0 elsewhere. Putαg = C0(CC0)−1eg so that bβg = bBCαg,uei= eU fi andubi= bU fi.
Furthermore, let us write τbi g− τi g = (bτi g−τei g) + (eτi g − τi g) and the MSE of τbi g can be obtained by MSE(bτi g) =MSE(eτi g) + 2E £¡ b τi g−eτi g ¢¡ e τi g− τi g ¢0 ¤ + E£¡bτi g−eτi g ¢¡ b τi g−τei g ¢0 ¤. (11)
The first term of the right hand side of (11) has the form
MSE(τei g) =E[(eτi g− τi g)(τei g− τi g) 0] =E[(K (bβg− βg) + L(γ − γ) + M(b uei− ui)) × (K (bβg− βg) + L(bγ − γ) + M(uei− ui)) 0] =K D[bβg]K0+ K cov[bβg,γ]Lb 0+ K cov[bβ g,uei]M 0+ LD[bγ]L0 + Lcov[bγ,uei]M 0+ MD[eu i, bβg]K0+ Lcov[bγ,B ]Kb 0 + Mcov[eui,γ]Lb 0+ MD[eu i− ui]M0. (12)
Observe that, from the definitions given to bβg= bBCαg,uei= eU fi andubi = bU fi the covariances presented in equation (12) are expressed by
cov[bβg,uei] = cov[ bBCαg,U fe i] = (α 0 g⊗ Ip)cov[bBC ,U ]( fe i⊗ Ip), cov[uei,γ] = cov[b U fe i,γ] = (fb 0 i⊗ Ip)cov[U ,e γ],b cov[bβg,γ] = cov[b BCb αg,γ] = (αb 0 g⊗ Ip)cov[BC ,b γ].b
Similarly, the dispersion matrices presented in equation (12) are expressed by D[βbg] = D[ bBCαg] = (α0g⊗ Ip)D[BC ](b αg⊗ Ip), D[eui] = D[ eU fi] = (f0i⊗ Ip)D[U ]( fe i⊗ Ip), and D[eui− ui] =(f0i⊗ Ip)D[U −U ](fe i⊗ Ip) =( f0i⊗ Ip)(D[U ] − 2cov( ee U ,U ) + D[U ])(fi⊗ Ip) =( f0i⊗ Ip)(D[U] − D[U ])( fe i⊗ Ip).
Altogether we can give the following theorem.
Theorem 6.1. The MSE of the linear prediction from small area meansτei g is given as
MSE(τei g) =K (α 0 g⊗ Ip)D[BC ](b αg⊗ Ip)K0+ K (α0g⊗ Ip)cov[BC ,b γ]Lb 0 + K (α0g⊗ Ip)cov[BC ,b U ]( fe i⊗ Ip)M0+ LD[bγ]L0 + Lcov[ eU ,γ]b 0( f i⊗ Ip)M0+ M( f0i⊗ Ip)cov[bBC ,U ]e 0(αg⊗ Ip)K0 + Lcov[ bBC ,γ]b 0(α g⊗ Ip)K0+ M( f0i⊗ Ip)cov[U ,e γ]Lb 0 + MΣuM0− M( f0i⊗ Ip)D[U ]( fe i⊗ Ip)M0, (13)
where the dispersion matricesD[γ], D[b BC ],b D[U ] and the covariance matrices cov[e U ,e γ], cov[b BC ,b U ],e cov[BC ,b γ] are presented in Theorem 5.2 and Theorem 5.3, respectively.b
6.2 Estimation of MSE(
b
τ
i g)
The second two terms of the right hand side of equation (11) are intractable and need to be approximated. It is important to note that in practice, the covariance matrixΣu is unknown.
As pointed out by different authors (see for example Das et al. (2004)), a naive estimator of MSE obtained by replacing the unknown covariance matrixΣu by its estimatorΣbuin (13) that ignores the variability associated withΣbu can lead to underestimation of the true MSE.
There-fore, we propose a parametric bootstrap method to estimate the MSE whenΣu is replaced by
its estimator.
The approximate estimator of MSE(bτi g) given in equation (11) can be decomposed as
MSE(τbi g) = G1i(Σu) +G2i(Σu) +G3i(Σu), (14) where G1i(Σu) +G2i(Σu) = MSE(eτi g) G3i(Σu) = 2E £¡ b τi g−eτi g ¢¡ e τi g− τi g ¢0 ¤ + E£¡bτi g−eτi g ¢¡ b τi g−eτi g ¢0 ¤.
For knownΣu, the quantity G1i(Σu) + G2i(Σu) is given in equation (13). WhenΣu is replaced
by its estimator, the quantity G1i(bΣu) + G2i(Σbu) introduces an additional bias related to Σbu, i.e.,E[Σbu] − Σu(see Datta and Lahiri (2000)). Following Butar and Lahiri (2003); Kubokawa and Nagashima (2012), we propose a parametric bootstrap method to estimate the first two terms of (14) by correcting the bias of G1i(bΣu) + G2i(bΣu) whenΣu is replaced by its estimatorΣbu in equation (13) and secondly for estimating the third term G3i(Σu) of equation (14).
Consider the bootstrap model
Y∗i | u∗i ∼ Np,ni ¡ b BCi+ 1pγb 0X i+ u∗iz0i,Σe, Ini¢, i = 1,...,m, (15) where u∗i ∼ Np(0,Σbu). If we put b τi g(Yi; bβg,γ,b Σbu) =K bβg+ Lbγ + Muei, (16) b τi g(Yi; bβ∗g,γb ∗, b Σu) =K bβ∗g+ Lbγ∗+ Mue ∗ i, (17) b τi g(Yi; bβ∗g,γb ∗, b Σ∗u) =K bβ ∗ g+ Lbγ∗+ Mub ∗ i, i = 1,...,m, g = 1,...,k. (18)
then the quantity G1i(Σu) +G2i(Σu) can be estimated by
G1i(bΣu) +G2i(bΣu) − E∗£G1i(bΣu∗) +G2i(bΣ∗u) −G1i(bΣu) −G2i(bΣu)
and the quantity G3i(Σu) estimated by 2E∗h³τbi g(Yi; bβ∗g,γb ∗, b Σ∗u) −τbi g(Y ∗ i; bβ ∗ g,γb ∗, b Σu) ´³ b τi g(Yi; bβ∗g,γb ∗, b Σu) −bτi g(Yi; bβg,γ,b Σbu) ´0i + E∗ h³ b τi g(Yi; bβ∗g,γb ∗, b Σ∗u) −bτi g(Yi; bβ ∗ g,γb ∗, b Σu) ´³´0i ,
whereE∗ is the expectation with respect to model (15), and the calculation ofγb ∗, bβ∗
g, Σb ∗
u is
performed similarly to that ofγ,b βbg, Σbu except thatγb∗, bβ ∗
g, Σb ∗
u are calculated based on Y∗i
instead of Yi.
Thus, when the unknown covarianceΣuis replaced by its estimator in MSE(τbi g) of equation (14), we obtain the proposed estimator of MSE whose results are summarized in the following theorem.
Theorem 6.2. Consider the model (3) and assume that the matrices X and Ko20K1are of full
rank. Consider the bootstrap model (15) and the estimated quantities (16)-(18). Then the estima-tor of MSE for predicted small area means given in (10) can be expressed by
MSE(bτi g) =2 h G1i(bΣu) +G2i(bΣu) i − E∗ h G1i(bΣ∗u) +G2i(Σb ∗ u) i + 2E∗ h³ b τi g(Yi; bβ∗g,γb ∗, b Σ∗u) −τbi g(Y ∗ i; bβ ∗ g,γb ∗, b Σu) ´ × ³ b τi g(Yi; bβ∗g,γb ∗, b Σu) −τbi g(Yi; bβg,γ,b Σbu) ´0i + E∗ h³ b τi g(Yi; bβ∗g,γb ∗, b Σ∗u) −bτi g(Yi; bβ ∗ g,γb ∗, b Σu) ´³´0i , (19)
where G1i(bΣu)+G2i(bΣu) and G1i(bΣ∗u)+G2i(bΣ∗u) are given by equation (13) withΣureplaced byΣbu
andΣb ∗
u, respectively. In addition, the dispersion matricesD[γ], D[b BC ],b D[U ] and the covariancee
matrices cov[U ,e γ], cov[b BC ,b U ], cov[e BC ,b γ] involved in equation (13) are presented in Theoremb
5.2 and Theorem 5.3, respectively.
Appendix
Proof of Theorem 5.1. First we show thatγ is an unbiased estimator. We haveb E[bγ] =1 p(X P1X 0)−1X P 1E[Y0]1p = 1 p(X P1X 0)−1X P 1¡C0B0+ X0γ10p¢1p =1 p(X P1X 0)−1X P 1X0γ10p1p = γ. Therefore,E[1pγb 0X ] = 1
pγ0X . Moreover, from the expression ofBC given in Theorem 3.1 andb (5), it follows that E[BC ] =b ³ E[Y ] −1 p1p1 0 pE[Y ]P1X0(X P1X0)−1X ´³ R2K02(K2K20)−+GR1K01P2 ´ C = ³ BC + 1pγ0X − 1 p1p1 0 p(BC + 1pγ0X )P1X0(X P1X0)−1X ´ × ³ R2K02(K2K20)−+GR1K01P2 ´ C =BC¡R2K02(K2K20)−+GR1K01P2¢C =B ³ K2K02(K2K20)−+ K1K01P2− K2K02(K2K02)−K1K01P2 ´ C =B ³ K2K02(K2K02)−+ Ik− K2 ³ K02(K1K01)−1K2 ´− K02(K1K01)−1 − K2K02(K2K02)−+ K2K02(K2K02)−K2 ³ K02(K1K01)−1K2 ´− K02(K1K01)−1 ´ C =BC .
Hence, alsoBC is an unbiased estimator. Finally, we have forb U given in (8) the meane E[U Z ] =(Σe e+ Σu)−1ΣuE[W ]QK0 1K o 2R 0 1Z0Z =(Σe+ Σu)−1ΣuB ³ Ik− K2K02(K2K02)− ´ K1QK0 1Ko2R 0 1Z0Z , (20)
since from the expression of W , it follows that E[W ] = BCGR1= B
³
Ik− K2K02(K2K02)− ´
K1.
Recall that from the proof of Lemma 3.2, we have
K1QK0 1K o 2= K2 ³ K02(K1K01)−1K2 ´− K02(K1K01)−1K1
and plugin this expression intoE[U Z ], given in (20), leads toe E[U Z ] =(Σe e+ Σu)−1ΣuB ³ Ik− K2K02(K2K02)− ´ × K2 ³ K02(K1K01)−1K2 ´− K02(K1K01)−1K1R01Z0Z = 0.
Hence, it is straightforward to see that
E[Y ] = E[ bb BC + 1pγb0X + eU Z ] = BC + 1pγ0X = E[Y ],
which completes the proof of the theorem.
In what follows, we use the notation (A)Q()0instead of (A)Q(A)0when it is possible and no confusions, in order to shorten the matrix expressions.
Proof of Theorem 5.2. Recall thatγ =b p1(X P1X0)−1X P1Y01p. Hence, it follows that
D[bγ] = 1 p2 ³ 10p⊗ (X P1X0)−1X P1 ´ D[Y0]³´0 = 1 p2 ³ 10p⊗ (X P1X0)−1X P1 ´³ Σu⊗ Z0Z + Σe⊗ IN ´³´0 =σ 2 e p ⊗ ³ (X P1X0)−1X P1X0(X P1X0)−1 ´ =σ 2 e p (X P1X 0)−1,
since P1Z0= 0. Moreover, replacing W by its expression in the Theorem 3.1, we can rewrite bBC by b BC =³Y −1 p1p1 0 pY P1X0(X P1X0)−1X ´³ R2K02(K2K20)−+GR1K01P2 ´ C
and therefore, the dispersion matrixD[BC ] is given byb D[BC ] =b ³ C0¡(K2K02)−K2R02+ P2K1R01G0¢ ⊗ Ip −1 pC 0¡(K 2K02)−K2R02+ P2K1R10G0¢X0(X P1X0)−1X P1⊗ 1p10p ´ × ³ Z0Z ⊗ Σu+ IN⊗ Σe ´³´0 .
Using the results from Lemma 4.1, we obtain D[BC ] =Cb 0P2C ⊗ ¡ Σu+ Σe¢ + ³ C0¡Ik− P2K1K01¢(K2K02)− סIk− K1K01P2¢C ´ ⊗ Σe+ σ2 e p C 0³(K 2K02)−K2R20 + P2K1R01G0 ´ × X0(X P1X0)−1X ³ R2K02(K2K20)−+GR1K01P2 ´ C ⊗ 1p10p. Moreover, fromU = Σe u(Σe+ Σu)−1W QK0 1Ko2R 0
1Z0, and replacing W by its expression, we can rewrite e U = Σu(Σe+ Σu)−1 ³ Y −1 p1p1 0 pY P1X0(X P1X0)−1X ´ GR1QK0 1K o 2R 0 1Z0. Given this, the dispersion matrixD[U ] has the forme
D[U ] =e ³ Z R1QK0 1Ko2R 0 1G0⊗ Σu(Σe+ Σu)−1− 1 pZ R1QK01Ko2R 0 1G0X0(X P1X0)−1X P1 ⊗ Σu(Σe+ Σu)−11p10p ´³ Z0Z ⊗ Σu+ IN⊗ Σe ´³´0 =Z R1QK0 1K o 2¡Im+ K1(K2K 0 2)−K1¢QK0 1K o 2R 0 1Z0⊗ Σu(Σe+ Σu)−1Σu +σ 2 e p Z R1QK01K o 2R 0 1G0X0(X P1X0)−1X GR1QK0 1K o 2R 0 1Z0 ⊗ Σu(Σe+ Σu)−11p10p(Σe+ Σu)−1Σu.
Proof of Theorem 5.3. By using the results from Lemma 4.1, we get
cov[U ,e γ] =b ³ Z R1QK0 1Ko2R 0 1G0⊗ Σu(Σe+ Σu)−1− 1 pZ R1QK01Ko2R 0 1G0X0(X P1X0)−1X P1 ⊗ Σu(Σe+ Σu)−11p10p ´³ Z0Z ⊗ Σu+ IN⊗ Σe ´³1 pP1X 0(X P 1X0)−1⊗ 1p ´ = −σ 2 e p Z R1QK01Ko2R 0 1G0X0(X P1X0)−1⊗ Σu(Σe+ Σu)−11p,
and cov[BC ,b U ] =e ³ C0¡(K2K02)−K2R02+ P2K1R01G0¢ ⊗ Ip − 1 pC 0¡(K 2K02)−K2R02+ P2K1R10G0¢X0(X P1X0)−1X P1⊗ 1p10p ´ × ³ Z0Z ⊗ Σu+ IN⊗ Σe ´³ GR1QK0 1K o 2R 0 1Z0⊗ Σu(Σe+ Σu)−1 − 1 pP1X 0(X P 1X0)−1X GR1QK0 1Ko2R 0 1Z0⊗ 1p10p(Σe+ Σu)−1Σu ´ =σ 2 e p C 0¡(K 2K02)−K2R02+ P2K1R01G0¢X0(X P1X0)−1X GR1QK0 1K o 2R 0 1Z0 ⊗ 1p10p(Σe+ Σu)−1Σu.
Similarly we have the last covariance as
cov[BC ,b γ] =b ³ C0¡(K2K02)−K2R02+ P2K1R01G0¢ ⊗ Ip −1 pC 0¡(K 2K02)−K2R02+ P2K1R10G0¢X0(X P1X0)−1X P1⊗ 1p10p ´ ׳Z0Z ⊗ Σu+ IN⊗ Σe ´³1 pP1X 0(X P 1X0)−1⊗ 1p ´ = −σ 2 e p C 0¡(K 2K02)−K2R02+ P2K1R01G0¢X0(X P1X0)−1⊗ 1p.
References
Baillo, A. and Molina, I. (2009). Mean-squared errors of small-area estimators under a unit-level multivariate model. Statistics, 43(6):553–569.
Battese, G. E., Harter, R. M., and Fuller, W. A. (1988). An error-components model for predic-tion of county crop areas using survey and satellite data. Journal of the American Statistical
Association, 83(401):28–36.
Butar, F. B. and Lahiri, P. (2003). On measures of uncertainty of empirical bayes small-area estimators. Journal of Statistical Planning and Inference, 112(1):63–76.
Das, K., Jiang, J., and Rao, J. (2004). Mean squared error of empirical predictor. The Annals of
Statistics, 32(2):818–840.
Datta, G. S. and Lahiri, P. (2000). A unified measure of uncertainty of estimated best linear unbiased predictors in small area estimation problems. Statistica Sinica, pages 613–627.
Datta, G. S., Lahiri, P., Maiti, T., and Lu, K. L. (1999). Hierarchical Bayes estimation of un-employment rates for the states of the US. Journal of the American Statistical Association, 94(448):1074–1082.
Henderson, C. R. (1973). Sire evaluation and genetic trends. Journal of Animal Science,
1973(Symposium):10–41.
Henderson, C. R. (1975). Best linear unbiased estimation and prediction under a selection model. Biometrics, 31(2):423–447.
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(388):853–862.
Kubokawa, T. and Nagashima, B. (2012). Parametric bootstrap methods for bias correction in linear mixed models. Journal of Multivariate Analysis, 106:1–16.
Ngaruye, I., Nzabanita, J., von Rosen, D., and Singull, M. (2017). Small area estimation under a multivariate linear model for repeated measures data. Communications in Statistics - Theory
and Methods, 46(21):10835–10850.
Nummi, T. (1997). Estimation in a random effects growth curve model. Journal of Applied
Pan, J.-X., Fang, K., and Fang, K.-T. (2002). Growth Curve Models and Statistical Diagnostics. Springer Science & Business Media.
Pfeffermann, D. (2002). Small area estimation - new developments and directions.
Interna-tional Statistical Review, 70(1):125–143.
Pfeffermann, D. (2013). New important developments in small area estimation. Statistical
Sci-ence, 28(1):40–68.
Prasad, N. and Rao, J. (1990). The estimation of the mean squared error of small-area estimators.
Journal of the American Statistical Association, 85(409):163–171.
Rao, J. N. K. (2003). Small Area Estimation. John Wiley and Sons, New York.
Rao, J. N. K. and Isabel, M. (2015). Small Area Estimation. John Wiley and Sons, New York, 2 edition.
Yokoyama, T. (1995). Statistical inference on some mixed manova-gmanova models with ran-dom effects. Hiroshima Mathematical Journal, 25(3):441–474.
Yokoyama, T. and Fujikoshi, Y. (1992). Tests for random-effects covariance structures in the growth curve model with covariates. Hiroshima Mathematical Journal, 22:195–202.
