Multivariate Small Area Estimation of Social Indicators: The Case of Continuous and Binary Variables

Abstract

Large-scale sample surveys are not designed to produce reliable estimates for small areas. Here, small area estimation methods can be applied to estimate population parameters of target variables to detailed geographic scales. Small area estimation for noncontinuous variables is a topic of great interest in the social sciences where such variables can be found. Generalized linear mixed models are widely adopted in the literature. Interestingly, the small area estimation literature shows that multivariate small area estimators, where correlations among outcome variables are taken into account, produce more efficient estimates than do the traditional univariate techniques. In this article, the author evaluate a multivariate small area estimator on the basis of a joint mixed model in which a small area proportion and mean of a continuous variable are estimated simultaneously. Using this method, the author “borrows strength” across response variables. The author carried out a design-based simulation study to evaluate the approach where the indicators object of study are the income and a monetary poverty (binary) indicator. The author found that the multivariate approach produces more efficient small area estimates than does the univariate modeling approach. The method can be extended to a large variety of indicators on the basis of social surveys.

Large-scale sample surveys are usually not designed to produce accurate and precise estimates for small population domains, for example, small geographic areas or small groups in the population (e.g., on the basis of cross-classification between ethnic groups and age bands). However, many social phenomena, such as poverty, well-being, and social exclusion show an important spatial heterogeneity at a small geographic level (Molina and Strzalkowska-Kominiak 2020; Moretti, Shlomo, and Sakshaug 2020; Pratesi 2016). Importantly, policymakers in charge of implementing subnational social policies ask for disaggregated estimates of social indicators. Thus, small area estimates of such phenomena should be produced and released to support decision makers (Pratesi 2016).

In recent years, there has been growing attention to the methodological development of small area estimation methods. Statistical models are now prominent in this field. Indeed, model-based estimators “borrow strength” from related areas using linking models, and they provide more efficient estimates than do direct estimators (see Pfeffermann 2013; Rao and Molina 2015).

The small area estimation approach using linear mixed models is now widely used. These models are powerful given that they include random area effects to take into account between-area variation, that is, the unexplained variability between the small areas. A pioneer work in this context is the Battese, Harter, and Fuller model (Battese, Harter, and Fuller 1988), with further work on mean squared error (MSE) estimation proposed by Prasad and Rao (1990). This model is applicable to continuous response variables, and normality is assumed for the individual error and random area effects. In case of failure of the model assumptions, extensions are proposed in the literature, for example, to account for heteroscedasticity and outliers. Readers may refer to Rao and Molina (2015) for an extensive review of methodological developments.

In the presence of binary outcomes, the generalized linear mixed model (GLMM) with logit link function, that is, a logistic mixed model, is widely adopted. Specifically, an empirical plug-in predictor (EPP) under a GLMM is used in small area estimation of proportions in official statistics (Chandra, Chambers, and Salvati 2012; Chandra, Kumar, and Aditya 2018; Molina and Strzalkowska-Kominiak 2020; Rao and Molina 2015; Salvati, Chandra, and Chambers 2012). For example, the U.K. Office for National Statistics and the Australian Bureau of Statistics use this method (Chambers, Salvati, and Tzavidis 2016; Chandra et al. 2018). Binary variables can be found in most surveys. For instance, some poverty and well-being indicators are based on binary variables (Betti and Lemmi 2013). Additionally, some social indicators estimated on Labour Force Surveys are also constructed on these types of variables (see Chambers et al. 2016; Molina and Strzalkowska-Kominiak 2020). Hence, there is a high demand for small area proportions.

Because many social phenomena are naturally multidimensional (Betti and Lemmi 2013), and therefore correlated, this property can be used to further improve small area estimates (Benavent and Morales 2016; Fabrizi, Ferrante, and Pacei 2005; Moretti et al. 2020; Moretti, Shlomo, and Sakshaug 2021; Ubaidillah et al. 2019). In this context, multivariate small area estimation methods can be applied. Moretti et al. (2020) proposed the use of multivariate small area estimation methods to estimate latent well-being indicators, under the use of the multivariate extension of the Battese, Harter, and Fuller model. This model was studied in detail by Datta, Day, and Basawa (1999). Recently, Moretti (2022) extended a GLMM with logit link function to the case of bivariate proportions, showing good results in terms of efficiency compared with its univariate setting.

In this article, I investigate the multivariate small area estimation problem of continuous and binary response variables jointly. In particular, I apply a joint mixed-modeling strategy to the small area estimation of the mean of a continuous variable and proportion of a binary variable. This is an important problem in sociology and related social sciences. In fact, when building poverty indicators, they can be based on income and measured by a binary indicator taking a value of 0 or 1, denoting whether the personal or household income falls below a poverty threshold (for a review and discussion of poverty indicators, see Pratesi 2016). Therefore, income and poverty indicators are correlated; this can be accounted for in the small area models by assuming a joint modeling strategy.

Notation and Multivariate Small Area Estimation Problem

I now present the notation used in the article and the small area estimation problem we are studying. Let us consider a finite target population $U$ with size $N$ that is partitioned into $D$ nonoverlapping (disjoint) small areas, $U_d$ , $d=1,...,D$ of size $N_d$ such that ${\cup }_{d=1}^D{U}_d=U$ and $\sum _{d=1}^D{N}_d=N$ . A random sample $s$ of size $n$ is selected from $U$ . $n_d$ denotes the sample size in small area $d$ , and $\sum _{d=1}^D{n}_d=n$ .

Let ${\mathbf{y}}_{\mathrm{di}}=(y_{\mathrm{di}1},y_{\mathrm{di}2}{)}^T$ denote a vector of the values of $k=1,2$ variables of interest $\mathbf{Y}$ for unit $i$ in area $d$ . We assume $y_{\mathrm{di}1}$ is continuous and normally distributed, whereas $y_{\mathrm{di}2}$ takes value 0 or 1 only. We are interested in estimating the vector of means of $\mathbf{Y}$ denoted by ${\overline{\mathbf{Y}}}_d=({\overline{Y}}_{d1},{\overline{Y}}_{d2})$ . ${\overline{Y}}_{d1}$ is the mean of the continuous variable $Y_1$ , and ${\overline{Y}}_{d2}$ is the proportion related to variable $Y_2$ , which is also a mean.

The generic element related to variable $k$ is given by equation (1):

$${\overline{Y}}_{\mathrm{dk}}=N_d^{-1}\sum _{i\in U_d}{y}_{\mathrm{dik}}=N_d^{-1}(\sum _{i\in s_d}{y}_{\mathrm{dik}}+\sum _{i\in r_d}{y}_{\mathrm{dik}}),$$

(1)

where $s_d$ denotes the sample elements and $r_d$ the out-of-sample elements in area $d$ . Looking at equation (1), we see that the population mean can be split into sample $(s_d)$ and out-of-sample elements $(r_d)$ in area $d$ .

The design-based direct estimator for the kth small area mean ${\overline{Y}}_{\mathrm{dk}}$ is given by

$${\hat{\overline{Y}}}_{\mathrm{dk}}^{\mathrm{DIR}}=\frac{\sum _{i\in s_d}{w}_{\mathrm{di}}{y}_{\mathrm{dik}}}{\sum _{i\in s_d}{w}_{\mathrm{di}}},$$

(2)

where $w_{\mathrm{di}}$ denotes the survey weight for unit $i$ in area $d$ .

It is well known that estimator (2) becomes unstable when $n_d$ is small, because it is based on area-specific sample quantities only (see Rao and Molina 2015; Särndal, Swensson, and Wretman 2003). Additionally, the estimator cannot be computed for small areas with $n_d=0$ . Therefore, model-based small area estimation methods that “borrow strength” across small areas via statistical models are adopted to produce reliable, that is, accurate and precise, small area estimates of the target parameter given by equation (1) (Rao and Molina 2015). For the estimator of the variance of equation (2), see Rao and Molina (2015).

Small Area Estimation under a Joint Mixed Model for Continuous and Binary Outcome Variables

To build indirect small area estimators, statistical models with random area-specific effects that account for between-area variability are used in small area estimation (Rao and Molina 2015). In the univariate small area estimation setting, the Battese, Harter, and Fuller model (Battese et al. 1988) can be used to obtain small area estimates of the population mean in case of continuous variables, that is, ${\overline{Y}}_{d1}$ . In the case of binary responses, the GLMM with logit link function is often used to estimate small area proportions, ${\overline{Y}}_{d2}$ (Chambers et al. 2016; Chandra et al. 2018). Prior literature shows that use of multivariate mixed models in small area estimation provides more efficient estimates than does the univariate case where separate models are estimated on each response variable (Datta et al. 1999; Moretti et al. 2020; Rao and Molina 2015).

Here, we assume that ${\mathbf{y}}_{\mathrm{di}}$ follows a joint mixed effects model. This is obtained by assuming a multivariate distribution for the random effects of responses $k=1$ and $k=2$ (Ivanova, Molenberghs, and Verbeke 2016). The model is presented below, for $i=1,...,N_d$ and $d=1,...,D$ :

$$\{\begin{array}{c}{y}_{\mathrm{di}1}={\mathbf{x}}_{\mathrm{di}1}^T{\mathbf{\beta }}_1+e_{\mathrm{di}1}+u_{d1}\hfill \\ \multicolumn{1}{c}{\mathrm{logit}\left[{\pi }_{\mathrm{di}2}\right]=\log \left[\frac{{\pi }_{\mathrm{di}2}}{1-{\pi }_{\mathrm{di}2}}\right]={\eta }_{\mathrm{di}2}={\mathbf{x}}_{\mathrm{di}2}^T{\mathbf{\beta }}_2+u_{d2}}\end{array},$$

(3)

where ${\mathbf{x}}_{\mathrm{di}1}$ and ${\mathbf{x}}_{\mathrm{di}2}$ denote the vectors of observed values of $p$ unit-level auxiliary variables for unit $i$ in area $d$ related to variables $k=1$ and $k=2$ . These can be the same for both responses, depending on data availability and specific modeling problems (Fuller and Harter 1987; Ivanova et al. 2016; Moretti et al. 2020).

Model 3 is a multilevel model for multivariate mixed response types, explaining two outcomes simultaneously (Goldstein 2011), in this case $Y_1$ and $Y_2$ . Goldstein et al. (2009) studied a similar modeling approach. Because we are aiming to predict the small area means of the $Y_1$ and $Y_2$ variables, we cannot account for the correlation structure between the two by including one of them as a predictor in a univariate model instead of a joint modeling approach. In fact, in the model-based small area estimation approach, we need auxiliary information of the predictors for all the units in the population to build the EPP. These are not available from census or administrative data for $Y_1$ and $Y_2$ (e.g., the income variable is not available as an auxiliary variable from external data sources).

We assume multivariate (bivariate in this case) normality of the random area effects, that is, ${\mathit{u}}_d=(u_{d1},u_{d2}{)}^T~N(0,{\mathit{\Sigma }}_u)$ , where ${\mathbf{\Sigma }}_u$ denotes a $2\times 2$ positive-definite variance-covariance matrix. The off-diagonal elements are the covariances between $u_{d1}$ and $u_{d2}$ . The matrix is given as follows:

$${\mathbf{\Sigma }}_u=\left[\begin{array}{cc}{\sigma }_{u1}^2\hfill & {\rho }_u·\sqrt{{\sigma }_{u1}^2·{\sigma }_{u2}^2}\hfill \\ \multicolumn{1}{c}{{\rho }_u·\sqrt{{\sigma }_{u2}^2·{\sigma }_{u1}^2}}& {\sigma }_{u2}^2\hfill \end{array}\right],$$

where ${\rho }_u$ denotes the correlation coefficient. In addition, we assume that $y_{\mathrm{di}2}|u_{d2}~\mathrm{Binomial}(1,{\pi }_{\mathrm{di}2})$ with ${\pi }_{\mathrm{di}2}=E(y_{\mathrm{di}2}|u_{d2})$ .

Model 3 can be written for $i=1,...,n_d$ without loss of generality. Therefore, the model parameters are estimated on a random sample s drawn from the population $U$ (Rao and Molina 2015). Here, I estimate model parameters via the maximum likelihood (ML) approach. For estimation details that are beyond the scope of the present article, readers can refer to McCulloch (1994, 1997) and Booth and Hobert (1999) for the theory, and Berridge and Crouchley (2011) for its implementation.

Because of the multivariate nature of the likelihood function, a Gaussian quadrature has been adopted in the literature. ML estimation of the parameters is computationally complex for this model (Berridge and Crouchley 2011; Coull and Agresti 2000). Readers interested in the computations can refer to Appendix A.3 in Berridge and Crouchley (2011). The implementation is fully available in R via the package sabre. ¹ This approach is also studied and evaluated in Rabe-Hesketh and Skrondal (2008) and Skrondal and Rabe-Hesketh (2004).

To provide a small area estimator for ${\overline{\mathbf{Y}}}_d$ , we consider an EPP. Empirical best predictors in the case of GLMMs are difficult to obtain analytically, and they are not available in a closed form. EPPs are easier to obtain and widely used by statistical agencies (e.g., the U.K. Office for National Statistics and the Australian Bureau of Statistics) (Chambers et al. 2016; Chandra et al. 2018). This problem becomes even more challenging in the case of joint modeling of mixed types responses.

Therefore, we can write the EPP of the small area means under model 3 for area $d$ as follows:

$$\{\begin{array}{c}{\hat{\overline{Y}}}_{d1}^{\mathrm{EPP}}=N_d^{-1}(\sum _{i\in s_d}{y}_{\mathrm{di}1}+\sum _{i\in r_d}{\mathbf{x}}_{\mathrm{di}1}^T{\hat{\mathbf{\beta }}}_1+{\hat{u}}_{d1})\hfill \\ \multicolumn{1}{c}{{\hat{\overline{Y}}}_{d2}^{\mathrm{EPP}}=N_d^{-1}(\sum _{i\in s_d}{y}_{\mathrm{di}2}+\sum _{i\in r_d}{\hat{\mu }}_{\mathrm{di}2})}\end{array},$$

(4)

with ${\hat{\mu }}_{\mathrm{di}2}=\hat{E}(y_{\mathrm{di}2}|u_{d2})=\exp ({\mathbf{x}}_{\mathrm{di}2}^T{\hat{\mathbf{\beta }}}_2+{\hat{u}}_{d2})[1+\exp ({\mathbf{x}}_{\mathrm{di}2}^T{\hat{\mathbf{\beta }}}_2+{\hat{u}}_{d2}){]}^{-1}$ , where ${\hat{\mathbf{\beta }}}_1$ , ${\hat{\mathbf{\beta }}}_2$ , ${\hat{u}}_{d1}$ , and ${\hat{u}}_{d2}$ denote the estimates of the regression coefficients and predictions of random effects, respectively.

In practical applications, the auxiliary variables are available at the unit level for the sample units only, and area-specific aggregates are available for the population from the census or administrative data. As a consequence, equation (4) cannot be applied; however, a modification of it can be derived, and its performance is studied in the literature (Chandra et al. 2018). This is given as follows:

$$\{\begin{array}{c}{\hat{\overline{Y}}}_{d1}^{\mathrm{EPP}1}=f_d{\hat{\overline{Y}}}_{d1}+(1-f_d)({\overline{\mathbf{X}}}_{r,d1}^T{\hat{\mathbf{\beta }}}_1+{\hat{u}}_{d1})\hfill \\ \multicolumn{1}{c}{{\hat{\overline{Y}}}_{d2}^{\mathrm{EPP}1}=f_d{\hat{\overline{Y}}}_{d2}+(1-f_d)\exp ({\overline{\mathbf{X}}}_{r,d2}^T{\hat{\mathbf{\beta }}}_2+{\hat{u}}_{d2})[1+\exp ({\overline{\mathbf{X}}}_{r,d2}^T{\hat{\mathbf{\beta }}}_2+{\hat{u}}_{d2}{)]}^{-1}}\end{array},$$

(5)

where $f_d=n_d/N_d$ is the sampling fraction in area $d$ and ${\hat{\overline{Y}}}_{\mathrm{dk}}=\sum _{i\in s_d}{y}_{\mathrm{dik}}/n_d$ . ${\overline{\mathbf{X}}}_{r,\mathrm{dk}}=(N_d-n_d{)}^{-1}(N_d{\overline{\mathbf{X}}}_{\mathrm{dk}}-n_d{\overline{\mathbf{x}}}_{\mathrm{dk}})$ denotes the means of the auxiliary variables for the out-of-sample units, with ${\overline{\mathbf{X}}}_{\mathrm{dk}}=N_d^{-1}\sum _{i\in U_d}{\overline{\mathbf{x}}}_{\mathrm{dik}}$ and ${\overline{\mathbf{x}}}_{\mathrm{dk}}=n_d^{-1}\sum _{i\in s_d}{\mathbf{x}}_{\mathrm{dik}}$ , denoting the means of the auxiliary variables in the population and sample, respectively, for area $d$ and $k=1,2$ .

When the sampling fractions $f_d$ are small (negligible) and ${\overline{\mathbf{X}}}_{\mathrm{dk}}\approx {\overline{\mathbf{X}}}_{r,\mathrm{dk}}$ , the EPPs given in equation (5) can be written as follows (Chandra et al. 2018; Moretti 2022):

$$\{\begin{array}{c}{\hat{\overline{Y}}}_{d1}^{\mathrm{EPP}1.1}=f_d{\hat{\overline{Y}}}_{d1}+(1-f_d)({\overline{\mathbf{X}}}_{d1}^T{\hat{\beta }}_1+{\hat{u}}_{d1})\hfill \\ \multicolumn{1}{c}{{\hat{\overline{Y}}}_{d2}^{\mathrm{EPP}1.1}=f_d{\hat{\overline{Y}}}_{d2}+(1-f_d)\exp ({\overline{\mathbf{X}}}_{d2}^T{\hat{\beta }}_2+{\hat{u}}_{d2})[1+\exp ({\overline{\mathbf{X}}}_{d2}^T{\hat{\beta }}_2+{\hat{u}}_{d2}{)]}^{-1}.}\end{array}$$

(6)

This is reasonable to assume in real data small area applications, given that $f_d$ may be very small. Hence, in this article, I consider and evaluate EPP1.1 only.

MSE Estimation via Parametric Bootstrap

The MSE of equation (6) can be estimated via a parametric bootstrap algorithm. The parametric bootstrap is widely studied and applied in the small area estimation literature, where simulation studies are designed to evaluate the algorithm’s performance (González-Manteiga et al. 2007; Hobza, Morales, and Santamaría 2018; Moretti, Shlomo, and Sakshaug 2020). Here, I extend the algorithm originally developed in González-Manteiga et al. (2007) to the case considered in this article. Note that Moretti et al. (2020) extended and evaluated the algorithm in González-Manteiga et al. (2007) to the bivariate small area estimation of means problem in the case of continuous variables.

The steps of the parametric bootstrap algorithm are as follows:

1. Estimate model given in equation (3) on the random sample s. The following estimates are obtained: ${\widehat{\mathbf{\Sigma }}}_u$ ${\hat{\mathbf{\beta }}}_1$ , ${\hat{\mathbf{\beta }}}_2$ , and ${\hat{\sigma }}_e^2$ .

2. Generate the bootstrap area-specific random effects as follows: ${\mathit{u}}_d^{*(b)}~N(0,\widehat{{\Sigma }_u})$ , and individual error term for response $k=1$ only, that is, $e_{\mathrm{di}1}^{*(b)}~N(0,{\hat{\sigma }}_e^2)$ . The asterisk denotes the bootstrap quantities, and $b$ denotes the bth bootstrap replication, $b=1,...,B$ .

3. Calculate the true means of the bootstrap population for variables $k=1,2$ and area $d$ :

   $$\{\begin{array}{c}{\overline{Y}}_{d1}^{*(b)}={\overline{\mathbf{X}}}_{d1}{\hat{\beta }}_1+u_{1d}^{*(b)}\hfill \\ \multicolumn{1}{c}{{\overline{Y}}_{d2}^{*(b)}=\exp ({\overline{\mathbf{X}}}_{d2}{\hat{\beta }}_2+u_{d2}^{*(b)})[1+\exp ({\overline{\mathbf{X}}}_{d2}{\hat{\beta }}_2+u_{d2}^{*(b)}{)]}^{-1}}\end{array}$$

   (7)

4. Generate the bootstrap responses $y_{\mathrm{di}1}^{*(b)}$ and $y_{\mathrm{di}2}^{*(b)}$ according to model 3 for $i\in s_d$ as follows:

   $$\{\begin{array}{c}{y}_{\mathrm{di}1}^{*(b)}={\overline{\mathit{x}}}_{\mathrm{di}1}^T{\hat{\beta }}_1+e_{\mathrm{di}1}^{*(b)}+u_{d1}^{*(b)}\hfill \\ \multicolumn{1}{c}{y_{\mathrm{di}2}^{*(b)}|u_{d2}^{*(b)}~\mathrm{Binomial}(1,{\pi }_{\mathrm{di}2}^{*(b)})}\end{array},$$

   (8)

   where ${\pi }_{\mathrm{di}2}^{*(b)}=\exp ({\overline{\mathbf{x}}}_{\mathrm{di}2}^T{\hat{\beta }}_2+u_{d2}^{*(b)})[1+\exp ({\overline{\mathbf{x}}}_{\mathrm{di}2}^T{\hat{\beta }}_2+u_{d2}^{*(b)}){]}^{-1}$ .

5. Estimate the joint mixed model given in equation (3) on the responses generated at step 4 and obtain the bootstrap EPP1.1 according to equation (6). These are denoted by ${\hat{\overline{Y}}}_{d1}^{\mathrm{EPP}1.1*(b)}$ and ${\hat{\overline{Y}}}_{d2}^{\mathrm{EPP}1.1*(b)}$ .

6. Repeat steps 2 to 5 $B$ times.

7. The bootstrap estimators for the MSE of ${\hat{\overline{Y}}}_{d1}^{\mathrm{EPP}1.1*(b)}$ and ${\hat{\overline{Y}}}_{d2}^{\mathrm{EPP}1.1*(b)}$ are given by

$$\{\begin{array}{c}{\widehat{\mathrm{MSE}}}_{\mathrm{boot}}\left({\hat{\overline{Y}}}_{d1}^{\mathrm{EPP}1.1*(b)}\right)=B^{-1}\sum _{b=1}^B{({\hat{\overline{Y}}}_{d1}^{\mathrm{EPP}1.1*(b)}-{\overline{Y}}_{d1}^{*(b)})}^2\hfill \\ \multicolumn{1}{c}{{\widehat{\mathrm{MSE}}}_{\mathrm{boot}}\left({\hat{\overline{Y}}}_{d2}^{\mathrm{EPP}1.1*(b)}\right)=B^{-1}\sum _{b=1}^B{({\hat{\overline{Y}}}_{d2}^{\mathrm{EPP}1.1*(b)}-{\overline{Y}}_{d2}^{*(b)})}^2,}\end{array}$$

(9)

for $d=1,...,D$ , with $B$ denoting the number of replicates, that is, $B=500$ (Hobza et al. 2018).

Simulation Study

In this section, I present the results of a design-based simulation study. Design-based simulation studies are important because they allow one to evaluate the performance of the estimators in case of repeated samples drawn from a fixed population that does not exactly follow the assumed model (Molina and Strzalkowska-Kominiak 2020).

Generating the Population

As a fixed target population $U$ , I use data from Lehtonen and Veijanen (2016). This is available online from Pratesi (2016). These data were derived from AMELIA data (see Burgard et al. 2017; Lehtonen and Veijanen 2016). AMELIA is a synthetic population that allows comparative and reproducible research on the basis of European Union Statistics for Income and Living Conditions variables (Burgard et al. 2017).

The population size in area $d$ ranges between 12,550 and 42,340, and the areas are $D=40$ . As variables of interest, I consider the logarithmic transformation of the income, $Y_1$ , and a binary variable taking value 1 if the unit $i$ in area $d$ is poor and 0 otherwise, $Y_2$ . The unit $i$ is poor if the value of the income is below the poverty line, calculated as 60 percent of the median of the income (Chatterjee 2011). The auxiliary variable, $X_1$ , is the age of unit $i$ . Because of the nature of the response variables, the random effects are highly correlated with $\hat{\rho }=-0.80$ . $\hat{\rho }$ is obtained via a joint mixed model estimated on the fixed population $U$ .

Figure 1 shows the histogram of the logarithmic transformation of income, $Y_1$ , and Figure 2 shows the bar chart of the binary variable, $Y_2$ , denoting whether the unit is poor or not. $Y_1$ is approximately normally distributed, but slightly skewed and with some outliers; hence, it mimics real data distributions.

Figure 1.

Histogram of log income, $Y_1$ .

Figure 2.

Bar chart of indicator denoting whether the unit’s income is below the poverty line, $Y_2$ .

Simulation Steps and Quality Measures

The simulation consists of the following steps, where $l=1,...,L$ ; $L=500$ denotes the repetition:

1. From the target population $U={\cup }_{d=1}^D{U}_d$ where $U_d=\{(y_{\mathrm{di}1},y_{\mathrm{di}2},x_{\mathrm{di}}),i=1,...,N_d\}$ , select a sample $s_d^{(l)}$ without replacement of size $n_d=3$ for $d=1,...,D$ . Note that $U$ is fixed over the simulations because we are following a design-based simulation setting.

2. Estimate the joint mixed model given in equation (3) and the univariate mixed models on each separate response in each sample $s_d^{(l)}$ and obtain the univariate and multivariate predictors for both small area means. Readers can refer to the Battese, Harter, and Fuller model (Battese et al. 1988) for the univariate setting and to Chandra et al. (2018) for the small area estimation problem of proportions. The multivariate small area predictors are given by equation (6). The predictors are denoted by ${\hat{\overline{y}}}_{d1}^{\mathrm{MEPP}1.1(l)}$ , ${\hat{\overline{y}}}_{d2}^{\mathrm{MEPP}1.1(l)}$ and ${\hat{\overline{y}}}_{d1}^{\mathrm{UEPP}1.1(l)}$ , ${\hat{\overline{y}}}_{d2}^{\mathrm{UEPP}1.1(l)}$ , for the multivariate and univariate case, respectively, and lth repetition. The direct estimates are also calculated via equation (2) and denoted by ${\hat{\overline{y}}}_{d1}^{\mathrm{DIR}(l)}$ , ${\hat{\overline{y}}}_{d2}^{\mathrm{DIR}(l)}$ .

3. The following measures of performance are calculated to evaluate the estimators for $k=1,2$ in both the univariate and multivariate case (here, ${\hat{\overline{y}}}_{\mathrm{dk}}^{(l)}$ denotes any estimator for the true value in the population, i.e., ${\overline{y}}_{\mathrm{dk}}$ , for kth variable and dth area):

Absolute relative bias (ARB)

$$\mathrm{ARB}({\hat{\overline{y}}}_{\mathrm{dk}})=\left|\frac{L^{-1}\sum _{l=1}^L({\hat{\overline{y}}}_{\mathrm{dk}}^{(l)}-{\overline{y}}_{\mathrm{dk}})}{{\overline{y}}_{\mathrm{dk}}}\right|,$$

(10)

Root MSE (RMSE)

$$\mathrm{RMSE}({\hat{\overline{y}}}_{\mathrm{dk}})=\sqrt{L^{-1}\sum _{l=1}^L{({\hat{\overline{y}}}_{\mathrm{dk}}^{(l)}-{\overline{y}}_{\mathrm{dk}})}^2},$$

(11)

Relative root MSE (RRMSE)

$$\mathrm{RRMSE}({\hat{\overline{y}}}_{\mathrm{dk}})=\frac{\mathrm{RMSE}({\hat{\overline{y}}}_{\mathrm{dk}})}{{\overline{y}}_{\mathrm{dk}}},$$

(12)

where ${\overline{y}}_{\mathrm{dk}}$ denotes the true values for $k=1,2$ .

% Relative reduction in terms of RMSE (RelRed%)

$$\mathrm{RelRed}({\hat{\overline{y}}}_{\mathrm{dk}})\mathrm{\%}=L^{-1}\sum _{l=1}^L\frac{\mathrm{RMSE}({\hat{\overline{y}}}_{\mathrm{dk}}^{\mathrm{MEPP}1.1(l)})-\mathrm{RMSE}({\hat{\overline{y}}}_{\mathrm{dk}}^{\mathrm{UEPP}1.1(l)})}{\mathrm{RMSE}({\hat{\overline{y}}}_{\mathrm{dk}}^{\mathrm{UEPP}1.1(l)})}\times 100.$$

(13)

To present summary statistics, the median across the small areas $D$ is shown as a robust central tendency measure that avoids the effect of extreme values in some small areas (Chambers, Chandra, and Tzavidis 2011; Giusti et al. 2014). In this case, the same notation as above is used but the index $d$ is dropped.

Results

I now present the results of the simulation study with a focus on the quality measures presented above. Table 1 shows the median across the small areas of the ARB and RRMSE of the estimates obtained via the direct, univariate, and multivariate estimators for $k=1$ and $k=2$ means.

Table 1.

Median Values across Small Areas of ARB and RRMSE of the Estimates Obtained via the Direct, Univariate, and Multivariate Estimators for $k=1$ and $k=2$ Means

Estimator	ARB	RRMSE
${\hat{\overline{y}}}_1^{\mathrm{DIR}}$	.010	.091
${\hat{\overline{y}}}_2^{\mathrm{DIR}}$	.030	.842
${\hat{\overline{y}}}_1^{\mathrm{UEPP}}$	.018	.030
${\hat{\overline{y}}}_2^{\mathrm{UEPP}}$	.100	.164
${\hat{\overline{y}}}_1^{\mathrm{MEPP}}$	.014	.020
${\hat{\overline{y}}}_2^{\mathrm{MEPP}}$	.077	.150

Note: ARB = absolute relative bias; RRMSE = relative root mean squared error.

As shown, the ARB is small (negligible) across all the estimators. The multivariate approach produces estimates with a smaller ARB compared with the univariate case. Thus, it does not introduce much bias in the small area estimates. Furthermore, the multivariate small area predictor provides estimates with a smaller RRMSE compared with the univariate and direct estimators. Hence, the multivariate approach returns more efficient estimates compared with the univariate approach. This latter point can be seen if we consider the percentage relative reductions in terms of RMSE. Indeed, these are satisfactory and equal to $\mathrm{RelRed}({\hat{\overline{y}}}_{d1})\mathrm{\%}=-35\mathrm{\%}$ and $\mathrm{RelRed}({\hat{\overline{y}}}_{d2})\mathrm{\%}=-10\mathrm{\%}$ .

Figures 3 and 4 plot the RRMSE across the small areas where we compare the RRMSE of the multivariate versus univariate approach and their direct estimators setting for $k=1$ and $k=2$ , respectively. Figures 5 and 6 show the plots of the ARB across the small areas for the three approaches for $k=1$ and $k=2$ , respectively. Descriptive statistics related to these figures are shown in Table 1.

Figure 3.

Relative root mean squared error (RRMSE) of small area estimates for direct (DIR; dotted line), univariate (UEPP; dashed line), and multivariate (MEPP; solid line) estimators for the mean of log income, k = 1, ordered by decreasing sampling fraction.

Figure 4.

Relative root mean squared error (RRMSE) of small area estimates for direct (DIR; dotted line), univariate (UEPP; dashed line), and multivariate (MEPP; continuous line) estimators for the poverty proportion, k = 2, ordered by decreasing sampling fraction.

Figure 5.

Absolute relative bias (ARB) of small area estimates for direct (DIR; dotted line), univariate (UEPP; dashed line), and multivariate (MEPP; solid line) estimators for the mean of log income, k = 1, ordered by decreasing sampling fraction.

Figure 6.

Absolute relative bias (ARB) of small area estimates for direct (DIR; dotted line), univariate (UEPP; dashed line), and multivariate (MEPP; solid line) estimators for the poverty proportion, k = 2, ordered by decreasing sampling fraction.

As expected, Figures 3 to 6 demonstrate that the direct estimates suffer from a large RRMSE but very small ARB across all the small areas. The multivariate approach provides estimates with a smaller RRMSE than both the univariate and direct approaches, which is in line with the multivariate small area estimation literature. In addition, the ARB of the multivariate estimates is negligible across all the areas, showing it does not introduce a large bias in the estimates. I did not find any relationship between the sampling fraction and performances of the multivariate approach over its univariate setting. Figures A1 and A2 and Table A1 in Appendix A provides additional simulation outputs that are helpful for interpreting the results.

In summary, the results of this simulation study are promising. Although the income variable, $Y_1$ , does not follow a normal distribution perfectly (see Figure 1), as it occurs with real data, the multivariate approach provides estimates with a smaller RRMSE than the univariate case, and the ARB is negligible across the small areas. Potential applications of this approach to social indicators are discussed in the next section.

Conclusion and Future Research Directions

In this article, I studied the multivariate small area estimation problem in the case of mixed-type variables. The multivariate small area estimation literature has paid particular attention to the case of continuous variables. However, noncontinuous variables are widely diffused in social surveys. In this work, I considered the unit-level small area estimation approach, assuming the auxiliary variables are known for all the sampled units. I did not investigate the area-level approach in the present article. For the area-level approach, the reader may refer to Fay (1987), and more recent studies by Benavent and Morales (2016, 2020).

I compared a multivariate EPP with its univariate setting, where two mixed models are estimated separately on each response variable in a design-based simulation study. The predictors based on the univariate setting are used by different statistical agencies due to their good properties. I found the multivariate approach provides more reliable small area estimates than does the univariate approach. In particular, the multivariate predictor produces small area estimates with a smaller RRMSE, as well as bias slightly smaller than the univariate approach. Larger gains in efficiency can be seen for the continuous case.

The modeling strategy evaluated in this article is flexible, because joint mixed models can account for variables measured on different scales at the same time. Thus, we can improve the efficiency of the small area estimates by “borrowing strength” across response variables in a model-based, small area estimation approach.

These types of variables are widely present in social surveys, where continuous variables are strongly related to binary variables. Here, I considered income and poverty indicators; however, future work will investigate other modeling scenarios, that is, variables measured on other scales, such as models for count data. Indeed, the use of multivariate mixed (multilevel) models is of particular relevance in sociology, where multivariate research questions can be formulated and answered by interpreting model parameter estimates, in particular the association parameters. These association parameters, such as correlations, can be estimated via the use of multivariate models. For example, the literature has shown a strong relationship between happiness and health (see Argyle 1997; Graham 2012). If one considers a continuous latent variable measuring happiness, and a binary variable measuring a health domain, the model object of this study can be applied to investigate relationships in the case of multilevel data. Pierewan and Tampubolon (2015) used a multivariate mixed model to investigate the relationships between happiness and health. The model I proposed in this small area estimation problem can be used to investigate subnational differences of such indicators, and to produce more efficient estimates than in the scenario where two separate multilevel models are estimated.

In the context of income and expenditures analysis, researchers might be interested in studying relationships between deprivation and expenditures; these show correlations (see Moretti and Shlomo 2023) that can be taken into account in the modeling stage. Income is also related to food poverty. Thus, a binary variable indicating whether a household is food deprived or at risk for food deprivation can be built, and this is strongly related to household income (for the link between poverty and food insecurity, see Wight et al. 2014). Therefore, similar to the approach followed in the simulation study, the method proposed here is appropriate to provide geographic understanding around food deprivation and monetary poverty dimensions.

It is important to acknowledge that the mixed models adopted in this article can be extended to a spatial dependence setting (Chandra, Chambers, and Salvati 2019; Chandra et al. 2018). Here, one strategy to incorporate spatial information into the models is to extend the random effects model allowing for spatially correlated random area effects. For instance, this can be achieved by using a simultaneous autoregressive model (Anselin 1988; Cressie 2015). As highlighted by Chandra et al. (2018), in practical applications, it is often reasonable to consider that the effects of neighboring geographic areas are correlated (defined by a contiguity matrix that is included in the modeling stage). Readers interested in simultaneous autoregressive models in small area estimation can refer to the following literature, which presents extended evaluations and applications: Singh, Shukla, and Kundu (2005), Pratesi and Salvati (2008), Pratesi and Salvati (2009), Molina (2009), Marhuenda, Molina, and Morales (2013), and Porter, Wikle, and Holan (2015). In addition, Bayesian approaches to small area estimation considering conditionally autoregressive models can model the spatial dependence (see Besag, York, and Mollié 1991; Leroux, Lei, and Breslow 2000; Mercer et al. 2014). These spatial extensions are widely used for small area estimation problems of health and social indicators in sociology, demography, and epidemiology. Dwyer-Lindgren et al. (2016) used the spatial Besag-York-Mollie model (Besag et al. 1991) to estimate mortality rates. The Centers for Disease Control and Prevention PLACES project adopted a similar strategy to estimate small area indicators of morbidity (https://www.cdc.gov/places/methodology/index.html). Regarding poverty indicators, readers may refer to Marhuenda et al. (2013) and Giusti, Masserini, and Pratesi (2017). In this article, I did not include a spatial process into the random area effects. Spatial extensions of the proposed methods are interesting and useful for users, and they will be the subject of future research given they require particular methodological attention.