Stat Papers DOI 10.1007/s00362-013-0517-y REGULAR ARTICLE

Spatial robust small area estimation Timo Schmid · Ralf T. Münnich

Received: 29 May 2012 / Revised: 27 February 2013 © Springer-Verlag Berlin Heidelberg 2013

Abstract The accuracy of recent applications in small area statistics in many cases highly depends on the assumed properties of the underlying models and the availability of micro information. In finite population sampling, small sample sizes may increase the sensitivity of the modeling with respect to single units. In these cases, area-specific sample sizes tend to be small such that normal assumptions, even of area means, seem to be violated. Hence, applying robust estimation methods is expected to yield more reliable results. In general, two robust small area methods are applied, the robust EBLUP and the M-quantile method. Additionally, the use of adequate auxiliary information may further increase the accuracy of the estimates. In prediction based approaches where information is needed on universe level, in general, only few variables are available which can be used for modeling. In addition to variables from the dataset, in many cases further information may be available, e.g. geographical information which could indicate spatial dependencies between neighboring areas. This spatial information can be included in the modeling using spatially correlated area effects. Within the paper the classical robust EBLUP is extended to cover spatial area effects via a simultaneous autoregressive model. The performance of the different estimators are compared in a model-based simulation study. Keywords

Small area estimation · M-quantile · Robust EBLUP · Spatial correlation

T. Schmid (B) Institute of Statistics and Econometrics, Freie Universität Berlin, Garystr. 21, 14195 Berlin, Germany e-mail: [email protected] R. T. Münnich Economic and Social Statistics Department, University of Trier, Universitätsring 15, 54286 Trier, Germany e-mail: [email protected]

123

T. Schmid and R. T. Münnich

Mathematics subject classification (2000)

62F10 · 62F35

1 Introduction Economic and political decision processes are increasingly based on specific indicators and other statistical information. Nowadays the necessity of developing regional indicator values or disaggregated values is evident in order to allow for regional or group-specific comparisons. Surveys which shall deliver the necessary information for these indicators, however, are generally constructed for larger areas, e.g. national or regional level (NUTS—Nomenclature des unités territoriales statistiques—1–2 levels). Hence, sample information on lower levels such as municipalities and provinces (LAU—Local Administrative Unit—1–2 levels) is rarely available which leads to high variances of classical estimates. Examples of interest are manifold. In 2011, the last European Census round took place. Several countries, such as Germany and Switzerland, used sample information in addition to register information. Finally, accurate numbers on NUTS3 or even community level may have to be derived (cf. Münnich et al. 2007). Within the seventh Framework Programme, two projects investigate poverty measures for small areas (AMELI, http://ameli.surveystatistics.net, and SAMPLE, http://www.sample-project. eu/). The aim is to develop a methodology which enables high quality estimates of measures for poverty and social exclusion for small areas and smaller sub-populations. And finally, reliable estimates on profits, turnovers, or even performance indicators are needed for business statistics. A challenging problem is to deliver high quality estimates on Nomenclature statistique des activités économiques dans la Communauté européenne (NACE) classes while having highly dominant enterprises (cf. BLUE-ETS project, http://www.blue-ets.eu). In all cases, classical estimates of smaller areas or sub-populations suffer from possibly very small subgroup-specific sample sizes which leads to unacceptably large variances. In the given cases, applying small area estimation methods may lead to highly improved accuracy of the estimates of interest. In poverty measurement and business statistics, variables of interest tend to be very skewed. Further, and especially in business statistics, outliers and even extreme outliers are observed. These variables in connection with small sample sizes lead to severe problems while applying standard models which are based on normal assumptions due to the high sensitivity of the model estimates towards these influential units. One way to overcome these peculiarities is the application of robust methods. Chambers (1986) and Welsh and Ronchetti (1998) studied robust estimation of finite population parameters in the presence of outliers whereas Fellner (1986), Richardson and Welsh (1995) and Richardson (1997) investigated outlier-robust inference for the linear mixed model. Afterwards, Sinha and Rao (2009) extended this approach to robust small area estimation. This proposed robust estimator is an extension of the classical EBLUP estimator introduced, e.g. in Rao (2003) or Jiang and Lahiri (2006). Another method for robust small area estimation was presented by Chambers and Tzavidis (2006), using an approach based on fitting outlier robust M-quantile models

123

Spatial robust small area estimation

to survey data. Both approaches use robust prediction where they replace parameter estimates in optimal, but outlier sensitive, estimators by outlier robust versions. The reliability of the above methods is generally based on the availability of accurate auxiliary information. In prediction approaches in small area statistics, this information has to be available on universe level which generally limits the number of available variables. In business data, the covariates are hardly available on individual level. Furthermore, only listed enterprises which excess specific thresholds by number of employees or turnover may be obligated to publish their results. In many cases, however, additional information such as distances between areas, similarities of groups may be available such that spatial information can be used in the models. The spatial correlation can be accounted for by extending the random effects model to allow for spatially correlated area effects. In this context, Singh et al. (2005) and Pratesi and Salvati (2008) proposed a spatial extension of the EBLUP estimator. Later, Salvati et al. (2012) introduced a spatial version of the M-quantile models based on geographically weighted regression. In the second Section, robust small area estimation methods are introduced, i.e. the robust EBLUP and the M-quantile method. Within the third Section, a spatial extension of the robust EBLUP is proposed (SREBLUP). The performance of the SREBLUP in contrast to the estimators in Sect. 2 is shown in a small model-based simulation study in Sect. 4.

2 Robust small area methods This section provides an overview regarding robust small area methods. The first part focuses on the robust EBLUP (REBLUP) introduced by Sinha and Rao (2009) whereas the second part describes an alternative approach of Chambers and Tzavidis (2006) to small area estimation based on M-quantile methods.

2.1 EBLUP-based small area methods We start from the general linear mixed model, which is given by y = Xβ + Zv + ε,

(1)

v ∼ N (0, G), ε ∼ N (0, R), where y is the n × 1 vector of a dependent variable, X the n × k matrix of known auxiliary variables with associated k × 1 regression vector β. Z denotes the n × q design matrix for the random effects. The q × 1 vector of random effects v and the n × 1 vector of random errors ε are assumed to be independently normally distributed with means 0 and covariance matrices G and R respectively depending on a vector of variance parameters θ , e.g. θ = {σv2 , σe2 }. Furthermore, the variance-covariance matrix of the variable of interest y is obtained via V = R + ZG Z T .

123

T. Schmid and R. T. Münnich

In the context of small area estimation, the parameters of interest are for instance the population mean or the population total of the dependent variable in specific areas (cf. Rao 2003, Chapt. 2.2). These quantities can be described as a linear combination μ(θ ) = l T β + m T v of the regression coefficient β and the random effect v for specific vectors l and m. The vector of variance components is in practice also unknown. Therefore, it is necessary to replace θ by a consistent estimator θˆ to obtain a two-stage estimator μ( ˆ θˆ ). The resulting estimator μ( ˆ θˆ ) = l T βˆ + m T vˆ ,

(2)

ˆ θˆ ) is the is referred to the empirical best linear unbiased predictors (EBLUP). βˆ = β( BLUE estimator of β and vˆ = vˆ (θˆ ) denotes the BLUP estimator of v. In the case of unit level models which was derived by Battese et al. (1988), the EBLUP for instance for the mean y i in area i is given by yˆ i = Ni−1

⎧ ⎨ ⎩

j∈si

yi j +

 j∈ri

⎫ ⎬

(x iTj βˆ + z iTj vˆi ) , ⎭

(3)

where si and ri denote the index sets of the sampled and non-sampled quantities respectively in the area i. More details regarding variance estimation methods and the EBLUP can be obtained from Harville (1977), Richardson and Welsh (1994), Rao (2003), Jiang and Lahiri (2006), Zadlo (2009) or Niemiro and Wesolowski (2010). From a robust viewpoint, the proposed estimators (2) or (3) can be very sensitive to outliers in the data. Therefore, it is useful to substitute the estimators βˆ and vˆ by robust alternatives. Sinha and Rao (2009) maximized the density of the random variable y with respect to β and θ and obtained outlier robust ML-estimators of β and θ by solving the following robustified ML-equations 1

α(β) = X T V −1 U 2 ψ(r ) = 0

1 1 (θl ) = ψ T (r )U 2 V −1 ∂∂θVl V −1 U 2 ψ(r ) − tr V −1 ∂∂θVl K = 0,

(4) (5)

where ψ(u) = (ψb (u 1 ), ψb (u 2 ), ...)T

(6)

is the Huber function with ψb (u) = u min(1, b/|u|) and b > 0 is a tuning constant which is generally set to 1.345 (cf. Hampel et al. 1986 or Dasiou and Moyssiadis 2001). θl is the lth component of the variance parameter θ = (θ1 , . . . , θq )T and r is given by 1 r = U − 2 (y−Xβ). U is a diagonal matrix with diagonal elements equal to the diagonal elements of the matrix V . K is also a diagonal matrix with K = E(ψb2 (t))Idn×n , where t follows a standard normal distribution and Idn×n denotes the identity matrix with dimension n. Sinha and Rao (2009) used the Newton-Raphson algorithm to

123

Spatial robust small area estimation

solve the robust equations and defined the estimates at convergence to be the robust ψ ψ ML-estimators βˆ and θˆ of β and θ respectively, where the superscript ψ stands for influence function ψ (e.g. Huber function). ψ ψ Sinha and Rao (2009) then used βˆ and θˆ for the estimation of the robust areaspecific random effects v by plugging them into Fellner’s equation (cf. Fellner 1986). ψ ψ The robust estimators βˆ , θˆ and vˆ ψ are then used to predict (2) as ψ

ψ

μ( ˆ θˆ ) = l T βˆ + m T vˆ ψ

(7)

ψ for specified vectors l and m. Sinha and Rao (2009) denote the estimator μ( ˆ θˆ ) as the robust EBLUP (REBLUP) of μ. Under the unit level model, the REBLUP corresponding to the EBLUP (3) for the area mean yˆ i is given by

ψ yˆ i

= Ni−1

⎧ ⎨ ⎩

yi j +



j∈si

ψ

(x iTj βˆ +

j∈ri

⎫ ⎬

ψ z iTj vˆi ) . ⎭

(8)

2.2 M-quantile-based small area methods An alternative approach to outlier robust small area estimation is the M-quantile regression-based method introduced by Chambers and Tzavidis (2006). The M-quantile regression is an extension of the quantile regression, which is described in Koenker and Hallock (2001) and Koenker (2005). For an application of quantile regression in the context of business data see Landajo et al. (2008) or Abberger (1997). Breckling and Chambers (1988) extended the approach from Koenker and Bassett (1978) to quantile regression based on general influence functions ψ. Chambers and Tzavidis (2006) assumed that a linear M-quantile regression model can be denoted as Q q (x, ψ) = x T β ψ (q). Hence, this model allows for different regression parameters β ψ depending on the M-quantile coefficient q. The regression vector β ψ can be computed using an iteratively reweighted least squares algorithm to solve nj 



ψq y j − x Tj β ψ (q) x j = 0.

(9)

j=1

Chambers and Tzavidis (2006) characterized the variability across the population by the M-quantile coefficients of the population units. That is, for each unit j in area i with values y j and x j , the M-quantile coefficients is determined such that Q q j (x j , ψ) = y j . They defined the area level M-quantile coefficients q i by suitably averaging the individual q j coefficients of each sample value, e.g. as qˆ i = mean{q j | j ∈ i}.

(10)

123

T. Schmid and R. T. Münnich

These qˆ i can be interpreted as pseudo-area effects similar to vi in the model (1). The M-quantile estimator of the small area mean, y i , is given as MQ yˆ i = Ni−1

⎧ ⎨ ⎩

j∈si

yi j +

 j∈ri

⎫ ⎬

(x iTj βˆ ψ (qˆ i )) , ⎭

(11)

where βˆ ψ (qˆ i )) denotes the estimated regression vector of the fitted M-quantile regression line at qˆ i . More information about M-quantile small area estimation can be found in Tzavidis et al. (2008, 2010). Salvati et al. (2012) defined a spatial extension to linear M-quantile regression which models spatial non-stationarity in the mean. 3 Spatial robust EBLUP The REBLUP estimator (8) by Sinha and Rao (2009) is robust to model misspecifications or outliers. Since it does not consider spatial dependencies in the data, a spatial extension of the REBLUP based on the simultaneously autoregressive (SAR) model will be introduced subsequently and denoted by SREBLUP. Further, a parametric bootstrap method for the MSE estimation will be presented. 3.1 Estimation and prediction In this part, we define a spatial extension of the robust EBLUP estimator given by (8). Therefore, we have to extend the general linear mixed model (1) to allow for spatial correlated area effects v. In general, there are two approaches to describe the spatial dependency in v: the conditional autoregressive (CAR) models and the SAR models. The CAR model regards only first-order dependencies in the study, whereas the SAR model additionally takes second-order dependencies into account. Due to the fact that the SAR model displays the practical situation, SAR models have been used most of the times especially in the frequentist framework. Similarly to Salvati (2004), Petrucci et al. (2005), Molina et al. (2009) and Salvati et al. (2012), we focus in this paper on the SAR model since the CAR model is more convenient in a Bayesian framework (cf. Banerjee et al. 2004). Detailed information on the CAR model can be found in Besag (1974), Sun et al. (1999), Vogt and Münnich (2009) or Vogt (2010). Therefore, the volatilities of the fixed part of the model X β are the results of a SAR model and these are given by v = (I − pW )−1 u,

(12)

where u is a vector of independent error terms with zero mean and variance σu2 . The parameter p denotes the spatial autoregressive parameter and W is the D × D proximity matrix (cf. Anselin 1992 or Cressie 1993). The matrix W describes the neighborhood structure between the small areas and p defines the strength of the

123

Spatial robust small area estimation

spatial dependencies among the random effects. (I − pW ) is invertible if the spatial  parameter p ∈ [min(λi )]−1 , [max(λi )]−1 , where λi are the eigenvalues of the matrix W . Beyond that, the matrix W is defined in a row standardized form, in which the row elements sum up to one. In this situation, W is not symmetric but it is row stochastic. The model with spatial correlated random effects is given by y = Xβ + Zv + e,

(13)

where v and e have covariance matrices

−1 G = σu2 (I − pW )(I − pW T )

(14)

R = diag(σe2 )

(15)

respectively and the covariance matrix of y is obtained by

−1 V = R + ZG Z T = diag(σe2 ) + Zσu2 (I − pW )(I − pW T ) ZT .

(16)

Similar to Sinha and Rao (2009), we maximize the density of y, but now with respect to β, θ and p leading to ML-equations for the estimation of β, θ and p: X T V −1 ( y − Xβ) = 0 

∂ V −1 ∂V =0 V ( y − Xβ)T − tr V −1 ( y − Xβ)V −1 ∂θl ∂θl 

∂ V −1 ∂V V ( y − Xβ)T − tr V −1 = 0. ( y − Xβ)V −1 ∂p ∂p The derivations

∂V ∂θl

and

∂V ∂p

(17) (18) (19)

are given by

∂V = Z L Z T with L = ((I − pW )(I − pW T ))−1 ∂σu2 ∂V =I ∂σe2

∂V = Zσu2 −L(2 pW W T − W − W T )L Z T . ∂p

(20) (21) (22)

For the estimation of the spatial outlier robust estimators β, θ and p, we robustify the above equations via 1

α(β) = X T V −1 U 2 ψ(r ) = 0

 1 ∂ V −1 1 ∂V V U 2 ψ(r ) − tr V −1 K =0 (θl ) = ψ T (r )U 2 V −1 ∂θl ∂θl 

1 ∂ V −1 1 ∂V V U 2 ψ(r ) − tr V −1 K = 0,

( p) = ψ T (r )U 2 V −1 ∂p ∂p

(23) (24) (25)

123

T. Schmid and R. T. Münnich

where r , U and K are defined according to Eqs. (4) and (5). We apply an iterative algorithm, e.g. the Newton-Raphson algorithm to solve the robust equations via

−1 α(β (m) ) β (m+1) = β (m) + α  (β (m) )

−1 θ (m+1) = θ (m) +  (θ (m) ) (θ (m) )

−1 p (m+1) = p (m) +  ( p (m) )

( p (m) ),

(26) (27) (28)

where  ( p (m) ) is the derivative of ( p) evaluated at p (m) . The complete algorithm for the calculation of the spatial robust estimates β, θ and p can be denoted as follows: 1. Choose initial values for β (0) , θ (0) and p (0) . 2. (a) Compute β (m+1) iteratively from Eq. (26). (b) Compute θ (m+1) iteratively from Eq. (27). (c) Compute p (m+1) iteratively from Eq. (28). (d) Set m = m + 1. 3. Continue step 2 until convergence is achieved. Denote the estimates at convergence ψ,sp ψ,sp to be the spatial robust ML-estimators βˆ , θˆ and pˆ of β, θ and p respectively, where the superscript ψ stands for influence function ψ (e.g. Huber function) and the superscript sp is an indicator of spatial influence to the estimators. ψ,sp ψ,sp Afterwards, we use βˆ , θˆ and pˆ for the estimation of the spatial robust area effects v by plugging in the spatial robust estimates in Fellner’s equation (cf. Fellner 1986) (29) Z T R−1/2 ψ(R−1/2 ( y − Xβ − Zv)) − G −1/2 ψ(G −1/2 v) = 0

and solve this equation with an iterative algorithm, e.g. the Newton-Raphson algorithm and set vˆ ψ,sp = v (m+1) as our spatial robust area effects estimator. The spatial robust ψ,sp ψ,sp ψ,sp estimates βˆ , θˆ , vˆ and pˆ are used to estimate (2) as μ( ˆ θˆ

ψ,sp

) = l T βˆ

ψ,sp

+ m T vˆ ψ,sp

(30)

ψ,sp

as the spatial robust for specified vectors l and m. We refer this estimator μˆ θˆ EBLUP (SREBLUP) of μ. Under the unit level model, the SREBLUP corresponding to the EBLUP (3) and REBLUP (8) for the area mean yˆ i is given by ψ,sp yˆ i

= Ni−1

⎧ ⎨ ⎩

j∈si

yi j +

 j∈ri

(x iTj βˆ

ψ,sp

+

⎫ ⎬

ψ,sp z iTj vˆi ) . ⎭

(31)

Further details regarding the SREBLUP estimator are available in Schmid (2011).

123

Spatial robust small area estimation

3.2 Mean squared error ψ,sp In the following, the quality of the SREBLUP estimator μ( ˆ θˆ ) from (30) is measured. Sinha and Rao (2008) noted that the MSE of the REBLUP estimator from (7) cannot be expressed in any closed analytic expression due to its complexity. Therefore, a resampling method is used for the MSE estimation of the spatial extension of the REBLUP estimator. According to Hall and Maiti (2006) and Sinha and Rao (2009), a parametric bootstrap technique is presented for the MSE estimation of the SREBLUP: ψ,sp ψ,sp , θˆ and pˆ sp generate the area-specific random effects v ∗ from 1. For a given βˆ ψ ψ N (0, G(θˆ , pˆ sp )) and e∗ from N (0, R(θˆ , pˆ sp )) to create a bootstrap sample via

y∗ = X βˆ

ψ,sp

+ Zv ∗ + e∗

(32)

and evaluate the corresponding bootstrap population parameter μ( ˆ βˆ

ψ,sp

, θˆ

ψ,sp

, pˆ sp , v ∗ ) = l T βˆ

ψ,sp

+ mT v∗

(33)

as defined by Eq. (30). 2. Generate B bootstrap samples ( y∗(1) , ..., y∗(b) , ..., y∗(B) ) from model (32) and compute the bootstrap population parameters via Eq. (30) ψ,sp ψ,sp

ψ,sp , θˆ , pˆ sp , v ∗(b) = l T βˆ + m T v ∗(b) . μˆ (b) βˆ

(34)

(b)ψ,sp (b)ψ,sp 3. Derive the spatial robust bootstrap estimates of βˆ , θˆ , pˆ (b)sp and (b)ψ,sp ∗(b) of the bth bootstrap sample y and compute the bootstrap estimate vˆ ψ,sp ˆ of the spatial robust EBLUP μ( ˆ θ ) from (30) by

(b)ψ,sp (b)ψ,sp

(b)ψ μˆ ∗(b) βˆ , θˆ , pˆ (b)sp , v (b)ψ,sp = l T βˆ + m T vˆ (b)ψ . 4. Obtain the bootstrap estimate of the MSE of the SREBLUP μ( ˆ θˆ by M S E(μˆ S R E B LU P ) =

ψ,sp

(35)

) from (30)

B 1  ∗(b) ˆ (b)ψ,sp ˆ (b)ψ,sp (b)sp (b)ψ,sp μˆ (β ,θ , pˆ ,v ) B b=1

2 ψ,sp ψ,sp − μˆ (b) (βˆ , θˆ , pˆ sp , v ∗(b) ) , (36)

where μˆ (b) is given in Eq. (34) and μˆ ∗(b) in (35). Further details regarding the MSE of the SREBLUP estimator are available in Schmid and Münnich (2011) and Schmid (2011).

123

T. Schmid and R. T. Münnich

4 Simulation study In order to compare the performance of the different estimators from the past sections, a model-based simulation was conducted. Distributions of point and MSE estimates for different estimators of interest were evaluated by drawing R = 1,000 samples from different theoretic distributions according to the underlying models. Further, different parameter setups were chosen as scenarios in order to allow for violations of the normal assumptions and to include spatial dependencies in the data. In a first step, the synthetic population for the universe was generated for m = 40 small areas of size Ni = 100. This universe includes the values of the covariates x i j and the geographical information (longitude and latitude) on individual level. The auxiliary variables were normally distributed with mean μx = 1 and variance σx2 = 1 xi j ∼ N (1, 1). Regarding the geographical information, it is assumed that all units inside the specific areas are centroid, i.e. they all have the same longitude and latitude within the respective area. These longitude and latitude information were generated independently from U [0, 1]. The sample of size n i = 5 was selected from the synthetic universe by simple random sampling without replacement within each area. The population and the samples were held fixed throughout the simulation study. According to Chambers et al. (2013), a spatial nested error regression model yi j = 100 + 4x i j + vi + ei j

(37)

is used as a core model for the generation of the dependent variable y, where the error term ei j and the random area effects vi follow specific assumptions depending on the respective scenario. In correspondence with the simulation study in Sinha and Rao (2009), the different scenarios were only generated under symmetric outlier contamination in the area and individual effects in the model (37). The focus in this model-based simulation study is put on symmetric outliers to keep the number of different scenarios, and hence the calculation time, on a moderate level. For more information regarding the calculation time of the estimators see Schmid (2011). The area effects vi and the error terms ei j were created using the model from vi ∼ (1 − γv ) · N (0, G) + γv · N (0, 20) ei j ∼ (1 − γe ) · N (0, 4) + γe · N (0, 150), where γi and γe control the amount of contamination in the data. Setting γe = γv = p = 0 yields the standard EBLUP model with no contamination and no spatial correlation. G is the SAR dispersion matrix defined by the equation from (14). The spatial autoregressive parameter p is set to 0.8 for the spatial scenarios and to 0 for the nonspatial settings. The focus is set here on positive correlation because this setting is more realistic in the case of practical situations. The neighborhood matrix W is given by the nearest neighbor approach. Each small area j is defined as a neighbour of area i, if the Euclidean distance between the two areas is less than 0.15. This approach

123

Spatial robust small area estimation

gives a neighborhood matrix containing the values 0 and 1. Afterwards, W is row standardized to fulfill the requirements of the SAR model (cf. Petrucci et al. 2005), leading to a non-symmetric but row stochastic neighborhood matrix. The value 0.15 is motivated because of the fact that every small area has on average 5 neighbors in this situation, which is a realistic assumption in practical situations. The W matrix was kept fixed for all settings. Within the study, γv = γe = 0.05 was chosen for the amount of contamination. Altogether four non-spatial and four spatial scenarios were chosen. The set {(0, 0), (0, e), (v, 0), (v, e)} indicates the four scenarios by (0, 0) no contamination, (0, e) contamination only in the individual error term, (v, 0) contamination only in the area term and (v, e) contamination in both components. The index p denotes the spatial scenarios {(0, 0) p , (0, e) p , (v, 0) p , (v, e) p }. Each scenario was independently simulated R = 1,000 times. The random effects vi were created according to specific settings and held fixed for the specific scenario, whereas the error term ei j is generated in each case for the R = 1,000 Monte-Carlo simulation runs independently according to the particular scenario. The fixed area effects reflect a finite population view indicating different interceps in the given areas which follow a normal distrution. First tests using further randomization of the random effects did not indicate a change of the essence and were omitted by reasons of severe computational effort. Five different estimators are used in the simulation study—the standard EBLUP defined in (3), which serves as the standard reference; the robust EBLUP of Sinha and Rao (2009) (8); the M-quantile estimator defined in (11) from Chambers and Tzavidis (2006); and the spatial extension of the robust EBLUP, the SREBLUP defined in (31). Furthermore, the generalized regression (GREG) estimator (cf. Särndal et al. 1992) is investigated in the simulation study. The spatial extension of the M-quantile estimator of Salvati et al. (2012) models spatial non-stationarity instead of spatial correlation between the area effects and hence it is not in the focus of our simulation study. The MSE estimation of the SREBLUP is implemented via the parametric bootstrap procedure from (36) and the MSE of the REBLUP is investigated with the bootstrap method proposed by Sinha and Rao (2009). According to Chambers et al. (2013), the bootstrap methods were conducted for R = 500 Monte-Carlo simulations with B = 100 bootstrap replicates. Next to the resampling methods, the MSE of the EBLUP and MQ are also evaluated with analytic approximations. The uncertainty of the EBLUP was measured with the Prasad-Rao estimator (cf. Prasad and Rao 1990) whereas the MSE of the MQ estimator is calculated as presented by Chambers and Tzavidis (2006). The MSE the GREG is estimated as proposed by Särndal et al. (1992). The following quality measures are used to evaluate the performance of the estimators:

123

T. Schmid and R. T. Münnich

– The relative root mean squared error (RRMSE) for each area i:  ⎛ ⎞2  A   R ˆ y − yi 1 ⎝ i,r ⎠ , R R M S E iA =  R yi r =1

where A represents any estimator proposed above and y i indicates the true value of the mean in area i. R denotes the number of Monte-Carlo replications. – The relative bias (RB) for each area i: R BiA

A R 1  yˆ i,r − y i = . R yi r =1

– The coverage rate for a 95 % confidence interval (CI): The coverage rate is the proportion of the R = 500 Monte-Carlo replications that the confidence intervals include the true value. The confidence intervals are defined by:  C Ii =

y iA

−z



M S E(y iA ), y iA

+z



M S E(y iA )

 ,

where z = 1.96 is the 97.5 % quantile of the normal distribution. Figure 1 shows the results of the RRMSE for the different non-spatial simulation scenarios and the different estimators. The plot indicates that the EBLUP-based methods as well as the MQ estimator outperform the GREG in all scenarios. It can be also observed that the robust EBLUP methods (REBLUP and SREBLUP) lead to better results in the scenarios with outlier contamination in the data especially in the settings (0, e) and (v, e). Table 1 reveals a more distinguished view of the RRMSE. As expected, the SREBLUP and the REBLUP reveals nearly identical results in the non-spatial scenarios. The GREG shows really good results in the scenarios (v, 0) and (v, e) in the areas with outlying random effects, but performs badly in the areas without outlier random effects. The reason behind is that the GREG uses as a modelassisted estimator the information of the slope rather than the one of the intercept in the areas and hence is less sensitive against area outliers. Regarding the RB in the non-spatial settings, Fig. 2 indicates that the EBLUP, the REBLUP, the SREBLUP and the MQ have no bias in all scenarios. Similar to the RRMSE, the REBLUP and the SREBLUP perform on the same level in term of RB in the non-spatial settings. The behavior of the MSE estimators in the non-spatial scenarios is summarized in Table 2. The table shows the mean value of coverage rates in [%] as well as the CI length divided by the true value for each scenario. It can be observed that the GREG provides well coverage rates at an expense of large CI. In contrast, the resampling methods for the REBLUP and SREBLUP provide smaller CI, but suffer from a little undercoverage. The Prasad-Rao method for the EBLUP provides as expected good coverage rates in the scenario without outliers. In the settings with individual outliers, the EBLUP reveals slightly overcoverage and large CI. That

123

Spatial robust small area estimation 0

1

2

(v,0)

(v,e)

(0,0)

(0,e)

3

4

MQ SREBLUP REBLUP EBLUP GREG

MQ SREBLUP REBLUP EBLUP GREG 0

1

2

3

4

Fig. 1 Boxplots of the RRMSE [%] in the non-spatial scenarios with 5 % outlier contamination

Table 1 RRMSE in % in the non-spatial scenarios with 5 % outlier contamination Mean values of the RRMSE No outliers

Indiv. outliers

Area outliers

Both

Scenarios

[0,0]

[v,0]

[0,e]

[v,e]

[v,0]

[v,e]

Areas

1–40

1–38

1–40

1–38

39–40

39–40

GREG

0.931

0.926

1.489

1.486

0.897

1.386

EBLUP

0.646

0.684

0.812

0.870

1.456

2.531

REBLUP

0.661

0.675

0.695

0.707

1.205

1.396

SREBLUP

0.644

0.674

0.677

0.706

1.194

1.376

MQ

0.694

0.698

0.747

0.747

1.385

1.595

indicates that individual outliers lead to higher variance estimates for the EBLUP. As expected, the parametric bootstrap for the REBLUP and SREBLUP perform similarly in the non-spatial scenarios. More sophisticated methods like block bootstrapping (cf. Chambers and Chandra 2012) or the evaluation of bootstrap confidence intervals (cf. Chatterjee et al. 2008) can yield slightly better coverage rates for the REBLUP and SREBLUP, however this was not the focus of this paper. Figure 3 shows the plots of the RRMSE in the spatial scenarios with 5 % outlier contamination. It can be observed that the EBLUP-based small area methods have lower RRMSE than the GREG especially in the scenarios (0, e) and (v, e). We can

123

T. Schmid and R. T. Münnich −4

−2

0

(v,0)

(v,e)

(0,0)

(0,e)

2

4

MQ SREBLUP REBLUP EBLUP GREG

MQ SREBLUP REBLUP EBLUP GREG −4

−2

0

2

4

Fig. 2 Boxplots of the RB [%] in the non-spatial scenarios with 5 % outlier contamination Table 2 Coverage rates and relative CI length (in brackets) in the non-spatial scenarios with 5 % symmetric outlier contamination Mean value of coverage rates in % Scenarios

[0,0]

[0,e]

[v,0]

[v,e]

GREG

95.5 (4.54)

95.9 (5.00)

95.6 (4.56)

96.6 (5.05)

EBLUP

94.9 (2.04)

97.3 (3.11)

91.1 (2.19)

95.3 (3.30)

REBLUP Boot

93.0 (1.92)

92.9 (1.98)

92.3 (2.10)

92.7 (2.16)

SREBLUP Boot

93.7 (1.86)

93.7 (1.97)

91.1 (2.05)

93.0 (2.12)

MQ

84.8 (2.04)

86.2 (2.33)

88.8 (2.44)

89.2 (2.36)

Estimators

see that that the individual contamination has much more influence on the quality of the estimates than the area-specific contamination. The Figure also shows, as expected, that the robust methods outperform the non-robust estimators in the scenarios with outliers. Table 3 reveals a positive effect of spatial modeling. The SREBLUP has in total a lower RRMSE compared to the EBLUP and the REBLUP in all scenarios. If we compare the spatial and the non-spatial settings, we see that the spatial information leads to a lower RRMSE of the SREBLUP compared to REBLUP. In the areas with outlying random effects in the settings (v, 0) and (v, e), the GREG reveals the lowest RRMSE compared to the other methods as in the non-spatial scenarios. The performance of the estimators regarding the RB in the spatial scenarios is shown in Fig. 4, which makes clear that all

123

Spatial robust small area estimation 0

1

2

(v,0)p

(v,e)p

(0,0)p

(0,e)p

3

4

MQ SREBLUP REBLUP EBLUP GREG

MQ SREBLUP REBLUP EBLUP GREG 0

1

2

3

4

Fig. 3 Boxplots of the RRMSE [%] in the spatial scenarios with 5 % outlier contamination Table 3 RRMSE in % in the spatial scenarios with 5 % outlier contamination Mean values of the RRMSE No outliers Scenarios

[0,0]

Indiv. outliers

Area outliers

Both

[v,0]

[0,e]

[v,e]

[v,0]

[v,e]

Areas

1–40

1–38

1–40

1–38

39–40

39–40

GREG

0.921

0.921

1.459

1.460

0.938

1.541

EBLUP

0.876

0.889

1.068

1.144

1.879

2.929

REBLUP

0.889

0.870

0.924

0.924

1.811

1.709

SREBLUP

0.819

0.836

0.848

0.885

2.029

1.927

MQ

0.903

0.890

0.953

0.948

2.032

2.115

estimators are unbiased in these settings. The results of the MSE estimation methods in the spatial scenario are shown in Table 4. In accordance with the non-spatial setting, the GREG exposes good coverage rates at the price of a larger CI. The investigated bootstrap methods perform similar. The spatial patterns in the data did not result in higher coverage rates or smaller confidence intervals in the parametric bootstrap of the SREBLUP in some areas. The presented parametric bootstrap methods for the REBLUP and SREBLUP suffer from undercoverage, whereas the average length of the CI is quite efficient just as in the non-spatial scenario. The MSE estimation method of the SREBLUP only indicates coverage rates of approximately 92 % in all spatial settings, similar to the bootstrap of the REBLUP. However, future research will have

123

T. Schmid and R. T. Münnich −4

−2

0

(v,0)p

(v,e)p

(0,0)p

(0,e)p

2

4

MQ SREBLUP REBLUP EBLUP GREG

MQ SREBLUP REBLUP EBLUP GREG −4

−2

0

2

4

Fig. 4 Boxplots of the RB [%] in the spatial scenarios with 5 % outlier contamination Table 4 Coverage rates and relative CI length (in brackets) in the spatial scenarios with 5 % symmetric outlier contamination Mean value of coverage rates in % Scenarios

[0,0]

[0,e]

[v,0]

[v,e]

95.7 (4.56)

96.2 (5.02)

96.0 (4.56)

96.3 (5.09)

EBLUP

92.0 (2.13)

95.5 (3.13)

89.5 (2.29)

92.4 (3.35)

REBLUP Boot

91.9 (2.09)

92.7 (2.22)

92.1 (2.22)

91.3 (2.38)

SREBLUP Boot

93.6 (1.94)

93.7 (2.08)

93.0 (2.09)

91.5 (2.31)

MQ

94.3 (2.73)

92.7 (2.62)

94.4 (2.70)

92.0 (2.65)

Estimators GREG

to show whether these results are valid in more general settings and in the case of design-based studies with sophisticated sampling designs. 5 Summary and outlook In this paper we propose a spatial extension of the robust EBLUP approach from Sinha and Rao (2009). This estimator combines outlier robust small area estimation with the modeling of spatial structures in a single model. The results of the simulation study suggest that the proposed spatial robust EBLUP has the potential to lead to considerably better results compared to other estimators. In terms of the measure of

123

Spatial robust small area estimation

uncertainty, the proposed parametric bootstrap of the SREBLUP performs on the same level than the resampling method for the REBLUP. First outcomes of design-based simulations on real business data show that the SREBLUP incorporates the spatial patterns in the data and leads to more precise results. However, the design-based simulation results are just a first picture and a larger number of simulations in a design-based environment is necessary to generalize the results. Next to that, improvements of the numerical algorithms should help to find proper applications in official statistics in the context of business data. Acknowledgements The research was conducted within the BLUE-ETS research project which is funded by the European Commission within the 7th Framework Programme. For more information on the project, we refer to the project page http://www.blue-ets.eu. The first author was supported by the Foundation of German Economy. The authors are grateful to the SAMPLE project (http://www.sample-project.eu/) and Chambers et al. (2013) for providing R-code used in the simulations in this paper. The authors would like to thank the Editor-in-Chief and two anonymous referees for their very valuable comments which helped to improve the paper.

References Abberger K (1997) Quantile smoothing in financial time series. Stat Pap 38(2):125–148 Anselin L (1992) Spatial econometrics: methods and models. Kluwer, London Banerjee S, Carlin B, Gelfand AE (2004) Hierarchical modeling and analysis for spatial data. Chapman & Hall, Boca Raton Battese GE, Harter RM, Fuller WA (1988) An error component model for prediction of county crop areas using survey and satellite data. J Am Stat Assoc 83(401):28–36 Besag J (1974) Spatial interaction and the statistical analysis of lattice systems (with discussion). J Royal Stat Soc Ser B 36:192–236 Breckling J, Chambers R (1988) M-quantiles. Biometrika 75(4):761–771 Chambers R (1986) Outlier robust finite population estimation. J Am Stat Assoc 81:1063–1069 Chambers R, Chandra H (2012) A random effect block bootstrap for clustered data. J Comput Graph Stat. doi:10.1080/10618600.2012.681216 Chambers R, Chandra H, Salvati N, Tzavidis N (2013) Outlier robust small area estimation. J Royal Stat Soc Ser B Stat Methodol. doi:10.1111/rssb.12019 Chambers R, Tzavidis N (2006) M-quantile models for small area estimation. Biometrika 93(2):255–268 Chatterjee S, Lahiri P, Li H (2008) Parametric bootstrap approximation to the distribution of eblup and related predition intervals in linear mixed models. Ann Stat 36(3):1221–1245 Cressie N (1993) Statistics for spatial data. Wiley, New York Dasiou D, Moyssiadis C (2001) The 50 % breakdown point in simultaneous m-estimation of location and scale. Stat Pap 42(2):243–252 Fellner WH (1986) Robust estimation of variance components. Technometrics 28(1):51–60 Hall P, Maiti T (2006) On parametric bootstrap methods for small area prediction. J Royal Stat Soc Ser B 68(2):221–238 Hampel FR, Ronchetti EM, Rousseeuw PJ, Stahel WA (1986) Robust statistics: the approach based on influence functions. Wiley, New York Harville DA (1977) Maximum likelihood approaches to variance component estimation and to related problems. J Am Stat Assoc 72(358):320–340 Jiang J, Lahiri P (2006) Mixed model prediction and small area estimation. Test 15(1):1–96 Koenker R (2005) Quantile regression. Cambridge University Press, New York Koenker R, Bassett G (1978) Regression quantiles. Econometrica 46:33–50 Koenker R, Hallock KF (2001) Quantile regression. J Econ Perspect 15(4):143–156 Landajo M, de Andrés J, Lorca P (2008) Measuring firm performance by using linear and non-parametric quantile regressions. J Royal Stat Soc Ser C 57(2):227–250 Molina I, Salvati N, Pratesi M (2009) Bootstrap for estimating the mean squared error of the spatial eblup. Comput Stat 24(3):441–458

123

T. Schmid and R. T. Münnich Münnich R, Gabler S, Ganninger M (2007) Some remarks on the register-based census 2010/2011 in Germany. In: Proceedings of the workshop innovative methodologies for censuses in the new millennium, Southampton Niemiro W, Wesolowski J (2010) Synthetic and composite estimation under a superpopulation model. Stat Pap 51(3):497–509 Petrucci A, Pratesi M, Salvati N (2005) Geographic information in small area estimation: small area models and spatially correlated random area effects. Stat Transit 7(3):609–623 Prasad NGN, Rao JNK (1990) The estimation of the mean squared error of small area estimators. J Am Stat Assoc 85(409):163–171 Pratesi M, Salvati N (2008) Small area estimation: the eblup estimator based on spatially correlated random area effects. Stat Methods & Appl 17:113–141 Rao JNK (2003) Small area estimation. Wiley, New York Richardson AM (1997) Bounded influence estimation in the mixed linear model. J Am Stat Assoc 92(437):154–161 Richardson AM, Welsh AH (1994) Asymptotic properties of restricted maximum likelihood (reml) estimates for hierarchical mixed linear models. Aust J Stat 36:31–43 Richardson AM, Welsh AH (1995) Robust restricted maximum likelihood in mixed linear models. Biometrics 51(4):1429–1439 Salvati N (2004) Small area estimation by spatial models: the spatial empirical best linear unbiased prediction (spatial eblup). Technical report. University of Florence, Department of Statistics Salvati N, Tzavidis N, Pratesi M, Chambers R (2012) Small area estimation via M-quantile geographically weighted regression. Test 21(1):1–28 Schmid T (2011) Spatial robust small area estimation applied on business data. Ph. D. thesis, University of Trier Schmid T, Münnich R (2011) Parametric bootstrap mean squared error estimation for the spatial robust eblup. In: Proceedings of the second ITACOSM conference, survey research methods and applications, pp 211–214 Singh B, Shukla G, Kundu D (2005) Spatio-temporal models in small area estimation. Surv Methodol 31(2):183–195 Sinha SK, Rao JNK (2008) Robust small area estimation under unit level models. In: Proceedings of the survey research methods section, American statistical association, vol 82 Sinha SK, Rao JNK (2009) Robust small area estimation. Can J Stat 37(3):381–399 Särndal C-E, Swensson B, Wretman J (1992) Model assisted survey sampling. Springer, New York Sun D, Tsutakawa RK, Speckman P (1999) Posterior distribution of hierarchical models using car(1) distributions. Technical report. National Institute of Statistical Sciences Tzavidis N, Marchetti S, Chambers R (2010) Robust estimation of small area means and quantiles. Aust N Z J Stat 52(2):167–186 Tzavidis N, Salvati N, Pratesi M, Chambers R (2008) M-quantile models with application to poverty mapping. Stat Methods Appl 17(3):393–411 Vogt M (2010) Bayesian spatial modeling. Ph. D. thesis, University of Trier Vogt M, Münnich R (2009) On the existence of a posterior distribution for spatial mixed models with binomial responses. Metron 67(2):199–207 Welsh A, Ronchetti E (1998) Bias-calibrated estimation from sample surveys containing outliers. J Royal Stat Soc Ser B 60:413–428 Zadlo T (2009) On mse of eblup. Stat Pap 50(1):101–118

123