Semi-Parametric Empirical Best Prediction for small area estimation of unemployment indicators

TL;DR

This paper introduces a semi-parametric approach for small area estimation of unemployment rates, improving accuracy and computational efficiency by avoiding strict distributional assumptions and heavy integral calculations.

Contribution

It proposes a semi-parametric empirical best predictor that estimates the distribution of random effects from data, enhancing small area unemployment estimation.

Findings

Effective in small sample areas with limited data

Reduces computational burden compared to traditional methods

Provides accurate estimates with bias correction

Abstract

The Italian National Institute for Statistics regularly provides estimates of unemployment indicators using data from the Labor Force Survey. However, direct estimates of unemployment incidence cannot be released for Local Labor Market Areas. These are unplanned domains defined as clusters of municipalities; many are out-of-sample areas and the majority is characterized by a small sample size, which render direct estimates inadequate. The Empirical Best Predictor represents an appropriate, model-based, alternative. However, for non-Gaussian responses, its computation and the computation of the analytic approximation to its Mean Squared Error require the solution of (possibly) multiple integrals that, generally, have not a closed form. To solve the issue, Monte Carlo methods and parametric bootstrap are common choices, even though the computational burden is a non trivial task. In this paper, we propose a Semi-Parametric Empirical Best Predictor for a (possibly) non-linear mixed effect model by leaving the distribution of the area-specific random effects unspecified and estimating it from the observed data. This approach is known to lead to a discrete mixing distribution which helps avoid unverifiable parametric assumptions and heavy integral approximations. We also derive a second-order, bias-corrected, analytic approximation to the corresponding Mean Squared Error. Finite sample properties of the proposed approach are tested via a large scale simulation study. Furthermore, the proposal is applied to unit-level data from the 2012 Italian Labor Force Survey to estimate unemployment incidence for 611 Local Labor Market Areas using auxiliary information from administrative registers and the 2011 Census.

Figures

a, we may clearly observe that observed data lead to the estimation of a random effect with a clear degree of skewness. If the standard Gaussian assumption had been reasonable, the NPML estimate of the random effect distribution would have been a symmetric distribution centered around zero. As a consequence, we may conclude that such an assumption may not be that adequate for the current application. Furthermore, by comparing Figures 1b and 1c, we may observe that the parametric assumption also affects the posterior mean of the area-specific intercepts, leading to a less skewed distribution th

Tables

Table 1. Table 1: Sample percentage distribution with standard errors (S.E.) of Unemployed status by Sex-Age

	Unemployed	S.E.	Employed	S.E.	Inactive	S.E.
M: 15-24	11.0	0.4	21.5	0.5	67.6	0.6
M: 25-34	11.7	0.4	71.6	0.5	16.7	0.4
M: 35-65	5.3	0.1	70.8	0.3	23.8	0.2
F: 15-24	8.8	0.3	14.1	0.4	77.1	0.5
F: 25-34	11.7	0.4	53.9	0.6	34.5	0.5
F: 35-65	4.2	0.1	48.3	0.3	47.5	0.3

Table 2. Table 3: Distribution of the optimal number of mixture components across simulations.

	Scenario 1					Scenario 2
m / k	2	3	4	5		2	3	4	5
100	0.988	0.012	0.000	0.000		0.968	0.032	0.000	0.000
200	0.962	0.038	0.000	0.000		0.950	0.049	0.001	0.000
500	0.858	0.132	0.001	0.000		0.943	0.056	0.001	0.000

Table 3. Table 4: Average coverage rate over areas and computational time (in minutes) of RMSE ^ ​ ( p ^ i Naive ) ^ RMSE superscript subscript ^ 𝑝 𝑖 Naive \widehat{\mbox{RMSE}}(\hat{p}_{i}^{\mbox{\scriptsize Naive}}) , RMSE ^ ​ ( p ^ i EBP ) ^ RMSE superscript subscript ^ 𝑝 𝑖 EBP \widehat{\mbox{RMSE}}(\hat{p}_{i}^{\mbox{\scriptsize EBP}}) , RMSE ^ ​ ( p ^ i sp-EBP ) ^ RMSE superscript subscript ^ 𝑝 𝑖 sp-EBP \widehat{\mbox{RMSE}}(\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}}) , and RMSE ∗ ^ ​ ( p ^ i sp-EBP ) ^ superscript RMSE superscript subscript ^ 𝑝 𝑖 sp-EBP \widehat{\mbox{RMSE}^{*}}(\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}}) , for m = 100 , 200 , 500 𝑚 100 200 500 m=100,200,500 .

	Coverage					Computational Time
m	$100$	$200$	$500$			$100$	$200$	$500$
Scenario 1
sp-EBP	0.888	0.920	0.940			0.059	0.122	0.332
sp-EBP*	0.890	0.923	0.944			0.528	1.042	2.986
EBP	0.864	0.912				16.907	41.385	161.609*
Naive	0.962	0.967	0.976			0.005	0.018	0.206
Scenario 2
sp-EBP	0.928	0.931	0.927			0.053	0.110	0.303
sp-EBP*	0.933	0.933	0.928			0.481	0.959	2.437
EBP	0.966	0.974				17.229	42.112	162.528*
Naive	0.903	0.906	0.910			0.004	0.018	0.219
∗ Results refer to a single Monte Carlo draw

Table 4. Table 5: Parameter estimates, standard errors and corresponding p-values for the mixed logistic model fitted to the ILFS data based on an unspecific (left) and a Gaussian (right) random effect distribution.

	Unspecific			Gaussian
	Estimate	SE	p-value		Estimate	SE	p-value
Intercept	-3.052	0.176	0.001		-3.002	0.150	0.001
M:25-34	0.118	0.057	0.038		0.122	0.054	0.024
M:35-65	-0.787	0.054	0.001		-0.778	0.048	0.001
F:15-24	-0.222	0.062	0.001		-0.221	0.058	0.001
F:25-34	0.113	0.05	0.023		0.118	0.054	0.029
F:35-65	-1.039	0.052	0.001		-1.027	0.050	0.001
Middle School	0.206	0.059	0.001		0.211	0.054	0.001
High School	0.235	0.062	0.001		0.239	0.053	0.001
University Degree or Beyond	-0.032	0.083	0.682		-0.029	0.064	0.653
$\log ($U-count)	0.113	0.012	0.001		0.104	0.022	0.001
$\ell$	-21833.599				-21837.357
AIC	43695.199				43696.714

Equations

$$f_{y\mid\alpha}(y_{ij}\mid\mbox{\boldmath$\alpha$}_{i};\mbox{\boldmath$x$}_{ij})=\exp\left\{\frac{y_{ij}\theta_{ij}-b(\theta_{ij})}{a(\phi)}+c(y_{ij},\phi)\right\},$$

$$\theta_{ij}=\eta_{ij}=\mbox{\boldmath$x$}_{ij}^{\prime}\boldsymbol{\beta}+\mbox{\boldmath$w$}_{ij}^{\prime}\mbox{\boldmath$\alpha$}_{i}.$$

$$f_{y\mid\alpha}(\mbox{\boldmath$y$}_{i}\mid\mbox{\boldmath$\alpha$}_{i};\mbox{\boldmath$X$}_{i})=\prod_{j=1}^{N_{i}}f_{y\mid\alpha}(y_{ij}\mid\mbox{\boldmath$\alpha$}_{i};\mbox{\boldmath$x$}_{ij})=\exp\left\{\sum_{j=1}^{N_{i}}\frac{y_{ij}\,\theta_{ij}-b(\theta_{ij})}{a(\phi)}+c(y_{ij},\phi)\right\},$$

$$\displaystyle f_{y}(\mbox{\boldmath$y$}_{i};\mbox{\boldmath$X$}_{i})$$

$${\tilde{\zeta }}^{BP}$$

$$f_{\alpha\mid y}(\mbox{\boldmath$\alpha$}\mid\mbox{\boldmath$y$})=\frac{\prod_{i=1}^{m}f_{y\mid\alpha}(\mbox{\boldmath$y$}_{i}\mid\mbox{\boldmath$\alpha$}_{i};\mbox{\boldmath$X$}_{i})f_{\alpha}(\mbox{\boldmath$\alpha$}_{i})}{\prod_{i=1}^{m}\int_{\mathbb{R}^{q}}f_{y\mid\alpha}(\mbox{\boldmath$y$}_{i}\mid\mbox{\boldmath$\alpha$}_{i};\mbox{\boldmath$X$}_{i})f_{\alpha}(\mbox{\boldmath$\alpha$}_{i})d\mbox{\boldmath$\alpha$}_{i}},$$

$$\displaystyle L(\boldsymbol{\Phi})=\prod_{i=1}^{m}f_{y}(\mbox{\boldmath$y$}_{i};\mbox{\boldmath$X$}_{i})=\prod_{i=1}^{m}\int_{\mathbb{R}^{q}}f_{y\mid\alpha}(\mbox{\boldmath$y$}_{i}\mid\mbox{\boldmath$\alpha$}_{i};\mbox{\boldmath$X$}_{i})f_{\alpha}(\mbox{\boldmath$\alpha$}_{i})d\mbox{\boldmath$\alpha$}_{i},$$

$$\text{MSE}({\hat{\zeta }}^{EBP})$$

$$\displaystyle e(\mbox{\boldmath$\Phi$})$$

$$\displaystyle d(\mbox{\boldmath$\Phi$})$$

$$\displaystyle=\int_{\mathbb{R}^{m\times q}}\zeta(\mbox{\boldmath$\beta$},\mbox{\boldmath$\alpha$},\mbox{\boldmath$\Sigma$})^{2}\,f_{\alpha}(\mbox{\boldmath$\alpha$})\,d\mbox{\boldmath$\alpha$}-E_{y}\left[\left(\int_{\mathbb{R}^{m\times q}}\zeta(\mbox{\boldmath$\beta$},\mbox{\boldmath$\alpha$},\mbox{\boldmath$\Sigma$})f_{\alpha\mid y}(\mbox{\boldmath$\alpha$}\mid\mbox{\boldmath$y$})\,d\mbox{\boldmath$\alpha$}\right)^{2}\right],$$

$$\widehat{\text{MSE}}({\hat{\zeta }}^{EBP})=\frac{1}{m}e(\hat{\mathbf{\Phi }})+d(\hat{\mathbf{\Phi }}).$$

$$L(\boldsymbol{\Phi})=\prod_{i=1}^{m}\int_{\mathbb{R}^{q}}f_{y\mid\alpha}(\mbox{\boldmath$y$}_{i}\mid\mbox{\boldmath$\alpha$}_{i};\mbox{\boldmath$X$}_{i})\,d\mbox{\boldmath$\alpha$}_{i}\simeq\prod_{i=1}^{m}\sum_{g=1}^{G}f_{y\mid\alpha}(\mbox{\boldmath$y$}_{i}\mid\mbox{\boldmath$\xi$}_{g};\mbox{\boldmath$X$}_{i})\pi_{g},$$

$$\theta_{ijg}=\eta_{ijg}=\mbox{\boldmath$x$}_{ij}^{\prime}\boldsymbol{\beta}+\mbox{\boldmath$w$}_{ij}^{\prime}\mbox{\boldmath$\xi$}_{g}.$$

$${\tilde{\zeta }}^{\text{mbox}sp-BP}$$

$$\tau_{ig}=\Pr\left(z_{ig}=1\mid\mbox{\boldmath$y$}_{i}\right)=\frac{\pi_{g}\,f_{y\mid\alpha}(\mbox{\boldmath$y$}_{i}\mid\mbox{\boldmath$\xi$}_{g};\mbox{\boldmath$X$}_{i})}{\sum_{l=1}^{G}\pi_{l}\,f_{y\mid\alpha}(\mbox{\boldmath$y$}_{i}\mid\mbox{\boldmath$\xi$}_{l};\mbox{\boldmath$X$}_{i})},$$

$$\text{MSE}({\hat{\zeta }}^{\text{mbox}sp-EBP})$$

$$\displaystyle d^{\mbox{\scriptsize sp}}(\mbox{\boldmath$\Phi$})$$

$$\displaystyle=\sum_{g_{1}\cdots g_{m}}\zeta(\mbox{\boldmath$\beta$},\mbox{\boldmath$\xi$}_{g_{1},\dots,g_{m}},\mbox{\boldmath$\pi$})^{2}\,\prod_{i=1}^{m}\pi_{g_{i}}-E_{y}\left[\left(\sum_{g_{1}\cdots g_{m}}\zeta(\mbox{\boldmath$\beta$},\mbox{\boldmath$\xi$}_{g_{1},\dots,g_{m}},\mbox{\boldmath$\pi$})\prod_{i=1}^{m}\tau_{ig_{i}}\right)^{2}\right].$$

$$\widehat{\text{MSE}}({\hat{\zeta }}^{\text{mbox}sp-EBP})=\frac{1}{m}{e}^{\text{mbox}sp}(\hat{\mathbf{\Phi }})+d^{\text{mbox}sp}(\hat{\mathbf{\Phi }}).$$

$$d^{\mbox{\scriptsize sp}}(\hat{\mbox{\boldmath$\Phi$}})=d^{\mbox{\scriptsize sp}}(\mbox{\boldmath$\Phi$}_{0})+\left(\frac{\partial d^{\mbox{\scriptsize sp}}}{\partial\boldsymbol{\Phi}}\right)^{\prime}\Bigg{\rvert}_{\mbox{\boldmath$\Phi$}_{0}}(\hat{\boldsymbol{\Phi}}-\boldsymbol{\Phi}_{0})+\frac{1}{2}(\hat{\boldsymbol{\Phi}}-\boldsymbol{\Phi}_{0})^{\prime}\left(\frac{\partial^{2}d^{\mbox{\scriptsize sp}}}{\partial\boldsymbol{\Phi}\boldsymbol{\Phi}^{\prime}}\right)(\hat{\boldsymbol{\Phi}}-\boldsymbol{\Phi}_{0})+o_{p}(m^{-1}),$$

$$E[d^{\mbox{\scriptsize sp}}(\hat{\mbox{\boldmath$\Phi$}})]=d^{\mbox{\scriptsize sp}}({\mbox{\boldmath$\Phi$}})+\frac{1}{m}b^{\mbox{\scriptsize sp}}(\mbox{\boldmath$\Phi$})+o_{p}(m^{-1}),$$

$$b^{\text{mbox}sp}(\mathbf{\Phi })$$

$$\displaystyle=b_{1}^{\mbox{\scriptsize sp}}(\hat{\mbox{\boldmath$\Phi$}})+b_{2}^{\mbox{\scriptsize sp}}(\hat{\mbox{\boldmath$\Phi$}}).$$

$$\displaystyle b_{1}^{\mbox{\scriptsize sp}}(\hat{\mbox{\boldmath$\Phi$}})$$

$$\displaystyle-\left(\frac{\partial d^{\mbox{\scriptsize sp}}}{\partial\boldsymbol{\Phi}}\right)^{\prime}\Bigg{\rvert}_{\mbox{\boldmath$\Phi$}_{0}}\frac{m}{2}I_{e}({\boldsymbol{\Phi}_{0}})^{-1}tr\left\{I_{e}(\mbox{\boldmath$\Phi$}_{0})^{-1}\left[\frac{\partial I_{e}^{k}(\mbox{\boldmath$\Phi$}_{0})}{\partial\mbox{\boldmath$\Phi$}^{\prime}}\right]_{1\leq k\leq K}\right\}.$$

$$b_{2}^{\mbox{\scriptsize sp}}(\hat{\mbox{\boldmath$\Phi$}})=\frac{{m}}{2}tr\left\{\left(\frac{\partial^{2}d^{\mbox{\scriptsize sp}}}{\partial\boldsymbol{\Phi}\partial\boldsymbol{\Phi}^{\prime}}\right)\Bigg{\rvert}_{\mbox{\boldmath$\Phi$}_{0}}V(\hat{\mbox{\boldmath$\Phi$}})\right\}.$$

$$\widehat{\text{MSE}}^{*}(\hat{\zeta}^{\mbox{\scriptsize sp-EBP}})=d^{\mbox{\scriptsize sp}}(\hat{\boldsymbol{\Phi}})+\frac{1}{m}\left[e^{\mbox{\scriptsize sp}}(\hat{\boldsymbol{\Phi}})-b^{\mbox{\scriptsize sp}}(\hat{\mbox{\boldmath$\Phi$}})\right].$$

$$\theta_{ij}=\log\frac{p_{ij}}{1-p_{ij}}=\eta_{ij}=\alpha_{i}+\mbox{\boldmath$x$}_{ij}^{\prime}\boldsymbol{\beta}.$$

$${\bar{Y}}_i=\frac{1}{N_i}\sum _{j=1}^{N_i}{Y}_{ij},$$

$$p_i=\frac{1}{N_i}\sum _{j=1}^{N_i}{p}_{ij}.$$

$$\theta_{ijg}=\log\frac{p_{ijg}}{1-p_{ijg}}=\eta_{ijg}=\xi_{g}+\mbox{\boldmath$x$}_{ij}^{\prime}\boldsymbol{\beta}.$$

$$f_{ig}=f_{y\mid\alpha}(\mbox{\boldmath$y$}_{i}\mid\xi_{g};\mbox{\boldmath$X$}_{i})=\exp\left\{\sum_{j\in s_{i}}\left[y_{ij}\theta_{ijg}-\log\left(1+e^{\theta_{ijg}}\right)\right]\right\}.$$

Full text

Key Words: Binary data; Exponential Family; Finite Mixture; General parameters; Mixed logistic model; Unit-level model.

1 Introduction

The Italian National Institute for Statistics (ISTAT) regularly provides estimates of unemployment indicators based on data obtained through the Italian Labor Force Survey (ILFS). The ILFS allows to obtain quarterly estimates of the main aggregates regarding the labor market; these are particularly important both at the local and the central government levels for the development of labor market policies. These estimates are planned to be reliable at a given, chosen a priori, geographical level, and may not be suitable to all needs. For example, direct estimates of unemployment indicators cannot be disseminated for Local Labor Market Areas (LLMAs). These are 611 unplanned domains obtained as clusters of municipalities, defined at the Census on the basis of daily working commuting flows. In this context, direct survey estimates of unemployment incidence cannot be computed and/or published for most LLMAs. This is due to the presence of out-of-sample areas and to many LLMAs having a small sample size which leads to estimates with an unacceptable large coefficient of variation. For these reasons, ISTAT has implemented the use of indirect, model-based, small area estimators to produce official, yearly, estimates of unemployment incidence for Italian LLMAs (D’Aló et al.,, 2012, 2017).

Small Area Estimation (SAE) has received considerable attention in the past decades in terms of theoretical developments and applications to Official Statistics. An updated appraisal of available approaches for SAE is given in Rao and Molina, (2015). In this context, Generalized Liner Mixed Models (GLMMs, Laird and Ware,, 1982) represent a typical tool of analysis. Area-specific random effects are used to account for sources of unobserved heterogeneity that are not captured by the covariates and describe correlation between units within the same small area. For Gaussian data, Battese et al., (1988) introduced and Prasad and Rao, (1990) developed an Empirical Best Linear Unbiased Predictor (EBLUP) to estimate small area characteristics. Tailored to the purpose of the ILFS, D’Aló et al., (2017) developed unit-level linear mixed models with area- and time-specific random effects, which, based on data from different survey cycles, implement estimation using aggregate data to manage a large number of records. In fact, the ILFS is a continuous survey that collects, every year, information on almost $250,000$ households in $1,400$ municipalities for a total of $600,000$ individuals. However, many survey variables, such as the unemployment status, are categorical in nature and, therefore, SAE methods based on linear mixed models may not be fully appropriate.

Jiang and Lahiri, (2001) developed an Empirical Best Prediction (EBP) method for the area-specific random effects under a mixed logistic model providing a second-order, bias-corrected, estimator for the corresponding Mean Squared Error (MSE). Jiang, (2003) extended this approach to deal with GLMMs for general responses in the Exponential Family. Several functions of area-specific model effects are also investigated by the author. More recently, Boubeta et al., (2016, 2017) derived the EBP and the corresponding (second-order) approximation to the MSE under an area-level mixed Poisson model for small area counts, while Hobza and Morales, (2016) specifically focused on the development of an EBP for small area proportions under the unit-level mixed logistic model according to Jiang, (2003) and investigated the empirical behavior of the proposal through a large-scale simulation study. An extension of this latter approach to deal with longitudinal responses was also recently proposed by Hobza et al., (2018).

In all of these approaches, the area-specific random effects are assumed to be iid draws from a Gaussian distribution. One of the drawbacks associated with this assumption entails the computational burden required to derive parameter estimates, compute the EBP and, in particular, provide the corresponding measure of reliability. For non-Gaussian responses, we need to deal with (possibly) multiple integrals that do not admit a closed form expression and, therefore, need to be approximated. Numerical approaches, based e.g. on (adaptive) Gaussian quadrature or Laplace approximations (see e.g. Pinheiro and Bates,, 1995), or using Monte Carlo approximations (see e.g. McCulloch,, 1997) are frequently used for this purpose. For this reason, ad-hoc alternatives, mainly based on plug-in predictors and Taylor linearizations, were proposed and are currently largely applied (Saei and Chambers,, 2003; González-Manteiga et al.,, 2007; Molina et al.,, 2007; López-Vizcaíno et al.,, 2013).

In this paper, we describe a further alternative and develop a Semi-Parametric EBP (sp-EBP) for the small area parameters of interest and a second-order, bias-corrected, approximation to the corresponding MSE. In particular, we propose to leave the distribution of the area-specific random effects (the mixing distribution) unspecified and estimate it from the observed data via a NonParametric Maximum Likelihood approach (NPML - Simar,, 1976; Laird,, 1978; Lindsay, 1983a, ; Lindsay, 1983b, ). This estimate known to be a discrete distribution defined over a finite number of locations leading to a (semi-parametric) finite mixture model with a conditional kernel in the Exponential Family. The proposed approach offers a number of advantages. First, it allows us to avoid unverifiable assumptions on the random effect distribution; second, since mixture parameters are directly estimated from the data and are completely free to vary over the corresponding support, extreme and/or asymmetric departures from the homogeneous model can be easily accommodated. Last and more important, the discrete nature of the mixing distribution allows us to avoid integral approximations and considerably reduces the computational effort. The gain with respect to the parametric alternatives is particularly evident when analyzing non-Gaussian responses.

We present the proposed approach for a general small area parameter, starting from a general response with density in the Exponential Family and, later, focusing on the relevant case of binary data. We compare our proposal to the EBP (Jiang,, 2003) and to the plug-in estimator (e.g. Saei and Chambers,, 2003; González-Manteiga et al.,, 2007) in terms of prediction accuracy and computational burden in a large scale simulation study. Then, we prove the benefits from using the proposed sp-EBP approach on data from the ILFS to estimate unemployment incidence for the 611 LLMAs using auxiliary information from administrative registers and the 2011 Census. We compare the proposed approach with direct estimates, and with the two aforementioned approaches based on parametric mixed logistic models.

The paper is organized as follows. Section 2 presents the Italian Labor Force Survey, the estimation problem, and the auxiliary information available. Section 3 introduces the notation and a brief review of the EBP and its MSE approximation. In Section 4, we describe the proposed approach: section 4.1 entails maximum likelihood estimation, while the proposed sp-EBP and its MSE approximation are detailed in Section 4.2. Section 5 focuses on the case of binary responses. Section 6 reports the results of the simulation study, while Section 7 entails the application of the proposed approach to the ILFS data. Last, Section 8 summarizes our findings and provides guidelines for future research.

2 The ILFS data

The ILFS is the most important statistical source of information on the Italian labor market. The target population includes the members of all Italian households who regularly live within the national borders, have Italian or foreign citizenship, and are regularly enrolled in the municipal lists. Households registered as resident in Italy who habitually live abroad and permanent members of collective facilities (hospices, children’s homes, religious institutions, barracks, etc.) are excluded. A two-stage, municipality-household, sampling design is used to collect data. Primary sampling units are stratified by province (LAU1) and population size. Secondary sampling units are selected with equal probabilities. All individuals with usual residence in the dwelling are interviewed.

The ILFS provides quarterly estimates of the main aggregates for the labor market, such as employment status, type of work and work experience, by gender, age, and region (NUTS2). Here, we focus on data from the first quarter of 2012 which consist of measurements taken on $93,217$ units aged 15-65 and distributed in $453$ LLMAs. LLMAs refer to 611 unplanned domains obtained as clusters of municipalities where the bulk of the labour force lives and works, and where establishments can find the largest amount of the labour force necessary to occupy the offered jobs. They respond to the need for meaningfully, comparable, sub-regional labour market areas for the reporting and the analysis of statistics. LLMAs are defined on a functional basis, the key criterion being the proportion of commuters who cross the LLMA boundary on their way to work. In $2011$, with the last Census, LLMAs were re-defined by the analysis of daily working commuting flows using a new allocation process, an evolution of the previous algorithm. Nearly half of the LLMAs stands in the size class from $10,000$ up to $50,000$ inhabitants, whereas the highest proportion of the population (35.0$\%$) lives in LLMAs with a dimension between $100,000$ and $500,000$ inhabitants. In 332 LLMAs (over 70$\%$ of the national population), more than three quarter of the labour force lives and works in the same LLMA, that is self-containment is more than 75%.

Figure 1a shows the distribution of the LLMAs by sample size. This plot does not include the 158 areas with zero sample size. Among the observed LLMAs, the sample size ranges between $13$ (Acqui Terme, Piedmont Region) and $3,301$ (Milan, Lombardy Region). The mean value is equal to $205.8$, while quartiles are $61\,(25\%)$, $122\,(50\%)$, and $223\,(75\%)$, respectively. That is, several LLMAs are characterized by a very small sample size that hinders reliability of direct estimates. Figure 1b reports the distribution of the (percent) coefficient of variation (CV) for the direct estimates of unemployment incidence. The vast majority of estimates have a CV that is larger than $33\%$ that is usually considered as a threshold for reliability.

Our main interest is on the Employment Status variable which can take one out of three different categories: employed ($53.6\%$), unemployed ($6.6\%$), and inactive ($39.8\%$). Together with information on employment status for sampled individuals, the following explanatory variables are also available. Sex-Age: a categorical variable with six categories corresponding to female or male (F/M) and three age groups (15-24, 25-34, and 35-65); Educational Level: a categorical variable with four categories corresponding to no education or primary school diploma, secondary school diploma, high school diploma, and university degree or beyond; U-count: a discrete variable measuring the number of unemployed in a given sex-age group for each LLMA according to the $2011$ Census.

To have a first insight on the data, we report in Tables 2-2 the sample distribution of the Employment Status by Sex-Age and Educational Level, respectively. From these tables, we may observe that unemployment incidence is generally higher for people aged 25-34 in the sample, regardless of gender; for the other age groups, unemployment is less frequent among females. By looking at the last column of the table, we notice that females are more frequently inactive when compared to males, regardless of the age group. This is likely due to their engagement in housekeeping and explains why unemployment incidence is lower in this group. Similarly, by looking at Table 2, we may observe that the percentage of unemployment is relatively higher for individuals with higher education. Also in this case, by looking at the last column, it is evident that such a finding is mainly related to job hunting. For instance, $70\%$ of individuals with a primary school diploma or less are out of the job market, as they are not actively looking for a job, and this can be explained by a relatively older age.

As highlighted before, the prediction of unemployment incidence for the Italian LLMAs cannot be based on direct survey estimation as direct estimates cannot be computed and/or published for most of the LLMAs. For these reason, unit-level SAE methods may provide a viable tool to obtain such estimates. In the following section, we introduce the EBP approach by Jiang, (2003) for the estimation of small area parameters, together with the approach to approximate the corresponding MSE. As stated in Section 1, one of the main drawbacks of such a method is the computational complexity we have face with non-Gaussian data and a large number of observations/small areas, as for the ILFS data. In Section 4, we develop a computationally efficient alternative based on a semi-parametric approach.

3 The Empirical Best Prediction

Let $U$ denote a finite population of size $N$, which can be partitioned into $m$ non-overlapping small areas/domains, with $U_{i}$ denoting the $i$-th small area with size $N_{i},i=1,\dots,m$. For a given small area $i$, data consist of $N_{i}$ measurements of a response variable $Y_{ij}$ and a $p$-dimensional vector of covariates \mbox{\boldmathx}_{ij}=(x_{ij1},\dots,x_{ijp})^{\prime}, with $j=1,\dots,N_{i}$. Also, let \mbox{\boldmath\alpha}_{1},\dots,\mbox{\boldmath\alpha}_{m} be iid, $q$-dimensional, vectors of area-specific random effects ($q\leq p$) with density $f_{\alpha}(\cdot)$, E_{\alpha}(\mbox{\boldmath\alpha}_{i})=0, and E_{\alpha}(\mbox{\boldmath\alpha}_{i}\mbox{\boldmath\alpha}_{i}^{\prime})=\mbox{\boldmath\Sigma} for all $i=1,\dots,m$. Last, let \mbox{\boldmathw}_{ij} denote a $q$-dimensional subset of \mbox{\boldmathx}_{ij} associated to \mbox{\boldmath\alpha}_{i}. We assume that a sample of size $n$ is drawn from the above population and denote by $s_{i}$ the set containing the $n_{i}$ population indexes of sample units belonging to small area $i$, with $n=\sum_{i=1}^{m}n_{i}$. On the other hand, the set $r_{i}\subseteq U_{i}$ contains the $N_{i}-n_{i}$ indexes for non-sampled units in small area $i$. For ease of notation, we assume that all areas are sampled, even though the presence of out of sample areas can be easily accommodated. We further assume that values of $Y_{ij}$ are known only for the sample ($i=1,\dots,m,j\in s_{i})$, while the values of \mbox{\boldmathx}_{ij} and \mbox{\boldmathw}_{ij}, are known for all units in the population ($i=1,\dots,m,j=1,\dots,N_{i}$). This assumption can be quite restrictive in some real-world applications, since it implies the availability of individual population information. However, when the auxiliary variables are categorical and/or take a finite number of values, the assumption can be relaxed. We will discuss this issue in more details in the application to the ILFS data. Last, we assume that sampling is non-informative for the small area distribution of Y_{ij}\mid\mbox{\boldmathx}_{ij}, allowing us to use population level models with sample data.

3.1 The model

According to a local independence assumption, we assume that, conditional on the area-specific random effects \mbox{\boldmath\alpha}_{i}, responses $Y_{ij}$ from the same small area $i$ are independent with density in the Exponential Family

[TABLE]

for $i=1,\dots,m$ and $j=1,\dots,N_{i}$. In the previous expression, $\phi$ is a dispersion parameter, $a(\cdot),b(\cdot)$ and $c(\cdot)$ are known functions and $\theta_{ij}$ is the canonical parameter for the chosen member of the family. Let $\beta$ denote a $p$-dimensional vector of fixed regression coefficients and let us assume that $\theta_{ij}$ is modeled via the following regression model

[TABLE]

The joint distribution of \mbox{\boldmathy}_{i}=(y_{i1},\dots,y_{iN_{i}})^{\prime} for the $i$-th small area, conditional on the vector of area-specific random effects \mbox{\boldmath\alpha}_{i}, is obtained by exploiting conditional independence

[TABLE]

where \mbox{\boldmathX}_{i} denotes the matrix of covariates associated to units in the $i$-th area. The marginal distribution of the area-specific sequence \mbox{\boldmathy}_{i} is obtained by integrating out \mbox{\boldmath\alpha}_{i}:

[TABLE]

Typically, a parametric specification for f_{\alpha}(\mbox{\boldmath\alpha}_{i}) is adopted, with a common choice being the zero mean, multivariate, Gaussian distribution. It is worth noticing that an implicit exogeneity assumption of observed covariates \mbox{\boldmathx}_{ij} is taken, that is f_{\alpha}(\mbox{\boldmath\alpha}_{i}\mid\mbox{\boldmathX}_{i})=f_{\alpha}(\mbox{\boldmath\alpha}_{i}) or E(\mbox{\boldmath\alpha}_{i}\mid\mbox{\boldmathX}_{i})=E(\mbox{\boldmath\alpha}_{i})=\mbox{\boldmath0}. When this assumption is not fulfilled, the auxiliary regression approach by Mundlak, (1978) can be adopted. This slightly modifies the linear predictor above and produces area-specific random effects that are (linearly) free of \mbox{\boldmathX}_{i}, see Neuhaus and McCulloch, (2006). In the following, we will assume that, if needed, such an approach is applied and that f_{\alpha}(\mbox{\boldmath\alpha}_{i}\mid\mbox{\boldmathX}_{i})=f_{\alpha}(\mbox{\boldmath\alpha}_{i}).

3.2 EBP and MSE approximation

We are interested in using sample data on responses $Y_{ij}$ ($i=1,\dots,m,j\in s_{i}$) and population data on covariates \mbox{\boldmathx}_{ij} ($i=1,\dots,m,j=1,\dots,N_{i}$) to predict a (possibly) non-linear function of fixed and random effects, say \zeta(\mbox{\boldmath\beta},\mbox{\boldmath\alpha},\mbox{\boldmath\Sigma}), with {\mbox{\boldmath\alpha}=(\mbox{\boldmath\alpha}_{1},\dots,\mbox{\boldmath\alpha}_{m})}. According to Jiang, (2003), the Best Predictor (BP) of $\zeta$ in terms of minimum MSE is given by

[TABLE]

where

[TABLE]

and f_{y\mid\alpha}(\mbox{\boldmathy}_{i}\mid\mbox{\boldmath\alpha}_{i};\mbox{\boldmathX}_{i})=\prod_{j\in s_{i}}f_{y\mid\alpha}(y_{ij}\mid\mbox{\boldmath\alpha}_{i};\mbox{\boldmathx}_{ij}). Since model parameters \boldsymbol{\Phi}=(\boldsymbol{\beta},\phi,\mbox{\boldmath\Sigma}) are unknown, they need to be estimated. Estimation can be accomplished by maximizing the observed data likelihood function:

[TABLE]

where, as before, f_{y\mid\alpha}(\mbox{\boldmathy}_{i}\mid\mbox{\boldmath\alpha}_{i};\mbox{\boldmathX}_{i}) refers to sample data only. To maximize equation (2), we need to evaluate an integral defined over the support of the area-specific random effects and this can be directly done in few cases, for instance when $f_{y\mid\alpha}(\cdot\mid\cdot)$ and $f_{\alpha}(\cdot)$ are conjugate. In all other cases, numerical approximations (e.g. Gaussian quadrature techniques) or simulation based methods (e.g. Monte Carlo integration) need to be used, often leading to a non trivial computational complexity. To overcome the issue, Jiang, (1998) suggested to derive estimates by exploiting the method of moments. A Penalized Quasi Likelihood (PQL) approach (e.g., Breslow and Clayton,, 1993) represents a further alternative which is less computationally demanding, even though it may provide inconsistent model parameter estimates (see e.g. Rodriguez and Goldman,, 1995).

Once parameters are estimated, we may compute the EBP of $\zeta$, that is \hat{\zeta}^{EBP}=\tilde{\zeta}^{BP}(\mbox{\hat{\mbox{\boldmath$\beta$}}},\mbox{\hat{\mbox{\boldmath$\alpha$}}},\hat{\mbox{\boldmath\Sigma}}). To evaluate the quality of such predictions, the second-order MSE estimator detailed by Jiang, (2003) can be considered. Under mild regularity conditions, the following decomposition holds:

[TABLE]

where

[TABLE]

and f_{\alpha}(\mbox{\boldmath\alpha})=\prod_{i=1}^{m}f_{\alpha}(\mbox{\boldmath\alpha}_{i}) denotes the joint density of the random effects \mbox{\boldmath\alpha}_{i},i=1,\dots,m. An estimator of $\text{MSE}(\hat{\zeta}^{EBP})$ can be obtained by replacing $\boldsymbol{\Phi}$ in equation (3) by a consistent estimator, that is

[TABLE]

However, as outlined by Jiang, (2003), while we get an error of order $o_{p}(m^{-1})$ when we replace $\Phi$ by \hat{\mbox{\boldmath\Phi}} into $e(\boldsymbol{\Phi})$, a bias correction is needed to obtain an unbiased estimator for $d({\boldsymbol{\Phi}})$. We discuss this issue in more detail in the following.

As it is clear, computing the MSE requires the solution of (multiple) integrals that may not admit a closed form expression. As stated before, Monte Carlo approximations or numerical integration techniques are required and this makes the computation extremely time-consuming. Bootstrap may represent a further alternative, particularly when dealing with a limited number of small areas. However, when $m$ is large, as in the case of the ILFS data, neither the analytic MSE approximation nor the bootstrap represent viable strategies due to computational issues. González-Manteiga et al., (2007) proposed a non-optimal Prasad-Rao-type MSE estimator derived from a Taylor series approximation. This estimator fails when sample sizes are too small, while its behavior is proved to be reliable in the case of large sample sizes.

4 The Semi-Parametric Empirical Best Prediction

As highlighted before, deriving the EBP of small area parameters and the corresponding MSE approximation as detailed by Jiang, (2003) is a non trivial task. In this section, we develop a computationally convenient alternative that allows us to avoid unverifiable parametric assumption on the random effect distribution. In Section 4.1, we present the proposed approach to derive model parameter estimates within a maximum likelihood framework. In Section 4.2, we detail the proposed Semi-Parametric Empirical Best Predictor (sp-EBP) and the corresponding second-order, bias-corrected, MSE estimator.

4.1 Model parameter estimation

When dealing with non-Gaussian responses and GLMMs with Gaussian random effects, maximum likelihood (ML) estimators, although optimal, can be time consuming as we need to approximate (possibly multi-dimensional) integrals that do not admit a closed form expression. An alternative may be based on leaving the distribution of \mbox{\boldmath\alpha}_{i} completely unspecified and follow the approach detailed by Aitkin, (1996, 1999). The area-specific random effects are treated as nuisance parameters and a NonParametric Maximum Likelihood (NPML) estimate of their distribution is derived. Different contributions to the theory of NPML can be found in the literature (Simar,, 1976; Laird,, 1978; Böhning,, 1982; Lindsay, 1983a, ; Lindsay, 1983b, ). Results by Lindsay, 1983a ; Lindsay, 1983b show that, as long as the (log-) likelihood function is bounded, it is maximized by a discrete distribution defined on, at most, as many support points as the number of distinct area profiles in the sample. In particular, the mixing distribution estimate is a discrete distribution which puts masses $\pi_{g}>0$ on locations \mbox{\boldmath\xi}_{g}=(\xi_{g1},\dots,\xi_{gq})^{\prime}, $g=1,\dots,G$, where the constraint $\sum_{g=1}^{G}\pi_{g}=1$ holds. In a regression context, the number of locations $G$ is bounded from above by the number of different profiles (\mbox{\boldmathy}_{i},{\mbox{\boldmathX}}_{i}) in the sample. That is, in the presence of categorical covariates, the number of locations does not necessarily grow with $m$.

Let $\Phi$ denote the global vector of model parameters, \boldsymbol{\Phi}=(\boldsymbol{\beta},\phi,\mbox{\boldmath\xi}_{1},\dots,\mbox{\boldmath\xi}_{G},\pi_{1},\dots,\pi_{G})^{\prime}; the observed data likelihood is approximated by

[TABLE]

where f_{y\mid\alpha}(\mbox{\boldmathy}_{i}\mid\mbox{\boldmath\xi}_{g};\mbox{\boldmathX}_{i})=\prod_{j\in s_{i}}f_{y\mid\alpha}(y_{ij}\mid\mbox{\boldmath\alpha}_{i}=\mbox{\boldmath\xi}_{g};\mbox{\boldmathx}_{ij}) denotes the product of densities in the Exponential Family with canonical parameter $\theta_{ijg}$ defined by the following (mixed) model:

[TABLE]

As it is clear, expression (6) resembles the likelihood of a finite mixture of distributions, with weights \pi_{g}=\text{Pr}(\mbox{\boldmath\alpha}_{i}=\mbox{\boldmath\xi}_{g}). That is, \mbox{\boldmath\alpha}_{i}\sim\sum_{g=1}^{G}\pi_{g}\delta({\mbox{\boldmath\xi}_{g}}), where $\delta(a)$ is a one-point distribution putting a unit mass at $a$. It is worth noticing that, while the discrete nature of the estimate for $f_{\alpha}(\cdot)$ may seem unappealing, most approximation techniques (e.g. based on Gaussian quadrature or Monte Carlo approaches) applied when a parametric specification is considered, are exactly of the type in equation (6). The only substantial difference is that locations \mbox{\boldmath\xi}_{g} and masses $\pi_{g}$ in the present proposal are estimated to best fit observed data.

To maximize the likelihood in (6), the EM algorithm (Dempster et al.,, 1977) can be employed. A drawback of such an algorithm is that it does not directly provide estimates for the covariance matrix of model parameters. A frequent solution to this issue is based on the use of the Oakes’ formula (Oakes,, 1999), as detailed in Sections 1 and 2 of the on-line Supplementary Material, where the EM algorithm is described. A crucial point in the proposed approach is the choice of the number of mixture components in (6). A simple and frequently used solution is as follows: parameter estimates are computed for varying values of $G$ and the model with the best fit, typically measured by penalized likelihood criteria (such as AIC or BIC), is retained. Typically, the optimal $G$ increases either (i) when the variability of the random effect distribution increases or (ii) when the number of small areas increases as, in this case, this may lead to a higher number of distinct area profiles in the sample. As long as convergence is entailed, the order for the mixing distribution estimate is $O_{p}(m^{-1/4})$, as compared to $O_{p}(m^{-1/2})$ for ML parameter estimates in regular models (see Chen,, 1995). However, according to Lindsay and Lesperance, (1995), some smooth functionals, such as the empirical Bayes estimates, can be estimated at the usual $O_{p}(m^{-1/2})$ rate. Furthermore, as shown by Redner and Walker, (1984), when the order of the mixture is finite and known, that is when \mbox{\boldmath\alpha}_{i}\sim\sum_{g=1}^{G}\pi_{g}\delta(\mbox{\boldmath\xi}_{g}) is the true mixing, with $G$ known, the usual ML asymptotics apply.

4.2 Semi-Parametric EBP and MSE approximation

Let us now turn to the main problem of interest, where we have a finite population of size $N$ which can be partitioned into $m$ non-overlapping domains or small areas. Furthermore, let \mbox{\boldmath\xi}=(\mbox{\boldmath\xi}_{1},\ldots,\mbox{\boldmath\xi}_{G})^{\prime} and \mbox{\boldmath\pi}=(\pi_{1},\ldots,\pi_{G})^{\prime} denote the vectors of locations and masses of the finite mixture, respectively. We aim at predicting a (possibly) non-linear function of fixed and random effects, \zeta(\mbox{\boldmath\beta},\mbox{\boldmath\xi},\mbox{\boldmath\pi}) by exploiting sample data on responses $Y_{ij}$ and populations data on covariates \mbox{\boldmathx}_{ij}. Under the proposed approach, the Semi-Parametric Best Predictor (sp-BP) of $\zeta$ is defined according to the following expression:

[TABLE]

where $\sum_{g_{1}\dots g_{m}}$ is a shorthand for $\sum_{g_{1}=1}^{G}\cdots\sum_{g_{m}=1}^{G}$, \mbox{\boldmath\xi}_{g_{1},\dots,g_{m}}=(\mbox{\boldmath\xi}_{g_{1}},\dots,\mbox{\boldmath\xi}_{g_{m}})^{\prime}, and $\tau_{ig}$ denotes the posterior probability for the $i$-th small area to belong to the $g$-th component of the finite mixture. In particular, denoting by $z_{ig},i=1,\dots,m,g=1,\dots,G,$ the component membership indicator for the $i$-th small area, $\tau_{ig}$ is defined by

[TABLE]

where, as before, f_{y\mid\alpha}(\mbox{\boldmathy}_{i}\mid\mbox{\boldmath\xi}_{g};\mbox{\boldmathX}_{i}) refers to sample data only. As it is clear, expression (7) denotes the expected value of \zeta(\mbox{\boldmath\beta},\mbox{\boldmath\xi},\mbox{\boldmath\pi}), with respect to the posterior distribution of the random effects $\alpha$. Since this is a discrete distribution, the integral approximation which is required in equation (1) directly translates into simpler summations.

An estimate of $\tilde{\zeta}^{\mbox{\scriptsize sp-BP}}$ can be obtained by replacing model parameters \mbox{\boldmath\beta},\mbox{\boldmath\xi}, and $\pi$ by consistent estimates. Here, we consider the estimates derived by the EM algorithm described in Section 1 of the on-line Supplementary Material. In the following, we will refer to such a quantity as the Semi-Parametric Empirical Best Predictor (sp-EBP) of $\zeta$, denoted by \hat{\zeta}^{\mbox{\scriptsize sp-EBP}}=\tilde{\zeta}^{\mbox{\scriptsize sp-BP}}(\mbox{\hat{\mbox{\boldmath$\beta$}}},\mbox{\hat{\mbox{\boldmath$\xi$}}},\mbox{\hat{\mbox{\boldmath$\pi$}}}).

To evaluate the quality of predictions, we develop an analytic approximation to the MSE of $\hat{\zeta}^{\mbox{\scriptsize sp-EBP}}$ based on the approach by Jiang, (2003), but considering a maximum likelihood estimator. Starting from equation (3), the MSE of the sp-EBP is given by:

[TABLE]

where the former term, e^{\mbox{\scriptsize sp}}(\mbox{\boldmath\Phi}), is defined according to expression (4) and can be derived by computing model derivatives with respect to $\beta$, $\alpha$, and $\pi$, together with the covariance matrix of model parameter estimates, V(\hat{\mbox{\boldmath\Phi}}). See Section 2 in the on-line Supplementary Material for computational details. On the other hand, d^{\mbox{\scriptsize sp}}(\mbox{\boldmath\Phi}) can be derived as follows

[TABLE]

The computational burden to obtain the above quantities is substantially lower than that required for the approach by Jiang, (2003). Intractable integrals appearing in equations (4) and (5) all translate into simple summations which can be solved analytically.

An estimator of $\text{MSE}(\hat{\zeta}^{\mbox{\scriptsize sp-EBP}})$ is obtained by replacing $\boldsymbol{\Phi}$ in (9) by a consistent estimator such as that obtained by maximizing the observed data likelihood in equation (6). That is,

[TABLE]

However, as we remarked before, this approach does not directly lead to an unbiased estimator of $\text{MSE}(\hat{\zeta}^{\mbox{\scriptsize sp-EBP}})$. When replacing \hat{\mbox{\boldmath\Phi}} in d^{\mbox{\scriptsize sp}}(\mbox{\boldmath\Phi}), we get an error of order $O_{p}(m^{-1/2})$ and a bias correction term needs to be considered. Jiang, (2003) provided an explicit expression for such a term when model parameters are estimated by the method of moments. Clearly, under the current approach, these results do directly not hold but, rather, need to be adapted.

Let \mbox{\boldmath\Phi}_{0} denote the “true” vector of model parameters and let us consider a second-order Taylor expansion of d^{\mbox{\scriptsize sp}}({\mbox{\boldmath\Phi}}) around \mbox{\boldmath\Phi}_{0} evaluated at \hat{\mbox{\boldmath\Phi}}:

[TABLE]

where $d^{\mbox{\scriptsize sp}}$ is a shorthand for d^{\mbox{\scriptsize sp}}(\mbox{\boldmath\Phi}). From expression (11), it is easy to see that

[TABLE]

where b^{\mbox{\scriptsize sp}}(\mbox{\boldmath\Phi}) denotes a bias correction defined as

[TABLE]

As it is shown in Section 3 of the on-line Supplementary Material, the former term on the right hand side of equation (4.2) is given by

[TABLE]

Here, I_{e}(\mbox{\boldmath\Phi}_{0}) denotes the expected information matrix, while S^{k}(\mbox{\boldmath\Phi}) and I_{e}^{k}(\mbox{\boldmath\Phi}_{0}) denote the $k$-th element of the score function S(\mbox{\boldmath\Phi}) and the $k$-th row of the expected information matrix I_{e}(\mbox{\boldmath\Phi}_{0}), respectively.

On the other hand, it can be shown that the latter term in equation (4.2) , b_{2}^{\mbox{\scriptsize sp}}(\hat{\mbox{\boldmath\Phi}}), can be computed as

[TABLE]

A second order, bias corrected, estimator of $\text{MSE}(\hat{\zeta}^{\mbox{\scriptsize sp-EBP}})$ is then given by

[TABLE]

We report the computational details required to derive b_{1}^{\mbox{\scriptsize sp}}(\hat{\mbox{\boldmath\Phi}}) and b_{2}^{\mbox{\scriptsize sp}}(\hat{\mbox{\boldmath\Phi}}) in Section 3 of the on-line Supplementary Material.

5 A special case: binary data

In this section, we focus on the relevant case of binary responses modeled via a mixed logistic model with random intercepts. Let $Y_{ij}$ denote the binary response associated to unit $j$ in the $i$-th small area ($i=1,\dots,m,j=1,\dots,N_{i}$), and let $\alpha_{i}$ denote an area-specific random effect. Again, let \mbox{\boldmathx}_{ij} denote a $p$-dimensional vector of covariates, and \mbox{\boldmathX}_{i} the matrix of covariates for the $i$-th small area. We assume that, conditional on $\alpha_{i}$, responses for units in the $i$-th small area are independent Bernoulli random variables with success probability $p_{ij}$, described by the following mixed logistic model:

[TABLE]

In the equation above, $\boldsymbol{\beta}$ is a $p$-dimensional vector of fixed model parameters that describes the effect of observed covariates on the logit transform of $p_{ij}$. We consider the practical problem of predicting small area proportions

[TABLE]

using the GLMM in equation (14). To this end, we will use the EBP for the quantity

[TABLE]

In fact, since $N_{i}$ is usually very large in most applications, as it is the case in the one at hand, the EBP for $p_{i}$ can also be used to predict the indicator $\bar{Y}_{i}$. Let us assume that responses $Y_{ij}$ are observed for sampled units only ($i=1,\dots,m,j\in s_{i}$), while covariates \mbox{\boldmathx}_{ij} are available at the population level ($i=1,\dots,m,j=1,\dots,N_{i}$). Following the approach detailed in the previous sections, we leave the distribution of the area-specific random effects in equation (14) unspecified and approximate it via a discrete distribution that puts masses $\pi_{g}>0$ on locations $\xi_{1},\dots,\xi_{G}$, with $\sum_{g=1}^{G}\pi_{g}=1$. By adopting a canonical link function, the logistic transform of the success probability for a generic area $i$ in the $g$-th component of the finite mixture is given by

[TABLE]

Using the standard notation for the Exponential Family, the joint conditional density for the observed responses in the $i$-th small area and the $g$-th component is

[TABLE]

Turning back to the problem of estimating $p_{i}$ in equation (15), the corresponding sp-BP is given by

[TABLE]

where $y_{i\cdot}=\sum_{j\in s_{i}}y_{ij}$ and $p_{ig}=N_{i}^{-1}\sum_{j=1}^{N_{i}}p_{ijg}$.

By letting

[TABLE]

the sp-BP of $p_{i}$ is given by

[TABLE]

The corresponding sp-EBP, denoted by $\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}}$, is obtained by substituting ML estimates of model parameters into expression (16):

[TABLE]

while the quality of predictions obtained via $\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}}$ can be evaluated through the following MSE expression:

[TABLE]

where

[TABLE]

and

[TABLE]

Here, $\tilde{p}_{i(h)}^{\mbox{\scriptsize sp-BP}}$ denotes the sp-BP of $p_{i}$ conditional on $y_{i.}=h$, that is

[TABLE]

The term \Pr\left(Y_{i\cdot}=h;\mbox{\boldmathX}_{i}\right) is obtained as

[TABLE]

where \Pr\left(Y_{i\cdot}=h\mid\xi_{g};\mbox{\boldmathX}_{i}\right) represents the probability of observing $h$ successes in $n_{i}$ independent, but non identically distributed, Bernoulli trials. This quantity can be obtained using the probability mass function of a Poisson-Binomial random variable (see Chen and Liu,, 1997) with parameter $(p_{i1g},\dots,p_{in_{i}g})$. The last term in equation (18) is obtained as

[TABLE]

where V(\hat{\mbox{\boldmath\Phi}}) is the covariance matrix of model parameter estimates and ${\partial\tilde{p}_{i(h)}^{\mbox{\scriptsize sp-BP}}}/{\partial\boldsymbol{\Phi}}$ is the vector of model derivatives conditional on $y_{i\cdot}=h$. Explicit formulas for these latter quantities are provided in Section 4 of the on-line Supplementary Material.

The second-order, bias-corrected, estimator of $\text{MSE}(\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}})$, that is $\widehat{\mbox{MSE}}^{*}(\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}})$, is obtained according to expression (13), after adapting the bias correction term to the binary case.

6 Model-based simulation study

In this section, we evaluate the empirical properties of the proposed approach via a large scale (model-based) simulation study. This consists of $T=1,000$ samples, where binary population data are generated under some model assumptions and sample data are selected from the simulated population. In particular, population data are generated considering $m=100,200,500$ small areas; then, samples are selected by simple random sampling without replacement within each area. The population and the sample sizes are constant across areas and are fixed to $N_{i}=100$ and $n_{i}=10$, respectively. According to the simulation study discussed by González-Manteiga et al., (2007), for each unit $j$ in small area $i$, we generate the target variable $Y_{ij},i=1,\dots,m,j=1,\dots,N_{i}$, from a Bernoulli distribution with success probability defined by

[TABLE]

with $\beta=1$, $x_{ij}\sim\mbox{Unif}(-1,b_{i})$, and $b_{i}=i/8,\leavevmode\nobreak\ i/16,\leavevmode\nobreak\ i/48$ for $m=100,200$ and $500$, respectively. To evaluate the impact of parametric assumptions on the distribution of the area-specific random effects, we considered two different scenarios. The first one (Scenario 1) uses area-specific random effects from a zero mean, Gaussian, distribution with standard deviation equal to $\sigma_{1}=0.5$. The second scenario (Scenario 2) involves area-specific random effects generated from a mixture of Gaussian distributions, $\alpha_{i}\sim\nu N(\mu_{1},\sigma_{2})+(1-\nu)N(\mu_{2},\sigma_{2})$, where $\nu$ represents a random draw from a Bernoulli distribution $\Pr(\nu=1)=0.7$, $\mu_{1}=0$, $\mu_{2}=3$, and $\sigma_{2}=0.05$. Based on this latter quantity, it is evident that, under this scenario, the random effect distribution closely resembles that of a discrete distribution putting masses $\nu$ and $1-\nu$ on locations $\mu_{1}$ and $\mu_{2}$. In this framework, the population is made by two separate sets of small areas having different baseline levels for the success probabilities. This may be reasonable e.g. for properly representing non-homogeneous unemployment rates typically observed in the North/South of Italy, as we will see in Section 7. Clearly, the chosen scenarios represent two extreme situations; we expect that, in real data applications, the random effect distribution lies in between them.

In this simulation study, our aim is that of evaluating the empirical behavior of the proposed approach. For each simulated sample, we estimated model parameters for a varying number of mixture components ($G=2,\dots,5$) and selected the optimal $G$ according to the AIC index. We report in Table 3 the distribution of the optimal number of mixture components $G$ across simulations. As it can be observed, in most of the cases the AIC index leads to selecting a model with $G=2$ components only. This reflects the reduced variability of the random effect distribution considered under both simulation scenarios. However, it is worth to highlight that, for higher sample sizes, the chance of selecting a higher $G$ slightly increases, especially when $\alpha_{i}$ is a random draw from a Gaussian density. This result is clearly related to the requirement of a higher number of components to properly approximate the “true”, continuous, distribution of the area-specific effects.

Starting from parameter estimates derived from the proposed approach, the sp-EBP for small area proportions was derived according to equation (17). The proposed predictor was then compared with the parametric EBP by Jiang and Lahiri, (2001) and the Naive predictor considered in González-Manteiga et al., (2007), both based on the assumption of Gaussian random effects. For the EBP, parameter estimates were derived via the ML approach based on a Laplace approximation available in the glmer function from the R lme4 package (Bates et al.,, 2015). Given the estimates, small area proportions and corresponding MSEs were derived by adopting the formulas detailed in Section 3. To evaluate the intractable integrals, we followed the approach suggested by Boubeta et al., (2016). That is, we started by generating $B=2,500$ replicates of the area-specific random effects $\alpha_{i}^{(b)}$ from a Gaussian density with zero mean and variance equal to the corresponding ML estimate. Then, we considered their antithetic transform $\alpha_{i}^{(B+b)}=-\alpha_{i}^{(b)}$ to obtain $2B$ random effect values. Finally, integrals were approximated by the corresponding empirical means. In the following, we will denote EBP estimates of small area proportions by $\hat{p}_{i}^{\mbox{\scriptsize EBP}}$. Although this approach is optimal, the computational complexity greatly limits its applicability. Via the current simulation study, we aim at understanding whether the sp-EBP approach we propose could represent an effective alternative, which is optimal in terms of minimum MSE and simpler from a computational point of view.

For completeness, we also included in the simulation study results from the Naive approach. In this case, parameter estimates were obtained using the PQL approach via the glmmPQL function from the R MASS package (Venables and Ripley,, 2002). To get predictions, parameter estimates were directly plugged into the expression for the area-specific proportions:

[TABLE]

The performance of the small area estimators were evaluated by computing, for each area $i=1,\dots,m,$ the bias and the Root Mean Squared Error (RMSE), defined as follows:

[TABLE]

[TABLE]

where $\hat{p}_{it}^{\mbox{\scriptsize Model}}$ denotes the model-based proportion estimate for the $i$-th small area in the $t$-th simulated sample obtained via either the EBP ($\hat{p}_{i}^{\mbox{\scriptsize EBP}}$), the sp-EBP ($\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}}$), or the Naive ($\hat{p}_{i}^{\mbox{\scriptsize Naive}}$) approach. For completeness, we also report in Section 5 of the on-line Supplementary Material the distribution of the Mean Absoute Error (MAE) across small areas for the EBP, the sp-EBP, and the Naive predictor under different experimental scenarios. Together with the bias, MAE is frequently used to evaluate the quality of predictions in the small area literature, even though it may not be a consistent tool to evaluate predictions obtained by posterior means (Gneiting,, 2011).

Figures 2 and 3 show the BIAS and the RMSE distribution across small areas for the three estimators under investigation for Scenario 1 and $m=100,200,500$, respectively; the red line denotes the corresponding mean values. As expected, when looking at the first two panels (i.e. $m=100,\leavevmode\nobreak\ 200$), the sp-EBP performs better than the Naive estimator and slightly worse than the EBP, with a gap that reduces as $m$ increases, both in terms of BIAS and RMSE. When $m=500$, performance values of EBP are not showed due to the computational burden required to get the estimates: for one replication, we needed 161.612 minutes on an Intel(R) I5-3330 architecture – 3.0 GHz, and, therefore, we couldn’t obtain results for $T=1,000$ replications in a reasonable amount of time.

Figures 4 and 5 show the performance of the estimators under Scenario 2. As before, results for the EBP approach for $m=500$ are not showed due to computational issues. As it is clear by looking at these plots, when the assumption of Gaussian random effects does not hold, parametric approaches seem to produce predictions with a reduced quality than those obtained via the semi-parametric alternative we propose. In particular, we may notice that $\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}}$ clearly outperforms the two competitors in terms of both bias and RMSE. Also, results for $\hat{p}_{i}^{\mbox{\scriptsize Naive}}$ and $\hat{p}_{i}^{\mbox{\scriptsize EBP}}$ seem to slightly worsen as $m$ increases. This may be possibly due to the higher information available and the stronger impact of the random effect distribution on the overall response variability when the number of small areas increases.

A further purpose of this simulation study is to investigate the performance of the MSE estimators to evaluate the accuracy of the predictions we discussed so far. In particular, for $\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}}$, we considered the proposed MSE estimator reported in equations (10) and (13); we will refer to the square root of these quantities as $\widehat{\mbox{RMSE}}(\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}})$ and $\widehat{\mbox{RMSE}^{*}}(\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}})$, respectively. For the estimator $\hat{p}_{i}^{\mbox{\scriptsize EBP}}$, we used the approximate MSE estimator proposed by Hobza and Morales, (2016); the corresponding square root will be denoted by $\widehat{\mbox{RMSE}}(\hat{p}_{i}^{\mbox{\scriptsize EBP}})$. Last, for the Naive predictor $\hat{p}_{i}^{\mbox{\scriptsize Naive}}$, we considered the approach suggested by González-Manteiga et al., (2007), based on linearizing the GLMM in equation (20) and, then, applying the Prasad-Rao MSE approximation for the corresponding linear mixed model; the square root of such an estimator will be denoted by $\widehat{\mbox{RMSE}}(\hat{p}_{i}^{\mbox{\scriptsize Naive}})$.

The performance of the RMSE estimators were evaluated by considering the ratio (R) between the estimated RMSE for the model-based estimates and the corresponding actual RMSE for each small area prediction, that is:

[TABLE]

The distribution over areas for such a ratio for varying $m$ and varying random effect distributions is reported in Figure 6. Under Scenario 1, $\widehat{\mbox{RMSE}}(\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}})$ and $\widehat{\mbox{RMSE}^{*}}(\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}})$ seem to perform generally better than the alternatives. In particular, simulation results suggest that the former estimator is more appropriate when a reduced number of small areas is available ($m=100,200$), while its precision decreases in case of larger $m$. On the other hand, $\widehat{\mbox{RMSE}^{*}}(\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}})$ shows slight overestimation of the actual Monte Carlo RMSE for $m=100$ and $m=200$, but it has to be preferred in the presence of a large number of small areas since the $R_{i}$ index is strongly concentrated around $1$.

The estimator $\widehat{\mbox{RMSE}}(\hat{p}_{i}^{\mbox{\scriptsize EBP}})$ underestimate the actual Monte Carlo RMSE, with a ratio which is always lower than $1$ for $m=100$. The quality of the results improves with $m$, even though it is always lower than that provided by the proposed approach. Such a finding may be possibly due to the estimation of the covariance matrix for parameter estimates which is not as accurate as expected with $B^{*}=250$ bootstrap resamples. In fact, it is worth noticing that Boubeta et al., (2016) highlighted the need of a very accurate estimate of the covariance matrix of parameter estimates to ensure high quality of the results. For this reason, in their simulation study, the authors suggested to estimate V(\hat{\mbox{\boldmath\Phi}}) by running a Monte Carlo experiment based on $10^{4}$ iterations in advance. In practice, when dealing with large sample sizes, such an approach is computationally very expensive and this is the reason why we considered a bootstrap approach based on $B^{*}=250$ iterations only. Last, the estimator suggested by González-Manteiga et al., (2007) overestimates the actual RMSE in all the scenarios we considered in this simulation study.

By looking at the bottom panel in Figure 6, we observe that, when dealing with non-Gaussian random effects, the MSE estimator of the sp-EBP has again good performances with an average ratio close to $1$ for all values of $m$. The effect of the bias correction term is less evident than before. When a reduced number of small areas is available, $\widehat{\mbox{RMSE}}(\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}})$ allows to estimate the actual RMSE with a higher precision than the corresponding bias-corrected version $\widehat{\mbox{RMSE}^{*}}(\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}})$. However, when $m=500$, the two estimators seem to perform similarly. The above results are not that surprising from our perspective. The bias correction term strongly relies on asymptotic results from ML theory. As a consequence, the quality of the approximation and, in turn, of the results, improves only when dealing with large sample sizes that render asymptotics more reliable. Considering that, in real applications, we expect the random effect distribution to lie in between the two “extreme” settings we considered in this simulation study and, also, that we generally need to deal with a large number of small areas, using $\widehat{\mbox{RMSE}^{*}}(\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}})$ seems to be generally more appropriate. From Figure 6, we may also notice that, under Scenario 2, the MSE estimator of the EBP works quite well (apart from being computationally prohibitive from large $m$), whereas that for the Naive estimator consistently underestimates the actual RMSE.

Furthermore, Table 4 shows the mean coverage rate (CR) for nominal $95\%$ Wald-type confidence intervals over simulations, that is

[TABLE]

As it is clear from the table, the proposed estimators show a good performance, with an average empirical coverage of approximately $92-94\%$ in all cases, except for $m=100$ under Scenario 1. On the other hand, both the EBP and, particularly, the Naive approach show a more erratic behavior. The former approach leads to under coverage for Scenario 1 and to over coverage for Scenario 2. This behavior is reversed for the Naive estimator.

To conclude, we also compared MSE estimators in terms of computational complexity. The last column of Tables 4 reports the computational time (averaged over simulations) required to get the estimates on an Intel(R) I5-3330 architecture (3.0 GHz) under each simulation setting. As it can be seen, the proposed MSE estimators show good performance also in this respect. When compared to the Naive approach, they clearly require a higher effort, which is, however, always under control. When compared to the EBP approach, the computational burden is considerably reduced. It is important to notice that, due to computational issues, results reported for the EBP approach when $m=500$ refer to a single Monte Carlo draw in place of being the average of $T=1,000$ draws as for the other methods. In this respect, it is clear that this approach does not represent an option for empirical applications with large $m$, as the one we discuss here.

When comparing the two MSE estimators we propose (with and without bias correction), we may observe that deriving $\widehat{\mbox{RMSE}^{*}}(\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}})$ requires a higher computational effort than that required for $\widehat{\mbox{RMSE}}(\hat{p}_{i}^{\mbox{\scriptsize sp-EBP}})$: this is clearly due to the computation of model derivatives in equation (4.2) which does not represent an easy task. However, such an effort is rewarded by the quality improvements we discussed so far, at least for large $m$.

7 Estimating unemployment incidence for LLMAs in Italy

In this section, we use ILFS data to estimate unemployment incidence for $611$ LLMAs in Italy. According to the simulation results in Section 6, the sp-EBP is a potentially useful approach as $(i)$ it performs better than the Naive predictor in terms of bias and efficiency; $(ii)$ it dramatically decreases the computational complexity of the MSE estimator for the parametric EBP which becomes unfeasible for a large number of small areas and/or large sample sizes. The use of the proposed approach is made easy by the availability of a (computationally efficient) algorithm for estimation and inference developed in R language from the authours. This is part of the on-line Supplementary Material at the publisher’s web-site, together with an example data set similar to the real one.

7.1 The model

To predict unemployment incidence in Italy, we considered a response variable $Y_{ij}$ taking value $1$ if unit $j$ in small area $i$ is unemployed and [math] otherwise. We followed an approach similar to that used by Molina et al., (2014) and considered the variables introduced in Section 2 and their transformations in the linear predictor, that is Sex-Age (reference = 15-24), Educational Level (reference = no education or primary school diploma), and the logarithmic transform of U-count. We ran the EM algorithm described in the Supplementary Material (Section 1) for different model specifications and a varying number of components ($G=2,\dots,6$) for the random effect distribution. The optimal solution, corresponding to the smallest AIC value, is based on $G=3$ components and includes in the linear predictor a random intercept and main covariate effects only. We report in Table 5 model parameter estimates, standard errors, and resulting p-values, together with the corresponding log-likelihood and AIC index. For comparison, we also report such quantities for the corresponding parametric model based on Gaussian random effects. Looking at this table, we may first observe that the AIC index suggests a better fit of the model based on an unspecified random effect distribution with respect to its parametric counterpart, even though differences in terms of parameter estimates are rather negligible. In particular, looking at the estimates for Sex-Age, we may notice that, when controlling for the effect of other explanatory variables in the model and for the effect of unobserved heterogeneity, the odds of being unemployed for younger people is higher than that for the older ones. For instance, the odds of being unemployed for a male in the 25-34 group are $e^{0.118}=1.125$ times those of males aged 15-24. On the other hand, for a male who is aged 35-65 years, the odds are $e^{-0.787}=0.455$, that is $54.5\%$ lower than those for the baseline category. Such differences are even stronger for females. Turning to the Educational Level, the odds of being unemployed for subjects with middle or high school diploma is higher than that of low educated subjects (parameter estimates for middle and high school diplomas are all positive). On the other hand, having a University degree or higher education has not a significant effect. These findings are in line with the results reported in the preliminary analysis: low educated females and relatively younger individuals (the reference category) are more frequent in the inactive category. Last, as expected, results reported in Table 5 suggest that the probability of being unemployed increases as the total number of unemployed registered at the $2011$ census increases.

Figure 7 shows the estimated prior (Figure 1a) and posterior distribution (Figure 1b) estimates for the random effects obtained using the proposed (semi-parametric) approach, together with the estimated posterior distribution deriving from the parametric approach (Figure 1c). In particular, in Figure 1b, we report the posterior mean of the area-specific random intercept calculated as

[TABLE]

where $\hat{\bar{\xi}}=\sum_{g}\hat{\xi}_{g}\hat{\pi}_{g}$ is the overall intercept estimatereported in Table 5. By focusing the attention on Figure 1a, we may clearly observe that observed data lead to the estimation of a random effect with a clear degree of skewness. If the standard Gaussian assumption had been reasonable, the NPML estimate of the random effect distribution would have been a symmetric distribution centered around zero. As a consequence, we may conclude that such an assumption may not be that adequate for the current application. Furthermore, by comparing Figures 1b and 1c, we may observe that the parametric assumption also affects the posterior mean of the area-specific intercepts, leading to a less skewed distribution than that obtained under the proposed approach.

7.2 Small area predictions

As highlighted in Section 3, we need the covariate values, \mbox{\boldmathx}_{ij}, to be known for all units in the population to predict the target variable. This would require access e.g. to census micro-data. However, in the important and special case where the components of \mbox{\boldmathx}_{ij} are all categorical, or take a finite number of values, the method described in this paper only requires the corresponding area level cross-tabulations to be available. This is the case of the ILFS data, where information on the covariates in the model are available at an aggregate level for the whole population. Figure 8 shows the map of unemployment incidence prediction for the $611$ LLMAs obtained using direct estimation, the proposed sp-EBP approach, the parametric EBP, and the Naive approach. Direct estimates are computed using Hájek-type estimators with adjusted weights that account for nonresponse and calibrate to population level information of demographic variables. The patterns of unemployment produced by the proposed approach are consistent with those obtained by all the other methods. As expected, model-based maps are smoother when compared to direct estimates; relatively larger values for unemployment incidences are mainly located in the South of Italy and in the Islands.

To assess the quality of predictions, we used a set of diagnostic tools based on the requirement that model-based small area estimates should be coherent with, in the sense of being close to, the corresponding unbiased direct estimates, albeit more precise. Figure 9 shows the estimates derived from the sp-EBP approach versus the direct, the EBP, and the Naive estimates, respectively. From this figure (first panel), we may observe that our approach leads to predictions which are close to those provided by a direct approach, with a correlation coefficient equal to $0.881$. From the remaining panels in Figure 9, it is evident that model-based estimates for unemployment incidence are all very close to each other, with correlation coefficients equal to $0.978$ (sp-EBP vs. EBP) and to $0.977$ (sp-EBP vs. Naive).

Coherence of direct and sp-EBP estimates can be also evaluated by computing a goodness-of-fit diagnostic (Brown et al.,, 2001), which is obtained from the following Wald-test statistic:

[TABLE]

where the estimated MSE of the sp-EBP is calculated by using formulas in Section 5. Considering the results of the simulation experiments with $m=500$, we decided to consider the bias-corrected MSE reported in equation (13). The above test is based on the idea that, should model-based estimates be close to the “true” small area parameters of interest, the unbiased direct estimates could be considered as random variables with expected value equal to the value of the corresponding model-based estimates. Here, $W=360.56$ and such a value needs to be compared to the $95$-th percentile of a $\chi^{2}$ distribution with $452$ d.f., $\chi^{2}_{452,0.95}=502.56$. In this respect, we may conclude that model-based estimates are not significantly different from direct estimates.

To assess the potential gain in precision we obtain by using the proposed sp-EBP approach in place of the direct one, we compare in Figure 10 the empirical cumulative density functions (ecdfs) of the vcoefficients of variation (CV) of both estimators. The first panel uses CV’s from all areas, while the second (third) one focuses on small areas with sample size smaller (higher) than $100$. As it is clear, by looking at the first panel, the ecdf corresponding to sp-EBP almost always dominates the one for the direct estimates, highlighting that CV values for the former approach are lower than those estimated with the latter. Only for very small CV values, the ecdfs show an inverse relation: CV’s for direct predictions are smaller than those for sp-EBP. This is more evident in the second panel and is related to the presence of some areas with a small sample size for which $\hat{p}_{i}^{\mbox{\scriptsize Direct}}$ is zero or is very close to zero, and so is $\widehat{\mbox{RMSE}}(\hat{p}_{i}^{\mbox{\scriptsize Direct}})$. However, also in this case, about the $60\%$ of the Italian LLMAs, CVs associated to the direct estimator are above the standard $33\%$ threshold which is typically considered for reliability in the SAE context. Such a percentage reduces to about $20\%$ when considering the proposed sp-EBP approach; in addition, less than 5% of the estimates have CV% larger than 40%. When we move to higher CVs, the sp-EBP approach always provides smaller CV values when compared to the direct approach and such CVs are always smaller than $40\%$. When focusing on the third panel in Figure 10, it is evident that, as the sample size gets larger, direct and model-based estimates tend to have quite similar CV values, although those associated with model-based estimates are still consistently smaller.

8 Conclusions

The paper described some tools to derive Best Predictions for responses with distribution in the Exponential Family in the presence of clustered data. In particular, we proposed a semi-parametric version of the EBP and the corresponding second-order, bias-corrected, MSE approximation using a NPML approach and leaving the distribution of the random effects unspecified. Motivated by a real application to data on unemployment incidence in LLMAs in Italy, we focused on a binary response modeled via a mixed logistic model with random intercepts, which represents a relevant case in the SAE framework.

Simulation experiments showed that the proposed estimator performs equally or better than the competitors. In particular, when moving far from the assumption of Gaussian distributed random effects, the proposed semi-parametric approach performed better than the corresponding parametric versions. Also, when compared to the parametric EBP, simulation results highlighted better performance of the proposed approach in terms of computational load required to get predictions and the corresponding MSE. The simulation study, where different sample sizes were considered, showed that the semi-parametric approach is always reliable, especially for large $m$. Such a gain comes from the discrete nature of the mixing distribution estimate which substantially simplifies the calculations to get the EBPs and the corresponding MSEs.

We illustrated the benefits of our proposal discussing the estimation of unemployment incidence for Italian LLMAs in Italy in $2012$. In this context, direct estimates cannot be published for most of the LLMAs given the unacceptable large value of the coefficient of variation for those areas with a small sample size. In this respect, model-based approaches represent a necessary strategy. Since the sample size and the number of small areas are particularly large in this application, the implementation of the EBP turns to be particularly cumbersome, and the evaluation of its precision prohibitive. This application indicated that the proposed methodology leads to estimates which are coherent with, but more efficient than, the direct estimates, still being comparable with alternative model-based estimates.

Although the approach we propose is presented for responses with density in the Exponential Family, we did not explore the behavior of the small area sp-EBPs for counts or multinomial responses. However, a possible extension to multi-category outcomes is quite straightforward. Also, we notice that suitable extensions of the proposed approach to allow for spatial correlation could be envisioned by properly modeling, for each small area, the prior mixture probabilities as a function of neighborhood components membership. Last, developing design-consistent small area estimators under the proposed methodology represent a topic of interest, especially for those researchers working in survey sampling from a design-based or a model-assisted perspective. More specifically, we could adopt a model-assisted approach, thereby the model is used only to motivate the predictors, but their properties are evaluated only with respect to the randomization distribution induced by the sampling design.

Acknowledgements

The work by Marino, Ranalli and Salvati has been developed under the support of the project PRIN-SURWEY (grant 2012F42NS8, Italy). The work of Alfó has been supported by the project “Mixture and latent variable models for causal inference and analysis of socio-economic data” (grant RBFR12SHVV, FIRB - Futuro in Ricerca). The work of Salvati was carried out with the support of the project InGRID 2 (Grant Agreement No 730998, EU).