COVARIATE SELECTION FOR SMALL AREA ESTIMATION IN REPEATED SAMPLE SURVEYS

van den Brakel, Jan A.; Buelens, Bart Article COVARIATE SELECTION FOR SMALL AREA ESTIMATION IN REPEATED SAMPLE SURVEYS Statistics in Transition New Series Provided in Cooperation with: Polish Statistical Association Suggested Citation: van den Brakel, Jan A.; Buelens, Bart (2015) : COVARIATE SELECTION FOR SMALL AREA ESTIMATION IN REPEATED SAMPLE SURVEYS, Statistics in Transition New Series, ISSN 2450-0291, Exeley, New York, NY, Vol. 16, Iss. 4, pp. 523-540, https://doi.org/10.21307/stattrans-2015-031 This Version is available at: https://hdl.handle.net/10419/207789 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/ STATISTICS IN TRANSITION new series and SURVEY METHODOLOGY 523 STATISTICS IN TRANSITION new series and SURVEY METHODOLOGY Joint Issue: Small Area Estimation 2014 Vol. 16, No. 4, pp. 523–540 COVARIATE SELECTION FOR SMALL AREA ESTIMATION IN REPEATED SAMPLE SURVEYS Jan A. van den Brakel1,Bart Buelens2 ABSTRACT If the implementation of small area estimation methods to multiple editions of a repeated sample survey is considered, then the question arises which covariates to use in the models. Applying standard model selection procedures independently to the different editions of the survey may identify different sets of covariates for each edition. If the small area predictions are sensitive to the different models, this is undesirable in official statistics since monitoring change over time of statistical quantities is of utmost importance. Therefore, potential confounding of true change and methodological alterations should be avoided. An approach to model selection is proposed resulting in a single set of covariates for multiple survey editions. This is achieved through conducting covariate selection simultaneously for all editions, minimizing the average of the edition-specific conditional Akaike Information Criteria. Consecutive editions of the Dutch crime victimization survey are used as a case study. Municipal estimates of three survey variables are obtained using area level models. The proposed averaging strategy is compared to the standard method of considering each edition separately, and to an elementary approach using covariates selected in the first edition. Resulting models, point estimates and MSE estimates are analyzed, indicating no substantial adverse effects of the conceptually attractive averaging strategy. Key words: area level models, cAIC, Hierachical Bayesian predictors. 1. Introduction At national statistical institutes, estimation procedures for surveys based on probability samples are traditionally based on design-based or model-assisted inference procedures. Well-known examples are the π-estimator (Narain, 1951; Horvitz and Thompson, 1952) and the generalized regression estimator (Särndal, Swensson and Wretman, 1992). These approaches are particularly appropriate in the case of large sample sizes. In the case of small sample sizes, however, design-based and modelassisted estimators have unacceptably large variances. This occurs when estimates 1Statistics Netherlands, Department of Statistical Methods and Maastricht University, Department of Quantitative Economics. E-mail: [email protected] 2Statistics Netherlands, Department of Statistical Methods. E-mail: bb[email protected] 524 Van den Brakel and Buelens: Covariate selection for small area estimation are required for detailed breakdowns of the population in subpopulations or domains according to various socio-demographic or geographic classification variables. In such cases, model-based estimation procedures are required to increase the effective sample size of the separate domains with sample information observed in other domains or preceding periods. This class of estimation procedures is known in the literature as small area estimation (SAE) (Rao, 2003; Pfeffermann, 2013) and offers promising opportunities for official statistics (Boonstra et al., 2008). A common approach to introducing SAE in an existing survey is to apply SAE methods to historic editions of the survey, producing small area estimates for multiple past editions at the same time. This article focuses on the selection of covariates to be used in the SAE models in this setting. In the literature, model selection procedures mostly focus on the selection of optimal models for one particular survey data set (Claeskens and Hjort, 2008). If in each edition of a repeated survey a separate and different model is selected, the question arises to what extent the small area predictions are comparable over time. In official statistics potential confounding of estimates of change over time of some statistic with variations in the inference procedures must be avoided. This article contributes to the existing literature by addressing the question how to select a single optimal model for the production of SAE predictions for independent, repeated editions of a sample survey. An approach is proposed in which the model selection criterion is averaged over all available editions, leading to a single set of covariates to be used in each edition. This novel approach is compared to the standard approach of selecting a set of covariates for each edition independently using four past editions of the Dutch crime victimization survey. In addition, a simple scenario is included whereby covariates are selected using only the first of a series of survey editions. In this paper models are considered that only use cross-sectional correlation. Alternative approaches that combine cross-sectional and temporal data are proposed by Rao and Yu (1994), Datta et al. (1999) and Pfeffermann and Tiller (2006). These approaches might also be considered to select one single optimal model for subsequent survey editions. These approaches are not considered for implementation in the Dutch crime victimization survey since they are considerably more complex and computationally intensive. The article continues in Section 2 with a presentation of the SAE methods used and details covariate selection procedures. Section 3 introduces the crime victimization survey and potential covariates. Results are presented and discussed in Section 4. Conclusions are drawn in Section 5. STATISTICS IN TRANSITION new series and SURVEY METHODOLOGY 525 2. Methods 2.1. Small Area Estimation In small area estimation multilevel models are used to improve the estimation of small domain parameters. These models use relevant auxiliary information as covariates. In this article the area level model is used (Fay and Herriot, 1979), where the input data for the model are the direct estimates for the domains. Approaches to covariate selection discussed below can be applied to unit level models (Battese, Harter and Fuller, 1988) as well. The area level model is considered, since it takes the complexity of the sample design into account as the dependent variables of the model are the design-based estimates derived from the probability sample and available auxiliary information used in the weighting model of the generalized regression (GREG) estimator. Let ˆ θidenote the GREG estimates of the target variables θifor the domains i=1,...,m. In the area level model, the direct domain estimates are modeled with a measurement error model, i.e. ˆ θi=θi+ei, where ei denotes the sampling error with design variance ψi. The unknown domain parameter is modeled with available covariates for the i−th domain, i.e. θi=z0 iβ+vi, with ziaK-vector with the covariates zi,kfor domain i,βthe corresponding K-vector with fixed effects and vithe random area effects with variance σ2 v. For each variable a separate univariate model is assumed. Combining both components gives rise to the basic area level model, originally proposed by Fay and Herriot (1979): ˆ θi=z0 iβ+vi+ei,(1) with model assumptions vi iid ∼N(0,σ2 v)and ei ind ∼N(0,ψi).(2) It is assumed that viand eiare independent and that ψiis known. Model(1) is a linear mixed model and estimation often proceeds using Empirical Best Linear Unbiased Prediction (EBLUP), where the between domain variance σ2 vis estimated with the Fay-Herriot moment estimator, maximum likelihood or restricted maximum likelihood, see Rao (2003), ch. 6 for details. A weakness of these methods is that in some situations the estimated model variance tends to zero, see e.g. Bell (1999) and Rao (2003). To avoid these problems, the Hierarchical Bayesian (HB) approach is followed in this article, Rao (2003), section 10.3. Therefore, the basic area level model is expressed as an HB model by (1) and (2) and a flat prior on βand σ2 v. The HB estimates for θiand its MSE are obtained as the posterior mean and variance of θi. To account for the uncertainty in the between 526 Van den Brakel and Buelens: Covariate selection for small area estimation domain variance, integration over the posterior density for σ2 vis conducted. Estimates for the design variances ψiare available from the GREG estimator but are used as if the true design variances are known, which is a standard assumption in small area estimation. Therefore, it is important to provide reliable estimates for ψi. The stability of the estimates for ψiis improved using the following ANOVA-type pooled variance estimator ψi=1−fi ni S2 p, S2 p=1 n−m m ∑ i=1 (ni−1)S2 i;GREG, with fithe sample fraction in domain i,nithe sample size in domain i,n=∑m i=1ni and S2 i;GREG the estimated population variance of the GREG residuals. 2.2. Conditional AIC The model selection procedures discussed here are optimization routines, minimizing the conditional Akaike Information Criterion (cAIC) proposed by Vaida and Blanchard (2005). The cAIC is applicable to mixed models where the focus is on prediction at the level of clusters or areas (Vaida and Blanchard, 2005). It is defined as cAIC = −2L+2p, where Lis the conditional log-likelihood and pa penalty based on a measure for the model complexity. In the case of a fixed effects model, pis the number of model parameters. The random part of a mixed model also contributes to the number of model degrees of freedom pwith a value between 0 in the case of no domain effects (i.e. ˆ σ2 v=0) and the total number of domains min the case of fixed domain effects (i.e. ˆ σ2 v→∞). In the expression of the cAIC, pis the effective degree of freedom of the mixed model and is defined as the trace of the hat matrix H, which maps the observed data to the fitted values, i.e. ˆy=Hy, see Hodges and Sargent (2001). When comparing models, the one with the lowest cAIC value is preferred. 2.3. Covariate selection procedures Covariate selection procedures are aimed at establishing a set of covariates – in the present setting the fixed effects – to use in models specified by equation (1). This boils down to finding an optimal subset from a larger set of available candidates. All three methods detailed below proceed along the same lines: they follow a step-forward covariate selection strategy which starts from an intercept-only model STATISTICS IN TRANSITION new series and SURVEY METHODOLOGY 527 adding covariates one-by-one until there is no improvement in terms of the selection criterion. This may result in sub-optimal models as the procedure converges to a local minimum of the selection criterion but not necessarily to the global minimum (Claeskens and Hjort, 2008). The focus here, however, is on establishing a single set of covariates for use in repeated survey editions. Alternative search routines converging to the global minimum of the selection criterion can be applied analogously to the step-forward routine used here. Some general notation is introduced. When Ccandidate covariates are available for inclusion as a fixed effect in a model specified by equation (1), the set of selected covariates is denoted by sand the set of remaining covariates by r. For ease of use the candidate covariates are assumed to be ordered in a fixed but arbitrary order, so that they can be referred to by their index. For example, a model containing the jth covariate – with 1 ≤j≤C– as a fixed effect, can be identified by s={j}. Consequently, in such case r={i}i6=j. Evidently, the equality s∪r={1,...,C} always holds. Sets of selected and remaining covariates that are specific to a survey edition tare denoted by stand rtrespectively. 2.3.1 Selecting an optimal set for each edition separately For a series of independent cross sectional surveys repeated at times t=1,...,T, a standard covariate selection routine consists in selecting covariates for each edition indepentently. Covariate selection procedure ’stnd’. Repeat for all editions t∈ {1,...,T}: Initialization Set rt={1,...,C}and st={}, obtain the corresponding cAIC value and call this cAIC0. Set i=0. Repeat Attempt extending the model with one covariate: a / Set i=i+1. b / Calculate cAIC for all models st∪{ j},∀j∈rt, and call these cAICj. c / If min(cAICj)<cAICi−1then set cAICi=min(cAICj), extend stto include the corresponding covariate j, remove that covariate from rt. Until The model is not extended or all candidate covariates are included in the model. The result are sets stof selected covariates for each edition t. In general, stand st0 can be different for t6=t0. 528 Van den Brakel and Buelens: Covariate selection for small area estimation The generic model specification given by equation (1) is adapted to reflect the repeated nature of the survey. ˆ θi,t=z[stnd]0 i,tβt+vi,t+ei,t,(3) for i=1,...,mand t=1,...,T, with model assumptions vi,t iid ∼N(0,σ2 v,t)and ei,t ind ∼N(0,ψi,t).(4) The vectors z[stnd] i,tconsist of covariates contained in stat the level of the domains i, with stestablished through the stnd covariate selection procedure. 2.3.2 Selecting one optimal set for all editions simultaneously Since the standard method may result in different sets of covariates for different survey editions, an alternative is proposed here, resulting in a single set of covariates for all editions. Formally, the following procedure enforces that st=st0for all t,t0∈ {1,...,T}. Covariate selection procedure ’avrg’. Consider all survey editions t=1. . . Tsimultaneously. Initialization Let r={1,...,C}and s={}. Use rand sfor all t, obtain the corresponding cAIC values, and call these cAIC0,t. Define cAIC0= 1 T∑tcAIC0,t. Set i=0. Repeat Attempt extending the model with one covariate: a / Set i=i+1. b / For all editions t∈ {1,...,T}, calculate cAIC for all models s∪{j}, ∀j∈r, and call these cAICj,t. c / Define cAICj=1 T∑tcAICj,t. d / If min(cAICj)<cAICi−1then set cAICi=min(cAICj), extend sto include the corresponding covariate j, remove that covariate from r. Until The model is not extended or all candidate covariates are included in the model. This strategy is based on averaging the model selection criterion cAIC and results in a single set sof covariates to be used in all editions t. The corresponding model specification, with a fixed set of covariates for repeated surveys, is written as (3) STATISTICS IN TRANSITION new series and SURVEY METHODOLOGY 529 and (4) where the vectors z[stnd] i,tare replaced by vectors, say z[avrg] i,t, that consist of covariates contained in sat the level of the domains iat the time periods t, with s established through the avrg covariate selection procedure. 2.3.3 Selecting an optimal set based on the first edition only An elementary approach also resulting in a single set of covariates is to use the first edition of a series of repeated surveys to establish the set of covariates and to use these in all subsequent editions. Covariate selection procedure ’frst’. Apply procedure stnd for t=1 to obtain s1. The set of covariates s1obtained based on the first edition is used at all times. The model takes the form of (3) and (4) where the vectors z[stnd] i,tare replaced by vectors, say z[f rst] i,t, that consist of covariates contained in s1at the level of the domains iat the time periods t. This strategy is included to assess and illustrate its performance. In other settings than the one discussed in the present article, statisticians may be in a situation where a survey is foreseen to be repeated in the future, but SAE estimates are required at the time of the first edition. The only option then is to use that edition for covariate selection. 3. Data 3.1. Crime victimization survey The Dutch crime victimization survey underwent several redesigns in the past, including in 2008 and 2012. In the period from 2008 through 2011 the survey is known as the Integrated Safety Monitor (ISM). These four editions of the ISM are used as a case study in the present article. The purpose of the ISM is to publish information on crime victimization, public safety and satisfaction with police performance, among others. Each annual ISM sample is obtained independently through stratified simple random sampling of persons aged 15 years or older residing in the Netherlands. The population register serves as the sampling frame. The country is divided into 25 police districts, which are used as the stratification variable in the sample design. The yearly sample size of about 19,000 respondents is divided equally over the strata. In addition to this national sample, local authorities such as municipalities and police districts can draw supplementary samples in their own regions on a voluntary basis, with the purpose to obtain precise local estimates. 530 Van den Brakel and Buelens: Covariate selection for small area estimation These supplementary samples are also based on stratified simple random sampling, but now with a more detailed geographical stratification variable, usually neighborhood. Table 1 gives an overview of the oversampling and the number of respondents for the years 2008 through 2011. Participation in the oversampling scheme by local authorities was encouraged in the years 2009 and 2011 resulting in much larger samples in these editions. Table 1: Overview of response and oversampling in ISM surveys 2008 - 2011. 2008 2009 2010 2011 Number of oversampled municipalities 77 239 21 225 Size response national sample 16,964 19,202 19,238 20,325 Size response supplemental sample 45,839 182,012 19,982 203,621 Percentage of population in oversampled areas 29% 65% 16% 66% Data collection is based on a sequential mixed mode design using internet (WI), paper (PAPI), telephone interviewing (CATI) or face-to-face interviewing (CAPI). For the data collection of the additional regional samples the WI, PAPI and CATI modes are mandatory. The use of the CAPI mode is recommended but not mandatory since this mode is very costly. Statistical inference for official publication purposes is based on the GREG estimator. The inclusion probabilities in the ISM are determined by the sampling design, accounting for stratification and oversampling at regional levels. The GREG estimator uses a complex weighting scheme that is based on the auxiliary variables age, gender, ethnicity, urbanization, household size, police district, and the strata used in the regional oversampling scheme. In addition, the weighting scheme contains a component that calibrates the response to a fixed distribution over the data collection modes with the purpose to stabilize the measurement error between the subsequent editions of the ISM, (Buelens and Van den Brakel, 2015). Variance estimates are obtained with the standard Taylor series approximation of the GREG estimator, see Särndal, Swensson and Wretman (1992), ch. 6. The GREG estimator can be used to produce reliable official statistics for regions with relatively large sample sizes. With the aforementioned sample design this implies that the GREG estimator can be used to produce official statistics at the level of police districts and in the regions where additional samples are drawn also at the level of municipalities. For regions where no additional samples are drawn, sample sizes are too small to produce reliable estimates at the level of municipalities with the GREG estimator. Since there is a growing demand for such figures, SAE procedures are developed to produce reliable official statistics on crime victimization at the municipal level. Three important ISM variables under study are listed in STATISTICS IN TRANSITION new series and SURVEY METHODOLOGY 537 0008−Haarlemmerliede en Spaarnwoude 0144−Rhenen 0281−Midden−Drenthe 0418−Amsterdam 0.0 0.1 0.2 0.3 0.4 0.0 0.1 0.2 0.3 0.4 2008 2009 2010 20112008 2009 2010 2011 year estimate type GREG AVRG STND FRST 0008−Haarlemmerliede en Spaarnwoude 0144−Rhenen 0281−Midden−Drenthe 0418−Amsterdam 2 3 4 2 3 4 2008 2009 2010 20112008 2009 2010 2011 year estimate type GREG AVRG STND FRST 0008−Haarlemmerliede en Spaarnwoude 0144−Rhenen 0281−Midden−Drenthe 0418−Amsterdam 0.0 0.2 0.4 0.6 0.0 0.2 0.4 0.6 2008 2009 2010 20112008 2009 2010 2011 year estimate type GREG AVRG STND FRST Figure 1: Time series of GREG and SAE estimates obtained through the stnd,avrg and frst approaches for four municipalities for victim (top), degen (middle) and contpol (bottom). 538 Van den Brakel and Buelens: Covariate selection for small area estimation Appendix A: Auxiliary variables defined for municipalities westimmi: share of western immigrants in the population nonwestimmi: share of non-western immigrants in the population prov: province density: housing density (number of dwellings per square kilometer) logdens: natural logarithm of density sqrtdens: square root of density meanvalue: mean house value (available from housing register) carsphh: average number of cars owned by households young: share of population aged 15-30 old: share of population aged 65+ rent: share of houses that are rented (as opposed to owned) lowincome: share of households with a low income (nationwide in lowest quintile) highincome: share of households with a high income (nationwide in highest quintile) unemployed: share of population registered at the employment agency as looking for work totcrime: number of crimes registered by the Police per 1.000 inhabitants propcrimedef1: number of property crimes registered by the Police per 1.000 inhabitants (definition CBS) propcrimedef2: number of property crimes registered by the Police per 1.000 inhabitants (definition Bureau Veiligheid) biketheft: number of bicycle thefts registered by the Police per 1.000 inhabitants violcrime: number of violent crimes registered by the Police per 1.000 inhabitants mode2: share of non-interviewer administered modes (paper and web) in the ISM survey oversampled: binary variable indicating whether the municipality took part in the ISM oversampling scheme REFERENCES BATTESE, G.E., HARTER, R.M., FULLER, W.A., (1988). An error components model for prediction of county crop areas using survey and satellite data. Journal of the American Statistical Association, 83, 28–36. BELL, W.R., (1999). Accounting for uncertainty about variances in small area estimation. Technical report. Bulletin of the International Statistical Institute. BOONSTRA, H.J., (2012). hbsae: Hierarchical Bayesian Small Area Estimation, Manual R package version 1.0.. Statistics Netherlands, Heerlen. STATISTICS IN TRANSITION new series and SURVEY METHODOLOGY 539 BOONSTRA, H.J., VAN DEN BRAKEL, J.A., BUELENS, B., KRIEG, S., SMEETS, M., (2008). Towards small area estimation at Statistics Netherlands. METRON International Journal of Statistics, LXVI, 21–49. BUELENS, B., VAN DEN BRAKEL, J.A., (2014). Model selection for small area estimation in repeated surveys. Discussion paper 201423, Statistics Netherlands, Heerlen. http://www.cbs.nl/NR/rdonlyres/308ED398- 714A-41A4-A57 C-9DCCC3F30D35/0/201423x10pub.pdf BUELENS, B., VAN DEN BRAKEL, J.A., (2015). Measurement error calibration in mixed mode surveys. Sociological Methods and Research, 44, 391– 426. CLAESKENS, G., HJORT, N.L., (2008). Model selection and model averaging, Cambridge series on statistical and probabilistic mathematics, Cambridge University Press. DATTA, G.S., LAHIRI, P., MAITI, T., LU, K.L., (1999). Hierarchical Bayes estimation of unemployement rates for the states of the U.S.. Journal of the American Statistical Association, 94, 1074–1082. FAY, R.E., HERRIOT, R.A., (1979). Estimation of income for small places: an application of James-Stein procedures to census data. Journal of the American Statistical Association, 74, 268–277. HODGES, J.S., SARGENT, D.J., (2001). Counting degrees of freedom in hierarchical and other richly parameterized models. Biometrika, 88, 367–379. HORVITZ, D.G., THOMPSON, D.J., (1952). A generalization of sampling without replacement from a finite universe. Journal of the American Statistical Association, 47, 663–685. LAHIRI, P., SUNTORNCHOST, J., (2015). Variable selection for linear mixed models with applications in small area estimation. Sankhya B, 1–9, doi = 10.1007/s13571-015-0096-0. NARAIN, R., (1951). On sampling without replacement with varying probabilities. Journal of the Indian Society of Agricultural Statistics, 3, 581–613. PFEFFERMANN, D., (2013). New important developments in small area estimation. Statistical Science, 28, 40–68. PFEFFERMANN, D., TILLER, R., (2006). Small Area Estimation with State- Space Models subject to Benchmark Constraints. Journal of the American Statistical Association, 101, 1387–1397. R DEVELOPMENT CORE TEAM, (2009). R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria. http://www.R-project.org. RAO, J.N.K., (2003). Small Area Estimation, New York: John Wiley. 540 Van den Brakel and Buelens: Covariate selection for small area estimation RAO, J.N.K., YU, M., (1994). Small-area estimation by combining time-series and cross-sectional data. The Canadian Journal of Statistics, 22, 511-528. SÄRNDAL, C-E., SWENSSON, B., WRETMAN, J., (1992). Model Assisted Survey Sampling, New York: Springer. SCHOUTEN, B., VAN DEN BRAKEL, J.A., BUELENS, B., VAN DER LAAN, J., KLAUSCH, T., (2013). Disentangling mode-specific selection bias and measurement bias in social surveys. Social Science Research, 42, 1555- 1570. VAIDA, F., BLANCHARD, S., (2005). Conditional Akaike information for mixedeffects models. Biometrika, 92, 351–370. YOU, Y., ZHOU, Q., (2011). Hierarchical Bayes small area estimation under a spatial model with application to health survey data. Survey Methodology, 37, 25–36.