Domain prediction with grouped income data

Walter, Paul; Groß, Marcus; Schmid, Timo; Tzavidis, Nikos Article — Published Version Domain prediction with grouped income data Journal of the Royal Statistical Society: Series A (Statistics in Society) Provided in Cooperation with: John Wiley & Sons Suggested Citation: Walter, Paul; Groß, Marcus; Schmid, Timo; Tzavidis, Nikos (2021) : Domain prediction with grouped income data, Journal of the Royal Statistical Society: Series A (Statistics in Society), ISSN 1467-985X, Wiley, Hoboken, NJ, Vol. 184, Iss. 4, pp. 1501-1523, https://doi.org/10.1111/rssa.12736 This Version is available at: https://hdl.handle.net/10419/284819 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. http://creativecommons.org/licenses/by/4.0/ J R Stat Soc Series A. 2021;184:1501–1523. | 1501 wileyonlinelibrary.com/journal/rssa Received: 29 April 2020 | Accepted: 25 June 2021 DOI: 10.1111/rssa.12736 ORIGINAL ARTICLE Domain prediction with grouped income data PaulWalter1 | MarcusGroß1 | TimoSchmid2 | NikosTzavidis3 This is an open access article under the terms of the Creat ive Commo ns Attri bution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2021 The Authors. Journal of the Royal Statistical Society: Series A (Statistics in Society) published by John Wiley & Sons Ltd on behalf of Royal Statistical Society. 1Institute of Statistics and Econometrics, Freie Universität Berlin, Berlin, Germany 2Institute of Statistics, Otto- Friedrich- Universität Bamberg, Bamberg, Germany 3Southampton Statistical Sciences Research Institute and Department of Social Statistics & Demography, University of Southampton, Southampton, UK Correspondence Timo Schmid, Institute of Statistics, Otto- Friedrich- Universität Bamberg, Bamberg, 96052, Germany. Email: [email protected] Abstract One popular small area estimation method for estimating poverty and inequality indicators is the empirical best predictor under the unitlevel nested error regression model with a continuous dependent variable. However, parameter estimation is more challenging when the response variable is grouped due to data confidentiality concerns or concerns about survey response burden. The work in this paper proposes methodology that enables fitting a nested error regression model when the dependent variable is grouped. Model parameters are then used for small area prediction of finite population parameters of interest. Model fitting in the case of a grouped response variable is based on the use of a stochastic expectation– maximization algorithm. Since the stochastic expectation– maximization algorithm relies on the Gaussian assumptions of the unitlevel error terms, adaptive transformations are incorporated for handling departures from normality. The estimation of the mean squared error of the small area parameters is facilitated by a parametric bootstrap that captures the additional uncertainty due to the grouping mechanism and the possible use of adaptive transformations. The empirical properties of the proposed methodology are assessed by using modelbased simulations and its relevance is illustrated by estimating deprivation indicators for municipalities in the Mexican state of Chiapas. KEYWORDS data confidentiality, intervalcensored data, nested error regression model, small area estimation, survey response burden 1502 | WALTER ET AL. 1 | INTRODUCTION Recent applications of small area estimation (SAE) methodologies have been concerned with the estimation of areaspecific income indicators, for example the median income, the head count ratio (HCR) and the Gini coefficient (Rao & Molina, 2015; Rojas- Perilla etal., 2020; Tzavidis etal., 2018). Popular SAE methods that have been used in this context include the socalled World Bank method (Elbers etal., 2003) and the empirical best predictor (EBP) method (Molina & Rao, 2010). In these papers, SAE is based on the use of a unitlevel nested error regression (random effects) model estimated with income or consumption as a response variable that is measured on a continuous scale. It is tempting for survey designers to reduce survey related costs by collecting information on income using income bands as opposed to detailed income information (Micklewright & Schnepf, 2010). Collecting data in bands may also help with reducing respondent burden, item nonresponse and microdata disclosure risk. On the other hand, it is also reasonable to expect that collecting grouped data may result in a loss of information compared to collecting on a continuous scale. The impact of this loss of information on the quality of official statistics estimates is of particular importance. There are several surveys and censuses that collect grouped income data, for example, the household and land survey (HLS) of Japan (Statistics of Japan, 2013), the German Microcensus (Statistisches Bundesamt, 2018) and the censuses of Australia (Australian Bureau of Statistics, 2011), Colombia (Departamento Administrativo Nacional De Estadística, 2005), and New Zealand (Statistics New Zealand, 2013). In the United Kingdom, the Office for National Statistics experimented with the collection of grouped income data in the lead up to the 2001 census (Collins & White, 1996). Using statistical methods for grouped data is not a problem specific to SAE. In particular, regression methods for grouped data have been studied in the econometric literature (Hsiao, 1983) but to the best of our knowledge, these methods have not been extended to include random effects. An alternative approach is to view the response as discrete and use generalized linear mixed models. For example, the response can be viewed as a multicategory or an ordinal outcome with cutoff points defined in the latter case by the grouping structure relevant to the dataset of interest. In this case one can motivate the model by assuming the existence of an underlying continuous latent variable, which although different has some similarities to the approach we propose in this paper. The emphasis in this paper is on modelbased small area inference more specifically, on estimating not only linear but also nonlinear indicators that are functions of the continuous (latent in the case of grouping) response variable. Hence, we propose an extension of the EBP method when the response variable is grouped. The methodology works by reversing the process of grouping, leading to an outcome measured on a continuous scale which is then used for areaspecific prediction. Estimation of the parameters of the unitlevel nested error regression model is implemented via a stochastic expectation– maximization (SEM) algorithm (Celeux & Diebolt, 1985). The method we propose in this paper is not only applicable in the case of SAE. The SEM algorithm can be used to estimate the model parameters when the dependent variable is grouped and interest is in using the model for drawing substantive conclusions about the relationship between the dependent variables and the explanatory variables (Walter, 2019a). The proposed methodology also allows for the use of datadriven transformations when diagnostic analyses indicate departures from the model assumptions. Using transformations as part of smallarea estimation when the continuous outcome is fully available has been already proposed in the literature. In the context of arealevel models, there are several papers discussing fixed transformations (e.g. Slud and Maiti (2006)) and datadriven transformations (e.g. Sugasawa and Kubokawa (2017)). Rojas- Perilla etal. (2020) presented theoretical and numerical justifications for the use of datadriven transformations with unitlevel SAE models. In particular, they propose an EBP approach with datadriven transformations where the datadriven transformation parameter is estimated by likelihoodbased methods. Following Gonzalez- Manteiga etal. (2008), the estimation of the mean squared error (MSE) of the small area estimates— when the response variable is grouped— is facilitated by a parametric bootstrap. | 1503WALTER ET AL. This incorporates the additional uncertainty due to grouping of the response variable assuming that the censoring mechanism is known. The proposed method assumes that there is no measurement error in reporting the group associated with the latent continuous variable. In this paper, we develop the methodology under a twolevel nested error regression model. However, an extension to threelevel structures— incorporating possible cluster effects— along the lines of the methodology proposed by Marhuenda etal. (2017) is feasible. Finally, as is the case with the EBP method or the World Bank method, we assume access to microdata for the model covariates from census or administrative data. The proposed methodology makes the use of SAE methods with grouped outcomes possible and therefore it enables survey organizations to consider collecting data in this form. The paper is organized as follows. Section 2 presents the survey data we use in this paper and defines the indicators of interest. The EBP approach and the nested error regression model when the response variable is available on a continuous scale are discussed in Section 3. Section 4 introduces the SEM algorithm that is used for the estimation of the model parameters when the response variable is grouped. In Section 5, the EBP method with grouped data is presented. In Section 6, modelbased simulations are carried out. In Section 7, the proposed methodology is used to estimate poverty and inequality indicators from grouped income data from Mexico. Finally, the main results are summarized in Section 8. 2 | ESTIMATING SMALL AREA DEPRIVATION INDICATORS FOR MUNICIPALITIES IN THE MEXICAN STATE OF CHIAPAS: DATA SOURCES AND INDICATORS We start with an initial discussion of the data and the poverty indicators of interest before presenting the methodological details. The aim of the proposed methodology is to enable the estimation of poverty and inequality indicators from survey data with a grouped income variable. To illustrate the proposed approach in this paper we use data from Mexico. Despite Mexico being the 15th largest economy in the world (International Monetary Fund, 2017), the fight against poverty and inequality is of great importance for the country since high poverty rates are omnipresent. During the Mexican peso crisis, extreme poverty increased from 21% in 1994 to 37% in 1996 (Pereznieto, 2010). Today, poverty rates remain at considerably high levels. According to the World Bank (2010), 33% of the population in the country experienced moderate poverty and 9% extreme poverty in 2013. This demonstrates the relevance of estimating and mapping poverty at local levels such that appropriate interventions can be designed. The poverty indicators we are interested in include the area HCR, the area poverty gap (PGAP) as defined in Foster etal. (1984) (income deprivation), the area average household income (Mean) and the Gini coefficient (Gini) in each area i. The indicators are defined as follows: HCRi=1 ni n i � j=1 I(yij ≤z), PGAP i=1 ni ni � j=1�z−yij z�I(yij ≤z ), Meani=1 ni ni � j=1 yij, Ginii= 2 ∑ ni j=1jyij ni ∑ ni j=1 yij −(ni+1) ni , 1504 | WALTER ET AL. where yij denotes the outcome variable, ni is the sample size, I(·) denotes the indicator function and z is the poverty threshold. In the simulations and application in this paper, z is set to 60% of the median of income, as defined by (Eurostat, 2014). When the income variable is measured on a continuous scale estimation of these indicators can be facilitated with standard SAE methods. However, when income is only available as grouped variable the methodology we propose in this paper can be used. For computing the incomedefined indicators of interest, one needs to have access to grouped income data that have been equivalized to account for different household sizes. If this has not been done, the secondary analyst will need access to household composition data in order to create equivalized grouped income data. Let us assume household i reports income in [2000,3000] and consists of two adults and one child leading to a weight of 2.5, then the household has an equivalized household income in the interval [800=2000/2.5,1200=3000/2.5]. This leads to household specific intervals depending on the weight relating to the household composition and reported interval. For the application in this paper we use the 2010 equivalized household income from the ENIGH (Encuesta Nacional de Ingreso y Gasto de los Hogares) survey and a large sample of the 2010 National Population and Housing census in Mexico. Both data sets are collected by the National Institute of Statistics and Geography (INEGI, Instituto Nacional de Estadística y Geografía) and they are provided to us by the National Council for the Evaluation of Social Development Policy (CONEVAL, Consejo Nacional de Evaluación de la Política de Desarrollo Social). Both the census and survey data sets include socioeconomic and regional information at household level. While the data cover all 31 states of Mexico, the application focuses on the state of Chiapas. Chiapas is one of the poorest states in Mexico with an average income of about 40% of the national median income (Levy etal., 2016). The state is located in the south of Mexico at the border to Guatemala. The survey covers 42 of the 118 municipalities in Chiapas. Hence, there are 76 out- of- sample municipalities for which no sample data are available. In order to derive reliable estimates at the level of municipality for all 118 municipalities, we rely on the use of modelbased methods and auxiliary information from the census and survey data. The sample size we analyse is n=2486 households and originally CONEVAL asked 3018 households in Chiapas. The response rate was 82%. The census sample size is N=96350 households. The regional distribution of the sample size is given in Table 1. The sample size of the insample municipalities varies between 13 and 651 households with a median sample size of 33 households. Since sample sizes are small in many municipalities, SAE methods can potentially improve the accuracy of direct small area estimates. In the next section a brief introduction to the EBP method when a continuous response variable is available. Then, the newly proposed methodology when a grouped response variable is available is introduced. 3 | EMPIRICAL BEST PREDICTION METHOD The target of inference are the small area parameters that include linear and nonlinear indicators which can be expressed as functions of an income variable, for example average and median equivalized income, the HCR, the PGAP and the Gini coefficient. Since in this paper we assume the availability of TABLE 1 Distribution of the sample and census household sizes across areas Min. 1st Qu. Median Mean 3rd Qu. Max. Sample 13.00 17.00 33.00 59.19 51.00 651.00 Census 82.00 399.50 617.50 816.50 839.00 7172.00 | 1505WALTER ET AL. unitlevel survey and census/administrative data, two methods for estimating nonlinear indicators in common use are available, the World Bank method (Elbers etal., 2003) and the popular EBP method (Molina & Rao, 2010). Although our focus is on the use of the EBP method, the proposed methodology can be applied in conjunction with the World Bank method too. The EBP method makes use of unitlevel nested error regression model and is summarized below. The response variable is income that is only available in the survey. The explanatory variables used for modelling the income variable are available both in the survey and in the census data sets. After the model is fitted using the survey data, the estimated model parameters are combined with census microdata to form unitlevel synthetic census predictions of the income variable. These synthetic values are then used for estimating the target parameters. Census predictions are generated by using the conditional predictive distribution of the out- of- sample data given the sample data. Although estimation of linear and nonlinear indicators can be also implemented with arealevel regression models (Fabrizi & Trivisano, 2016; Schmid etal., 2017), we focus on unitlevel models which can be used to produce estimates of a wide range of parameters as a byproduct of fitting the model. With arealevel models the focus is on one target parameter at the time. In addition, approaches to direct estimation with grouped data need to be carefully considered. Possible approaches to doing this are briefly outlined in the concluding remarks. Consider a finite population U of size N, divided into D areas/domains. The terms areas and domains are used interchangeably in this paper. The population size of each of the D- domains U1,U2,…,UD is given by N1,N2,…,ND. Let us for now assume that the response variable denoted by yij is measured on a continuous scale, where j=1, 2, …,ni denotes the jth unit belonging to the ith domain, with i=1,2,…,D. The vector x is defined as xT ij =(x1ij,…,xpij ) , where p denotes the number of explanatory variables. For each area i, the sample size is ni with n = ∑D i=1 n i and the population vector yi for area i comprises sampled and nonsampled units yT i =(y T is ,y T ir) . A nested error linear regression model is used for modelling the relationship between the variable of interest and auxiliary information with the unexplained variation being captured by the random effects term, ui and the residuals eij. In the simplest case, a twolevel nested error regression model as defined in Battese etal. (1988) is given by Assuming normality for the unitlevel error terms and the domain random effects, the conditional distribution of the out- of- sample data given the sample data is also normal. Predictions for the entire population of area i are generated from the following model, where ui=E(ui|yis) is the conditional expectation of ui given the sample data yis. Implementation of (2) requires replacing the unknown quantities 𝛽,𝜎u,𝜎e , with estimates and simulating L synthetic populations of the income variable, y∗ ij . Linear and nonlinear indicators are computed in each domain i for each replication and the estimates are averaged over the number of Monte Carlo simulations L. Following Molina and Rao (2010) and Rojas- Perilla etal. (2020) this number is usually set equal to L=50 or L=100 but higher numbers are also possible. (1) yij =x T ij 𝛽+ui+eij,(j=1, …,ni), (i=1, …,D ), ui iid ∼N(0, 𝜎2 u), e ij iid ∼N(0, 𝜎2 e ). (2) y∗ ij =x T ij 𝛽+ui+u ∗ i+e ∗ ij, u ∗ i iid ∼N(0, 𝜎2 u×(1−𝛾i)), e∗ ij iid ∼N(0, 𝜎2 e), 𝛾i=𝜎2 u 𝜎2 u+𝜎2 e ni , 1506 | WALTER ET AL. For the estimation of the unknown quantities 𝛽,𝜎u,𝜎e when the response variable is grouped we propose the use of a SEM algorithm. 4 | THE NESTED ERROR LINEAR REGRESSION MODEL WITH A GROUPED RESPONSE VARIABLE In the case of grouped data, yij is unobserved and the only observed information concerning the dependent variable is whether it falls within an interval. The continuous scale is divided into K intervals, where the kth interval is given by (Ak−1,Ak). The variable kij ∈{1, …,K} indicates in which of the intervals the dependent variable falls into. The first and K- th interval are allowed to be open ended, therefore A0= −∞ and AK= +∞ are possible. Situations in which both or none of the outer intervals are open ended can also be handled by the proposed methodology. Furthermore, the interval length is allowed to be arbitrary and can vary between intervals. Since the underlying distribution of yij is unknown, the aim is to reconstruct the conditional distribution f(yij |xij,kij,ui,𝜃) , where 𝜃 =( 𝛽 , 𝜎2 e , 𝜎2 u) are the unknown model parameters, β is a p×1 vector of regression coefficients and the random effects ui and the unitlevel error terms eij are assumed to be independent and normally distributed. Estimation methods such as maximum likelihood (ML) or restricted maximum likelihood (REML) are used for estimating θ when yij is observed on a continuous scale (Lindstrom & Bates, 1990). However, when the response variable is grouped, estimation of the parameters of interest is more challenging. The likelihood, ∏i∏j f(k ij � x ij ,u i , 𝜃) , cannot be derived directly, but can be expanded to include the latent yij into ∏i∏j f(k ij � y ij ,x ij ,u i , 𝜃 )×f(y ij � x ij ,u i , 𝜃) . While the second part is well known and can be maximized by the aforementioned methods, the first part, f(kij |yij,xij,ui,𝜃) , demands a more sophisticated estimation procedure as the latent part yij needs to be integrated out. In this section, an SEM algorithm for fitting the model is proposed and datadriven transformations are also considered for handling potential departures from the model assumptions. Before presenting the model and estimation method in detail, we review alternative approaches to dealing with grouped response variables and compare the SEM algorithm to alternative fitting methods. Different approaches for dealing with grouped response variables in regression modelling that assume independent observations have been proposed in the literature. A naive approach uses ordinary least squares on the midpoints of the intervals. While this approach is easy to implement (Thompson & Nelson, 2003), it has two major drawbacks. The uncertainty associated with the value of each observation within each interval is not accounted for and dealing with openended intervals is not easy. Nevertheless, the naive approach can provide results of acceptable quality if the grouping is very fine (Fryer & Pethybridge, 1972). An alternative approach is to view the response as discrete and use a generalized linear mixed model. Approaches to modelling multicategory discrete outcomes have been proposed in the small area literature (Lopez- Vizcaino etal., 2015; Molina etal., 2007). Nevertheless, one difficulty with the use of discretetype models remains. In our application we are not only interested in estimating the proportion of units in a category but also interested in estimating indicators such as the PGAP, the HCR and the Gini coefficient. Therefore, if we decide to use a discretetype model we also need to develop a method for recovering estimates of target parameters that are usually computed by using a continuous outcome. To overcome these drawbacks, linear regression models for leftcensored (Tobin, 1958), rightcensored (Rosett & Nelson, 1975) and grouped (or intervalcensored) (Stewart, 1983) variables have been proposed. Stewart (1983) proposes an expectation– maximization (EM) algorithm for estimating the model parameters of a linear regression model with a grouped response variable. While the original paper introducing the EM algorithm (Dempster etal., 1977) proposed maximizing the likelihood within the M- step, it mentioned that the EM algorithm can | 1507WALTER ET AL. be also used to obtain REML estimates. Literature related to this includes Kim and Taylor (1995) and Foulley etal. (2000). To estimate the parameters of the nested error regression model when the outcome is grouped, we propose the use of a SEM algorithm (Celeux & Diebolt, 1985; Celeux etal., 1996). A similar SEM algorithm is proposed in Groß etal. (2017) for kernel density estimation on aggregated data. SEM can be regarded as a middle ground between the EM and full MCMC. With the EM algorithm we alternate between calculating the expectation of the conditional distribution f(yij |xij,kij,ui,𝜃) (E- Step) and obtaining θ via maximizing the complete data likelihood (M- Step). However, with fixed intervals (Ak−1,Ak) it can be seen that a bias in θ will be introduced, for example, 𝜎2 e would be underestimated as the overall variance of the expectations of yij is much smaller than that of the (latent) variable yij. SEM and full MCMC (Gelman etal., 2013) replace the E- Step by drawing from the conditional distribution of yij and therefore do not suffer from this drawback. This approach can be also viewed as part of the literature about measurement error models (Carroll etal., 2006), where the latent values yij are regarded as model parameters or partially observed data (Carpenter etal., 2012). In addition, MCMC also replaces the M- Step by sampling from the conditional distribution of θ. In summary, compared to EM, SEM avoids or reduces biases in the estimation of the parameters of interest, while compared to MCMC, SEM is considerably faster due to faster convergence because only the values yij are drawn. Using the SEM also saves time with the implementation because the users can make use of existing estimation algorithms for the M- Step, while it is easy to make draws of yij. Considering the assessment of convergence, SEM should be treated similarly to MCMC with its broad variety of convergence measures. Related to the SEM approach is also the simulated maximum likelihood (SML, Gouriéroux and Monfort, 1990) method. The SML also samples yij values but uses these samples to estimate the expectation of the density f(kij |yij,xij,ui,𝜃)×f(yij |xij,ui,𝜃) which is then maximized with respect to θ. However, SML is not unbiased, but it is consistent (Gouriéroux & Monfort, 1990), and not as straightforward to implement as SEM as one needs to deal with the numerical aspects of the optimization procedure. Let us now consider the model we use in this paper. To reconstruct the unknown distribution f(yij |xij,kij,ui,𝜃) we use the Bayes theorem and express the target distribution as follows: To avoid confusion, note that in the notation we use here we treat 𝜎2 u as part of θ. However, when writing the distribution of the random effects we make it explicit that this distribution depends on 𝜎2 u . Since f(kij |yij,xij,ui,𝜃)=f(kij |yij) , the conditional distribution of kij is given by and under the nested error regression model (1), 4.1 | The SEM algorithm Because yij and ui are unobserved, one approach to fitting the model defined above is to use the SEM algorithm. Generally speaking, the algorithm works by replacing the unobserved response data yij in the complete data likelihood by generating pseudo samples of the unobserved response data given the f(yij |xij,kij,ui,𝜃)∝f(kij |yij,xij,ui,𝜃)f(yij |xij,ui,𝜃). f(kij | yij)= { 1 if Akij−1≤yij ≤Akij , 0 else, f(yij | xij,ui, 𝜃 )∼N(x T ij 𝛽 +ui, 𝜎2 e ), f(ui |𝜎2 u )∼N(0, 𝜎2 u). 1508 | WALTER ET AL. observed data and the current values of θ (S- step) and then maximizes the complete data likelihood for updating θ and the predicted random effects in the M- step. The iterations stop after B+M steps. Assuming θ is known, pseudo samples, yij , are drawn from the following conditional distribution where I(·) denotes the indicator function. The conditional distribution of yij has the form of a twosided truncated normal distribution given by with 𝜇 ij =x T ij 𝛽 +u i , ϕ(·) is the probability density function of the standard normal distribution and Φ(·) is its cumulative distribution function. By definition Φ(A kij −𝜇 ij 𝜎 e) = 1 if A k ij = +∞ and Φ(A kij−1 −𝜇 ij 𝜎 e) = 0 if A k ij −1 = −∞ . For each observation with explanatory variables xij the corresponding yij is randomly drawn from N (x T ij 𝛽 +ui, 𝜎2 e) within the given interval A k ij −1 ≤y ij ≤A k ij . This is the S- step of the SEM algorithm. The M- step comprises fitting the nested error regression model using the newly generated (yij,xij) . The steps of the SEM algorithm are as follows: 1. Estimate  𝜃 =(  𝛽,𝜎 2 e ,𝜎 2 u) and u from (1) using the midpoints of the intervals as a substitute for the unknown yij. The parameters are estimated using REML and using the empirical best linear unbiased predictor (EBLUP) for u. 2. S- step: For j=1, …,ni and i=1,…,D sample from the conditional distribution f(yij |xij,kij,ui,𝜃) by drawing randomly from N (xT ij  𝛽+ui,𝜎 2 e) within the given interval A k ij −1 ≤y ij ≤A k ij obtaining (yij,xij) . The drawn pseudo yij are used as replacement for the unknown yij. 3. M- step: Reestimate the model parameters and the predicted random effects using (1) and the pseudo samples (yij,xij) from Step 2. The parameters are estimated as in Step 1. 4. Iterate Steps 2– 3 B+M times, with B burnin iterations and M additional iterations. 5. Discard the burnin iterations and estimate  𝜃 by averaging the derived M estimates. For openended intervals A0= −∞ and AK= +∞ , the midpoints M1 and MK in Step 1 are computed as follows: where Note that Step 1 is purely for the initialization of the algorithm. Empirical results show that using the midpoints of the intervals as a substitute for the unknown yij and the procedure for handling openended intervals in the first iteration step has little impact on the estimates. Empirical results are provided in table 9 in the online supplementary material (OSM). f(yij | xij,kij,ui, 𝜃 )∝I(Ak ij −1≤yij ≤Ak ij )×N(x T ij 𝛽 +ui, 𝜎2 e ), f(yij | xij,kij,ui,𝜃)= 𝜙 ( yij − 𝜇 ij 𝜎e ) 𝜎e ( Φ ( Akij −𝜇ij 𝜎 e) −Φ ( Akij−1−𝜇ij 𝜎 e)), M1=(A1−D)∕2, MK=(AK−1+D)∕2, D =1 (K−2) K−1 ∑ k=2| Ak−1−Ak |. | 1515WALTER ET AL. estimation results of λ in table 3 in the OSM. Finally, Table 3 shows that the proposed bootstrap MSE estimator has reasonably low relative bias. As expected, the MSE under the Box- Cox version of the SEM is somewhat more unstable than the corresponding MSE for the Log SEM. This is due to the fact that in the case of the Box- Cox method the transformation parameter is estimated with each bootstrap sample whereas for the Log method the transformation is held fixed. In order to further evaluate the impact of the grouping on the performance of the SEM estimators, we report as part of the OSM two additional Logscale scenarios. For the first one we used 14 equally spaced intervals leading to a large proportion of observations in the upper openended interval. For the second scenario we increased the interval size with increasing y values. The results are reported in table 12 in the OSM. 7 | ESTIMATING SMALL AREA DEPRIVATION INDICATORS FOR MUNICIPALITIES IN THE MEXICAN STATE OF CHIAPAS: AN APPLICATION OF THE SEM BOXCOX METHOD In our application the response variable, equivalized household income, is measured on a continuous scale. In order to assess the performance of the proposed methodology, we group equivalized TABLE 3 Performance of the bootstrap root mean squared error (MSE) estimator over areas Mean HCR PGAP Gini Indicator: Median Mean Median Mean Median Mean Median Mean Normal scenario 1 (14 intervals) rel.Bias[%] SEM 7.37 7.05 5.91 5.07 2.31 3.07 3.90 3.88 SEM Box- Cox 7.56 7.33 5.61 5.47 −6.86 −5.58 −3.67 −4.41 rel.RMSE[%] SEM 9.50 10.50 10.59 11.38 12.05 13.34 8.68 9.87 SEM Box- Cox 9.92 10.85 10.78 11.42 13.03 14.01 8.81 10.53 Normal scenario 2 (7 intervals) rel.Bias[%] SEM 5.30 5.84 4.71 3.65 −0.18 0.30 2.24 1.94 SEM Box- Cox 5.46 6.10 4.59 3.91 −15.77 −14.82 −10.22 −11.67 rel.RMSE[%] SEM 8.59 9.91 10.22 10.98 12.22 13.29 8.84 9.59 SEM Box- Cox 8.82 10.30 10.30 11.07 19.27 19.16 12.09 14.79 Logscale scenario (7 intervals) rel.Bias[%] SEM Log 7.22 6.56 6.73 7.58 6.74 7.13 0.88 0.77 SEM Box- Cox 13.17 26.00 6.78 7.65 7.10 7.57 6.54 6.40 rel.RMSE[%] SEM Log 33.49 34.78 13.19 14.25 21.02 21.63 7.95 8.36 SEM Box- Cox 42.19 60.85 13.23 14.36 21.33 21.93 16.49 17.05 1516 | WALTER ET AL. household income to 14 and 8 intervals. The distribution of the grouped equivalized household income is presented in tables 16 and 17 of the OSM. The variables in Table 4 were identified as possible covariates that predict equivalized household labour income well. The variables in the working model are selected by using the coefficient of determination proposed by Nakagawa and Schielzeth (2013). The conditional R2 c , interpreted as the variance explained by the whole model, is R2 c,lme = 0.61 when estimating the model with the observed continuous response variable on the transformed scale (Box- Cox transformation). When estimating the model with a grouped response variable on the transformed scale (Box- Cox transformation) using the SEM algorithm the R2 c,sem(14) is 0.61 and R 2 c,sem(8) is 0.62 for the 14 and 8 interval scenario respectively. The Box- Cox transformation is used as the preferred transformation method because it is datadriven. This is crucial when working with grouped data as response variable, because the normality assumption of the residuals cannot be checked. The estimated transformation parameters are  𝜆lme =0.16 for the continuous response,  𝜆(F) sem(14) = 0.18 and  𝜆(F) sem(8) = 0.17 for the 14 and 8 grouping scenarios respectively. The results indicate that the use of a logarithmic transformation or the use of the untransformed response variable may lead to erroneous results. Rojas- Perilla etal. (2020) and Tzavidis etal. (2018) show that even if λ is estimated to be close to 0 the EBP estimates using the Box- Cox transformation may outperform the EBP estimates using the logarithmic transformation. Estimates of the mean equivalized household labour income, HCR and PGAP for each of the 118 municipalities are obtained by using the SEM Box- Cox method based on 14 and 8 intervals, and by using LME Box- Cox based on the observed continuous response variable. The mean and median TABLE 4 Variables used in the nested error regression working model Variable type Description Response variable: Grouped equivalized household labour income Auxiliary variables: Value of all household goods Value of household communication equipment Share of employees in the household Educational level of head of household Social class of head of household Municipalities of Chiapas FIGURE 2 Estimated head count ratio (HCR) for municipalities in the state of Chiapas based on different estimation methods. The empirical best predictor (EBP) under the model with continuous response and Box- Cox transformation is abbreviated by LME Box- Cox, the EBP under the model with grouped data and Box- Cox transformation is abbreviated by SEM Box- Cox [Colour figure can be viewed at wileyonlinelibrary.com] LMEBox-Cox − HCR 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 SEM Box-Cox 14 Intervals − HCR 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 SEMBox-Cox 8 Intervals − HCR 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 | 1517WALTER ET AL. TABLE 5 Point estimates and corresponding relative efficiencies of the point estimates (EFF) over municipalities using the stochastic expectation– maximization (SEM) Box- Cox algorithm Box- Cox Mean HCR PGAP Gini Median Mean Median Mean Median Mean Median Mean Point est. LME 814.4 872.6 0.426 0.432 0.220 0.233 0.535 0.539 Point est. SEM 14 intervals 812.1 870.8 0.426 0.427 0.224 0.233 0.531 0.536 EFF 0.983 0.995 1.018 1.024 1.022 1.035 1.057 1.064 Point est. SEM 8 intervals 810.5 863.9 0.421 0.428 0.219 0.232 0.529 0.534 EFF 1.028 1.013 1.029 1.038 1.043 1.046 1.071 1.064 The empirical best predictor (EBP) under the model with continuous response and Box- Cox transformation is abbreviated by LME Box- Cox, the EBP under the model with grouped data and Box- Cox transformation is abbreviated by SEM Box- Cox. HCR stands for the head count ratio and PGAP for the poverty gap 1518 | WALTER ET AL. averaged over all municipalities are given in Table 5 and plotted in Figures 2 and 4. The results show that the point estimates from all three estimation methods are very similar. Interval censoring does not appear to impact significantly on the estimation results. Additionally, the relative efficiencies of the estimators (EFF) defined as EFF ( I EBP i )=RMSE sem ( I EBP i )∕RMSE lme ( I EBP i) is reported in the Table 5. It is notable that the efficiency loss is small even when the response variable is grouped to only eight intervals. In the 14 interval scenario the point estimates of the mean are even more efficient, but this result is only due to the Monte Carlo variability. The spatial distributions of the HCR in municipalities in Chiapas are shown in Figure 2 for all three estimation methods. The figure supports the previously mentioned results that the estimates obtained by using the different methods are very close. A possible way to further validate the estimation results is by comparing the direct estimates, where available, to the modelbased estimates. In Figure 3 the direct estimates of the mean (based on the observed continuous data) are compared to the modelbased estimates (SEM Box- Cox) of the mean using a grouped response variable with 14 intervals. As expected, the left panel shows a positive linear correlation between the estimates. However, there is a disparity between the intersection line (the identity) and the regression line. As anticipated, the modelbased estimates are less extreme compared to the direct estimates FIGURE 3 The left panel shows a scatter plot and the right panel a line plot of the direct and the modelbased estimates for each insample domain (municipality). The empirical best predictor under the model with grouped data and Box- Cox transformation is abbreviated by SEM Box- Cox [Colour figure can be viewed at wileyonlinelibrary. com] 0 1000 2000 3000 0 1000 2000 3000 Direct SEM Box−Cox Intersection Reg. line Mean 0 1000 2000 3000 010203040 Municipalities (ordered by sample size) Value Method Direct SEM Box−Cox Mean FIGURE 4 Estimated mean and poverty gap for municipalities in the state of Chiapas [Colour figure can be viewed at wileyonlinelibrary.com] SEM Box-Cox 14 intervals Mean 500 1000 1500 2000 2500 SEM Box-Cox 14 intervals PGAP 0.1 0.2 0.3 0.4 0.5 | 1519WALTER ET AL. for municipalities with very small and very high mean estimates. The right panel plots the value of the estimates for both estimation methods for each insample domain. The pattern shows that as the sample size increases the direct estimates and the SEM Box- Cox estimates are almost identical. Figure 4 presents municipal estimates of mean income and PGAP for the SEM Box- Cox algorithm based on 14 intervals. The plots for the other estimation methods are omitted because the results are comparable. We observe that municipalities in the middle and in the east of Chiapas exhibit high rates of HCRs and PGAPs and low levels of mean equivalized household labour income and are thus more adversely affected by poverty. These regions are characterized by high mountains, the Chiapas Highlands and a large concentration of indigenous population. There are, however, two regions in the centre of the state with relatively high mean income and low rates of poverty. These are the regions where the capital Tuxtla Gutiérrez and the larger city San Cristóbal de las Casas are located. Also the coastal region— especially in the south— where the most important city economically Tapachula is located, is less affected by poverty. The analysis shows that even though Chiapas is one of the poorest states in Mexico, there are spatial variations between the municipalities. These differences can be revealed by using SAE methods designed for grouped data. The proposed SEM Box- Cox method is, to the best of our knowledge, the first approach that allows the use of the popular EBP method in conjunction with a grouped response variable. This enables the estimation of spatially disaggregated target indicators with small sample sizes when confidentiality restrictions or decisions about the survey design require the use of relatively limited information for the response variable. 8 | CONCLUDING REMARKS The paper proposes SAE methodology when working with a response variable that is grouped. The novel aspects of the paper include the estimation of a nested error regression model when the response is grouped, the estimation both of linear and nonlinear indicators for small areas, the use of datadriven transformation with the nested error regression model and the estimation of the MSE of the small area target parameters that accounts for the fact that we are working with limited information compared to standard small area models. The proposed methods are evaluated using modelbased simulations under different scenarios for the unitlevel error terms. The results show that the proposed methods work well and in most scenarios the loss of accuracy is small when compared to the use of EBPs that are estimated by assuming the availability of full information for the response variable. As expected, the loss of accuracy also depends on the number of intervals used for grouping the data and the proposed methodology appears to work well even when the number of groups used is fairly small. The results also show that the use of an adaptive transformation works satisfactorily and the transformation parameter is estimated well in the presence of limited information for the response variable. Finally, the proposed MSE estimator appears to capture the different sources of variability and appropriately tracks the empirical MSE. The new methodology is used to estimate disaggregated poverty and inequality indicators for municipalities in Chiapas, a southern state of Mexico, using grouped income banded in 8 and 14 intervals. In order to evaluate the proposed methodology estimates of the target parameters are also obtained when income is fully available, that is, not grouped. The Box- Cox transformation is applied to ensure that the model assumptions are met. The estimates from the continuous and grouped responses are very close, indicating the validity of the proposed methodology. The plotted poverty maps enable policy makers to get a spatial overview of the distribution of poverty in Chiapas and to target poorer regions more precisely. The proposed methodology for estimating nonlinear indicators with grouped response data requires access to unitlevel census microdata for the covariates. Access to such data may be very 1520 | WALTER ET AL. challenging due to confidentiality constraints. Although the proposed methodology assumes access to unitlevel census microdata, it is important to discuss briefly an alternative when such data are not available. An alternative approach would be to use arealevel models which are based on direct estimates of the linear or nonlinear indicator of interest. Methods for direct estimation with grouped data can be mainly categorized in three groups: (1) Direct estimation based on the midpoints of the intervals, (2) parametric (Chen, 2017; Reed & Wu, 2008) and (3) nonparametric (Kakwani & Podder, 2008) modelling of the distribution function. How these different direct estimation methods can be combined with arealevel models is an open area for further research. The proposed modelbased small area methodology does not incorporate survey weights in the estimation. Conventionally, SAE methods are model based and in most cases the survey weights are not used in model fitting. However, not including the survey weights carries risks. One example is when the assumption of a noninformative sample selection mechanism does not hold, even after conditioning on auxiliary variables, hence wrongly assuming that the model for the sample also holds for the population. An approach to accounting for the survey weights in EBP was recently proposed by Guadarrama etal. (2018). Although not implemented in this paper, the pseudo EBP approach can be adapted to the setting of the present paper. Doing so requires fitting the nested error regression model and estimating the fixed effects and the variance components in each step of the SEM algorithm by using the methods in You and Rao (2002). The issue of missing data in SAE has received some attention in the literature. Similarly to the use of survey weights, most small area literature assumes a missing at random mechanism hence after conditioning on covariates, the probability to respond is assumed not to depend on the response. This is the assumption we are making in the present paper. Provided that the survey weights adjust for nonresponse, one approach to account for nonresponse is by incorporating the weights in model fitting as proposed by Guadarrama etal. (2018). In a recent paper Sverchkov and Pfeffermann (2018) proposed an alternative approach to modelling nonmissing at random mechanisms in SAE. However, to the best of our knowledge this has not been extended yet for estimating the more general parameters that are of interest in this paper. Current research focuses on extending the SEM method for fitting nested error regression models for more complex structures, for example models with random coefficients. In future research, we also plan to focus on the case where grouping also affects some of the auxiliary variables. This is a more challenging problem but perhaps more realistic if interest is in protecting data confidentiality. Finally, there are three further aspects that we do not discuss in this paper and are left open for future research. First, the proposed methodology does not adjust for the effects of heaping. This can be resolved by following the methods proposed by Groß and Rendtel (2016). Second, we also acknowledge that another type of measurement error may exist when respondents report their income in the wrong interval. However, this is a type of misclassification error that cannot be solved unless we are willing to impose additional assumptions or use results from a validation sample, which can be treated as a gold standard. Third, the proposed methodology assumes normality for the random effects. The assumption may be relaxed by leaving the distribution of the random effects unspecified and use nonparametric methods (Marino etal., 2019). This leads to a discrete mixture distribution which avoids the need to impose parametric assumptions like normality for the random effects. However, the extension to the case of grouped response variables is a topic for further research. ACKNOWLEDGEMENTS The research has received funding from the European Unions Horizon 2020 research and innovation programme under grant agreement No 730998, InGRID- 2, Integrating Research Infrastructure | 1521WALTER ET AL. for European expertise on Inclusive Growth from data to policy. Furthermore, Schmid and Tzavidis gratefully acknowledge support by grant ES/N011619/1 - Innovations in Small Area Estimation Methodologies from the UK Economic and Social Research Council and funding under the Data and Evidence to End Extreme Poverty (DEEP) research programme. DEEP is a consortium of the Universities of Cornell, Copenhagen, and Southampton led by Oxford Policy Management, in partnership with the World Bank's Development Data Group and funded by the UK Foreign, Commonwealth & Development Office. The authors are grateful to the National Council for the Evaluation of Social Development Policy (CONEVAL, Consejo Nacional de Evaluación de la Política de Desarrollo Social) for providing the data used in empirical work. The views set out in this paper are those of the authors and do not reflect the official opinion of CONEVAL. The numerical results are not official estimates and are only produced for illustrating the methods. Finally, the authors are indebted to the Joint Editor, Associate Editor and two referees for comments that significantly improved the paper. ORCID Timo Schmid http://orcid.org/0000-0002-7217-2501 REFERENCES Australian Bureau of Statistics (2011) Census household form. Available from https://unsta ts.un.org/unsd/demog raphi c/sourc es/censu s/quest/ AUS20 11en.pdf. Accessed data 2018- 04- 05. Battese, G.E., Harter, R.M. & Fuller, W.A. (1988) An error component model for prediction of county crop areas using survey and satellite data. Journal of the American Statistical Association, 83(401), 28– 36. Box, G.E.P. & Cox, D.R. (1964) An analysis of transformations. Journal of the Royal Statistical Society: Series B, 26(2), 211– 252. Carpenter, J.R., Goldstein, H. & Kenward, G.M. (2012) Statistical modelling of partially observed data using multiple imputation: principles and practice, Netherlands, Dordrecht: Springer, 15– 31. Carroll, R.J., Ruppert, D., Stefanski, L. A., & Crainiceanu, C. M. (2006) Measurement error in nonlinear models: A modern perspective. Boca Raton: CRC Press. Celeux, G. & Diebolt, J. (1985) The sem algorithm: A probalistic teacher algorithm derived from the em algorithm for the mixture problem. Computational Statistics Quarterly, 2:73– 82. Celeux, G., Chauveau, D. & Diebolt, J. (1996) Stochastic versions of the EM algorithm: An experimental study in the mixture case. Journal of Statistical Computation and Simulation, 55(4), 287– 314. Chen, Y.T. (2017) A unified approach to estimating and testing income distributions with grouped data. Journal of Business & Economic Statistics, 36(3), 1– 18. Collins, D. & White, A. (1996) In search of an income question for the 2001 census. Survey Methodology Bulletin, 39(7), 2– 10. Dempster, A., Laird, N. & Rubin, D. (1977) Maximum likelihood from incomplete data via the em algorithm. Journal of the Royal Statistical Society: Series B (Methodological), 39(1), 1– 38. Departamento Administrativo Nacional De Estadística (2005) Censo general 2005. Available from https://www.dane. gov.co/files/ censo s/libro Censo 2005n acion al.pdf?&. Accessed date 2018- 04- 05. Draper, N.R. & Cox, D.R. (1969) On distributions and their transformation to normality. Journal of the Royal Statistical Society: Series B (Methodological), 31(3), 472– 476. Elbers, C., Lanjouw, J.O. & Lanjouw, P. (2003) Microlevel estimation of poverty and inequality. Econometrica, 71(1), 355– 364. Eurostat (2014) Statistics explained: Atrisk- ofpoverty rate. Available from http://ec.europa.eu/euros tat/stati stics - expla ined/index.php/Gloss ary:Atrisk- ofpover ty_rate. Accessed date 2018- 05- 30. Fabrizi, E. & Trivisano, C. (2016) Small area estimation of the gini concentration coefficient. Computational Statistics & Data Analysis, 99:223– 234. Foster, J., Greer, J. & Thorbecke, E. (1984) A class of decomposable poverty measures. Econometrica, 52(3), 761– 766. 1522 | WALTER ET AL. Foulley, J.- L., Jaffrézic, F. & Robert- Granié, C. (2000) Emreml estimation of covariance parameters in gaussian mixed models for longitudinal data analysis. Genetics Selection Evolution, 32(2), 129. Fryer, J.G. & Pethybridge, R.J. (1972) Maximum likelihood estimation of a linear regression function with grouped data. Journal of the Royal Statistical Society: Series C, 21(2), 142– 154. Gelman, A., Carlin, J. B., Stern, H.S., Dunson, D. B., Vehtari, A. & Rubin, D. B. (2013) Bayesian data analysis. Boca Raton: CRC Press. Gonzalez- Manteiga, W., Lombardia, M.J., Molina, I., Morales, D. & Santamaria, L. (2008) Analytic and bootstrap approximations of prediction errors under a multivariate fayherriot model. Computational Statistics & Data Analysis, 52 (12), 5242– 5252. Gouriéroux, C. & Monfort, A. (1990) Simulation based inference in models with heterogeneity. Annals of Economics and Statistics, 20– 21, 69– 107. Groß, M. & Rendtel, U. (2016) Kernel density estimation for heaped data. Journal of Survey Statistics and Methodology, 4(3), 339– 361. Groß, M., Rendtel, U., Schmid, T., Schmon, S. & Tzavidis, N. (2017) Estimating the density of ethnic minorities and aged people in Berlin: Multivariate kernel density estimation applied to sensitive georeferenced administrative data protected via measurement error. Journal of the Royal Statistical Society: Series A, 180(1), 161– 183. Guadarrama, M., Molina, I. & Rao, J. (2018) Small area estimation of general parameters under complex sampling designs. Computational Statistics & Data Analysis, 121, 20– 40. Gurka, M.J., Edwards, L.J., Muller, K.E. & Kupper, L. (2006) Extending the boxcox transformation to the linear mixed model. Journal of the Royal Statistical Society: Series A, 169(2), 273– 288. Hsiao, C. (1983) Studies in econometrics, time series and multivariate statistics, chapter regression analysis with a categorized explanatory variable, Cambridge, Massachusetts: Academic Press, 93– 129. International Monetary Fund (2017) World economic outlook database. Available from http://www.imf.org/exter nal/ pubs/ft/weo/2017/01/weoda ta/index.aspx. Accessed date 2017- 10- 14. Kakwani, N.C. & Podder, N. (2008) Efficient estimation of the Lorenz curve and associated inequality measures from grouped observations Lorenz curve and associated inequality measures from grouped observations. In: Chotikapanich, D. (Ed.) Modeling income distributions and lorenz curves. New York: Springer, pp. 57– 70. Kim, D.K. & Taylor, J.M. (1995) The restricted em algorithm for maximum likelihood estimation under linear restrictions on the parameters. Journal of the American Statistical Association, 90(430), 708– 716. Kreutzmann, A.- K., Pannier, S., Perilla, N., Schmid, T., Templ, M. & Tzavidis, N. (2019) The R package emdi for estimating and mapping regionally disaggregated indicators. Journal of Statistical Software, 91(7), 1– 33. Levy, D., Hausmann, R., Santos, M. A., Espinoza, L. & Flores, M. (2016) Why is chiapas poor? Center for International Development at Harvard University Working Paper, (300) Lindstrom, M. J. & Bates, D. M. (1990) Nonlinear mixed effects models for repeated measures data. Biometrics, 46(3), 673– 687. Lopez- Vizcaino, E., Lombardia, M. J. & Morales, D. (2015) Small area estimation of labour force indicators under a multinomial model with correlated time and area effects. Journal of the Royal Statistical Society: Series A, 178(3), 535– 565. Marhuenda, Y., Molina, I., Morales, D. & Rao, J. N. K. (2017) Poverty mapping in small areas under a twofold nested error regression model. Journal of the Royal Statistical Society: Series A, 180(4), 1111– 1136. Marino, M. F., Ranalli, M. G., Salvati, N. & Alfò, M. (2019) Semiparametric empirical best prediction for small area estimation of unemployment indicators. Annals of Applied Statistics, 13(2), 1166– 1197. Micklewright, J. & Schnepf, S. (2010) How reliable are income data collected with a single question? Journal of the Royal Statistical Society: Series A, 173(2), 409– 429. Molina, I. & Rao, J.N.K. (2010) Small area estimation of poverty indicators. The Canadian Journal of Statistics, 38(3), 369– 385. Molina, I., Saei, A. & Jose Lombardia, M. (2007) Small area estimates of labour force participation under a multinomial logit mixed model. Journal of the Royal Statistical Society: Series A, 170(4), 975– 1000. Nakagawa, S. & Schielzeth, H. (2013) A general and simple method for obtaining r2 from generalized linear mixedeffects models. Methods in Ecology and Evolution, 4(2), 133– 142. Pereznieto, P. (2010) The case of mexico's 1995 peso crisis and argentina's 2002 convertibility crisis: Including children in policy responses to previous economic crises. UNICEF: Social and economic policy. Rao, J. & Molina, I. (2015) Small area estimation. Hoboken: John Wiley & Sons, Inc. | 1523WALTER ET AL. Reed, W.J. & Wu, F. (2008) New four- and fiveparameter models for income distributions. In: Chotikapanich, D. (Ed.) Modeling income distributions and lorenz curves. New York: Springer, pp. 211– 224. Rojas- Perilla, N., Pannier, S., Schmid, T. & Tzavidis, N. (2020) Datadriven transformations in small area estimation. Journal of the Royal Statistical Society: Series A, 183:121– 148. Rosett, R.N. & Nelson, F.D. (1975) Estimation of the twolimit probit regression model. Econometrica, 43(1), 141– 146. Schmid, T., Bruckschen, F., Salvati, N. & Zbiranski, T. (2017) Constructing sociodemographic indicators for national statistical institutes by using mobile phone data: Estimating literacy rates in Senegal. Journal of the Royal Statistical Society: Series A, 180(4), 1163– 1190. Slud, E. & Maiti, T. (2006) Meansquared error estimation in transformed fayherriot models. Journal of the Royal Statistical Society: Series B, 68(2), 239– 257. Statistics New Zealand (2013) New Zealand census of population and dwellings. Available from https://unsta ts.un.org/ unsd/demog raphi c/sourc es/censu s/quest/ NZL20 13enIn.pdf. Accessed date 2018- 05- 13. Statistics of Japan (2013) Survey results of the housing and land survey. Available from https://www.stat.go.jp/engli sh/ data/jyuta ku/resul ts.html. Accessed date 2021- 01- 13. Statistisches Bundesamt (2018) Der Mikrozensus stellt sich vor. Available from https://www.desta tis.de/DE/Zahle nFakt en/Gesel lscha ftSta at/Bevoe lkeru ng/Mikro zensus.html. Accessed date 2020- 12- 20. Stewart, M. (1983) On least square estimation when the dependent varaible is grouped. The Review of Economic Studies, 50(4), 737– 753. Sugasawa, S. & Kubokawa, T. (2017) Transforming response values in small area prediction. Computational Statistics and Data Analysis, 114:47– 60. Sverchkov, M. & Pfeffermann, D. (2018) Small area estimation under informative sampling and not missing at random nonresponse. Journal of the Royal Statistical Society: Series A, 181(4), 981– 1008. Thompson, M.L. & Nelson, K. (2003) Linear regression with type i interval- and leftcensored response data. Environmental and Ecological Statistics, 10(2), 221– 230. Tobin, J. (1958) Estimation of relationships for limited dependent variables. Econometrica, 26(1), 24– 36. Tzavidis, N., Zhang, L.- C., Luna, A., Schmid, T. & Rojas- Perilla, N. (2018) From start to finish: A framework for the production of small area official statistics. Journal of the Royal Statistical Society: Series A, 181(4), 927– 979. Walter, P. (2019a) A selection of statistical methods for intervalcensored data with applications to the german microcensus. Available from https://refub ium.fuberlin.de/bitst ream/handl e/fub18 8/23841/ Disse rtati on_Paul_Walter. pdf?seque nce=3&isAll owed=y. Walter, P. (2019b) smicd: Statistical Methods for Interval Censored Data. R package version 1.0.3. World Bank (2010) Poverty & equity data portal. Available from http://pover tydata.world bank.org/pover ty/count ry/ MEX/. Accessed date 2017- 10- 14. You, Y. & Rao, J.N.K. (2002) A pseudoempirical best linear unbiased prediction approach to small area estimation using survey weights. The Canadian Journal of Statistics/La Revue Canadienne de Statistique, 30(3), 431– 439. SUPPORTING INFORMATION Additional supporting information may be found online in the Supporting Information section. How to cite this article: Walter P, Groß M, Schmid T, Tzavidis N. Domain prediction with grouped income data. J R Stat Soc Series A. 2021;184:1501–1523. https://doi.org/10.1111/ rssa.12736