Small area estimation of proportions under area-level compositional mixed models

This paper introduces area-level compositional mixed models by applying transformations to a multivariate Fay–Herriot model. Small area estimators of the proportions of the categories of a classification variable are derived from the new model, and the corresponding mean squared errors are estimated by parametric bootstrap. Several simulation experiments designed to analyse the behaviour of the introduced estimators are carried out. An application to real data from the Spanish Labour Force Survey of Galicia (north-west of Spain), in the first quarter of 2017, is given. The target is the estimation of domain proportions of people in the four categories of the variable labour status: under 16 years, employed, unemployed and inactive.


Introduction
Statistical offices are interested on the estimation of socio-economic indicators, like proportions or counts, for the whole population or for subsets called domains. Sampling designs are developed for obtaining precise estimators on their target (planned) domains. Statisticians are also asked to provide estimates for unplanned domains (small areas), where sample sizes are too small to carry out such estimations. Small Area Estimation (SAE) deals with this kind of problems by combining tools of survey sampling and statistical modelling at the unit or at the area level. The monographs of Rao (2003) and of Rao and Molina (2015) give a general description of SAE.
In Galicia (north-west of Spain), the Spanish Labour Force Survey (SLFS) provides information about labour market indicators. The territory of Galicia is hierarchically divided into counties and municipalities. As the sampling design of the SLFS is stratified with strata defined by the size of the municipalities and most municipalities are not represented in the sample, the direct estimates at the municipal or county level have a low accuracy. In this context, estimating labour force indicators is thus a SAE problem. The objective of this paper is to estimate proportions of under 16 years, employed, unemployed and inactive people in counties of Galicia crossed by sex under an area-level model-based approach.
Under the unit-level approach, Chambers et al. (2016), Hobza and Morales (2016) and Hobza et al. (2018) have derived estimators of domain proportions and counts based on M-quantile or binomial logit models for binary outcomes. Under the arealevel approach, Chambers et al. (2014), Dreassi et al. (2014), Tzavidis et al. (2015) and Boubeta et al. (2016Boubeta et al. ( , 2017 applied binomial, negative binomial or Poisson regression models for estimating the domain counts or proportions. Esteban et al. (2012), Marhuenda et al. (2013Marhuenda et al. ( , 2014 and Morales et al. (2015) derived estimators of proportions based on area-level linear mixed models. Concerning the Bayesian approach to small area estimation of proportions, the contributions of Farrell (2000), Larsen (2003), Chen and Lahiri (2012) and Liu and Lahiri (2017) are relevant. These approaches are based on univariate models that do not consider the possibility of jointly estimating the counts or proportions of more than two categories.
In labour force statistics, some indicators of interest are the totals or proportions of the categories of a classification variable. This is to say, they are domain compositional parameters summing up to one or to a known integer number. The seminal book of Aitchison (1986), the more recent book edited by Pawlowsky-Glahn and Buccianti (2011) and the papers Egozcue et al. (2003) and Egozcue and Pawlowsky-Glahn (2019) are basic references for an introduction to compositional data analysis.
The estimation of compositional parameters requires using multivariate models. At the unit level, Scealy and Welsh (2017) have applied a directional mixed effects model for predicting the proportions of total weekly expenditure on food and housing costs for households in a chosen set of domains. At the area level, we can find methodologies for contingency tables where the cell counts are explained by categorical auxiliary variables or by regression-based inference procedures allowing for continuous auxiliary variables. Within the first approach, Zhang and Chambers (2004) developed a class of log-linear structural models that is suited to estimation of small area crossclassified counts based on survey data. Berg and Fuller (2014) gave a SAE procedure for a two-way table of proportions with predictors based on a nonlinear mixed model. Under the regression-based setup, Ferrante and Trivisano (2010) proposed a multivariate SAE approach for count data based on the multivariate Poisson-lognormal distribution an derived hierarchical Bayes predictors. Souza and Moura (2016) and Fabrizi et al. (2016) deal with multivariate Beta regression models in SAE. Saei and Chambers (2003), Molina et al. (2007) and López-Vizcaíno et al. (2013 have applied multinomial logit mixed models to category counts for estimating domain totals of labour status categories. This paper introduces three area-level compositional mixed models to obtain small area estimates of labour force proportions in Galicia. The first model is an additive logistic transformation of a multivariate Fay-Herriot (MFH) model with no restrictions on the covariance matrix of the vector of random effects. The second and third models are defined similarly by using centred and isometric logratio transformations, respectively. The new models give higher flexibility than the multinomial model for dealing with the covariance structure of the categories. The MFH model was suggested by Fay (1987) and studied by Datta et al. (1991Datta et al. ( , 1996, Ghosh et al. (1996), González-Manteiga et al. (2008a), Benavent and Morales (2016) and Arima et al. (2017).
We propose a trivariate Fay-Herriot (TFH) model for analyzing the SLFS data, where the vector of random effects (model errors) has an unstructured covariance matrix with unknown components and the vector of sampling errors has a known covariance matrix. As far as we know, this model is studied and applied for the first time to SAE problems. We do not implement a further parametric modelling of the last covariance matrix as Berg and Fuller (2012) do. The estimates of the TFH model parameters are obtained by using the residual maximum likelihood (REML) estimation method. The fitted model is then used to estimate the proportions of under 16 years, employed, unemployed and inactive people in Galician counties.
The estimation of the mean squared error (MSE) of a model-based predictor is an important issue that has no easy solution. Under nonlinear models, the problem is even more difficult. In this paper, we follow the resampling approach appearing in González-Manteiga et al. (2007, 2008a to introduce a parametric bootstrap procedure. Further research would be needed to approximate the MSEs of predictors as it is done by Slud and Maiti (2006) in transformed univariate Fay-Herriot models.
This paper introduces statistical methodology that is new in three main aspects: (1) the employment of three transformations of area-level compositional survey data, (2) the use of TFH models (as a particular case of the general multivariate case) with unstructured covariance matrix for modelling the transformed data and capturing the sample correlations and (3) the derivation of domain-level predictors of proportions and counts based on the TFH model fitted to the transformed data.
The remainder of the paper is organized as follows. Section 2 gives an introduction to the labour force data and to the SAE problem of interest. Section 3 introduces the additive, centred and isometric logratio area-level model-based approaches for estimating domain compositional parameters. Section 4 describes the considered TFH model. Section 5 develops the proposed compositional predictors and the corresponding MSE estimation procedures. Section 6 applies the proposed methodology to data from the SLFS of the first quarter of 2017 in Galicia. Section 7 gives some conclusions.
The paper contains four appendixes in a supplementary material file. Appendix A presents four simulation experiments. The target of Simulation 1 is to check the behaviour of the REML algorithm for fitting the TFH model. Simulation 2 investigates the performance of the compositional and the multinomial predictors of the category proportions. Simulation 3 analyses the parametric bootstrap estimator of the MSEs. Simulation 4 studies the behaviour of the predictors of proportions when the target data are generated by a multinomial mixed model. Appendix B reviews the competitor predictors based on multinomial mixed models. Appendix C gives the REML Fisher scoring algorithm and derives the REML score vector and Fisher information matrix. Appendix D describes the centred and the isometric logratio transformations of compositions.

The problem of interest
The SLFS in Galicia is a quarterly survey following a stratified two-stage random sampling design. Primary sampling units are composed by census sections, which are geographical areas with a minimum of 500 dwellings or about 3000 people. Secondary sampling units are composed by main family dwellings and permanent accommodations. Subsampling is not carried out in secondary sampling units, and data are collected on all persons who regularly live in the same dwelling. The SLFS gets information about the labour market.
Galicia is divided into four provinces. They are Coruña, Lugo, Ourense and Pontevedra, coded by the Spanish Statistical Office as 15, 27, 32 and 36, respectively. Each province is hierarchically partitioned in comarcas (counties) and municipalities. Our domains of interest are the counties crossed by sex. As there are 51 counties in the SLFS of Galicia, we have 102 domains. Our goal is to estimate domain proportions of people in the categories of the variable "labour status" by using SLFS data from the first quarter of 2017.
The SLFS is designed to obtain precise direct estimates of labour indicators at the province level. For the considered SLFS data, the minimum domain sample size is 10, the first quartile is 43 and the median is 79. Therefore, obtaining reliable estimates for target domains is a small area estimation problem and borrowing strength from auxiliary data is recommended.
In mathematical terms, Galicia is a population U = ∪ D d=1 U d partitioned in D domains U d . Each domain is partitioned in subsets U dk , k = 1, . . . , q, defined by the classification variable "labour force status" that classifies units into a finite number of categories. The q = 4 categories are ≤ 15 years (k = 1), employed (k = 2), unemployed (k = 3) and inactive (k = 4). Let N and N d be the sizes of U and U d , respectively. Consider the study variables taking the values z dk j = 1 if the unit j from the domain U d is in the category k and z dk j = 0 otherwise. The target parameters are the domain means (proportions) and the domain totals (counts), i.e.
(2.1) AsZ d1 +· · ·+Z dq = 1 and Z d1 +· · ·+ Z dq = N d , with N d known, we are interested in estimating domain compositions with q = 4 categories. For any sample s ⊂ U extracted from the population, s d denotes the subsample from U d of size n d and w d j represents the SLFS sampling weight of the unit j from the subsample s d . The sample proportions and counts are as: Direct estimators ofZ dk and Z dk are as: As the SLFS sampling weights w d j are derived from the inverses of the inclusion probabilities after non-response correction and calibration at the province and at the regional level, the direct estimators ( Appendix B describes these MLM models for a response vector ξ d = (ξ d1 , . . . , ξ dq−1 ) and the corresponding predictors of the multinomial category probabilities p dk . In both cases, ξ dk = z dk. or ξ dk =Ẑ dir dk , the covariances, given in (2.2) of Appendix B, under MLM models are negative. Further, under the multinomial distribution, it holds that cov M (ξ dk 1 , ξ dk 2 ) = 0 if and only if p dk 1 = 0 or p dk 2 = 0, which implies that cov M (ξ dk 1 , ξ dk ) = 0 ∀k = k 1 or cov M (ξ dk 2 , ξ dk ) = 0 ∀k = k 2 . Therefore, the multinomial correlation structure is rather rigid. More concretely, if the number of categories is q = 4, then the six variance components (three variances and three covariances) depend only on three category parameters through the formulas (2.2) of Appendix B. In practice, the multinomial covariance structure does not necessarily fit to the true covariances cov(ξ dk 1 . , ξ dk 2 . ), k 1 = k 2 , k 1 , k 2 = 1, . . . , q, of sample counts or direct estimators of totals.
For the categories k 1 , k 2 = 1, 2, 3, 4, Table 1 presents the correlations of the set of values {(Ẑ dir dk 1 ,Ẑ dir dk 2 ) : d = 1, . . . , D} (left) and {(z dk 1 . ,z dk 2 . ) : d = 1, . . . , D} (right). These correlations are calculated from the direct estimates (2.3) and the simple estimates (2.2) of the category proportions along the domains. As they are calculated with aggregated data at the domain level, we call them domain-level empirical correlations between categories or, in short, domain-level correlations. Table 1 shows that most, but not all, correlations are negative. Nevertheless, we find the positive correlations 0.21 (right) or 0.04 (right) and 0.07 (left) that are close to zero.  The design-based within-domain covariance cov π (Ẑ dir dk 1 ,Ẑ dir dk 2 ), k 1 , k 2 = 1, . . . , q− 1, can be estimated bŷ where the case k 1 = k 2 = k denotes estimated variance, i.e.v ar π (Ẑ dir dk ) = cov π (Ẑ dir dk ,Ẑ dir dk ). The last formulas are obtained from Särndal et al. (1992), pp. 43, 185 and 391, with the simplifications w d j = 1/π d j , π d j,d j = π d j and π di,d j = π di π d j , i = j in the second-order inclusion probabilities. By applying formula (2.4), Fig. 1 (left) plots the estimated design-based covariances s k 1 k 2 =ĉov π (Ẑ dir dk 1 ,Ẑ dir dk 2 ), k 1 , k 2 = 1, 2, 3, 4, k 1 = k 2 , for the direct estimators of all the category proportions. As these covariances are calculated from unit-level data, they are called unit-level covariances. Figure 1 Under the multinomial mixed model, this fact implies (as explained above) that the domain proportions of people in the categories "under 16 years" and "unemployed" should be close to zero, which contradicts the observed sampling proportions. Figure 1 (right) plots the corresponding unit-level correlations r k 1 k 2 . In the application to SLFS data, these unit-level covariances are used to set the error covariance matrix of the fitted TFH model.
The exploratory data analysis shows that the covariance structure of the MLM models M1 and M2, described in Appendix B, cannot take into account simultaneously both types of covariances and variances (unit level and domain level). This is to say, the multinomial mixed model does not fit well to the domain-level correlations appearing in Table 1 or to the within-domain covariances shown in Fig. 1. Further, the definition of the model itself does not allow the introduction of both types of covariances and variances as MFH models do. This is why we propose transforming the direct estimators of domain proportions and fitting the transformed data to a more flexible compositional mixed model.

Transformations of compositions
This section considers three transformation of q-compositions onto R q−1 . In what follows, we use the simpler notation z dk and with inverse the additive logistic (alogist) transformation For i, j = 1, . . . , q − 1, the first partial derivatives of the alogist transformation are as: This to say, under the alogist transformation z i is an increasing function of y i and a decreasing transformation of y j , i, j = 1, . . . , q − 1, j = i. Alternatively, we may consider the centred or the isometric logratio (clr or ilr) transformations described in Appendix D. For ease of exposition, this paper deals mainly with the alr transformation. The mathematical developments for the clr or ilr transformations can be done similarly.
In applications to real data, we can have domain compositions (z d1 , . . . , z dq ) with some zero components z dk = 0. Since the logarithm of zero is −∞, we cannot apply alr, clr and ilr logratio transformations for analyzing compositional data. In our application to real data, we have no zeros. However, in other domain-level data we might find a low number of zeros. This can occur if a domain sample size is very small and none of the sampled units belong to a given category k. A practical solution for this problem is replacing zeros by a small numerical value ε > 0 and making a rounding-off adjustment process. For example, if (z d1 , . . . , z dq ) has m zeros, we can apply the recommendation appearing in Section 11.5 of Aitchison (1986). This is to say, we can substitute the zeros by ε = (m + 1)(q − m)δq −2 and we can subtract m(m + 1)δq −2 to each positive component, where δ is the maximum rounding-off error of the positive z dk 's. We can take var π (z dk ) 1/2 as rounding-off error of z dk .
In any case, applied statisticians should do a sensitivity analysis of the results to the particular method employed to deal with zeros and a model diagnostic study. The presence of structural zeros or a non-negligible amount of sample zeros is a severe issue when using alr, clr and ilr logratio transformations. We are not in favour of using the three introduced transformations in those cases. Other transformations and models for dealing with compositional data could be alternatively applied. For example, the directional mixed effects models applied by Scealy and Welsh (2017) overcome this problem.
The first derivatives of the alr transformation h k are as: Let 1 a and 0 a be the a × 1 vectors with all the components equal to one and to zero, respectively. In matrix form, we have where I a is the a × a identity matrix. Alternatively, we can take z 0 = 1 D D d=1 z d or we can select z 0 depending on d and close to z d for a better Taylor expansion approximation. For the clr and ilr transformations, Appendix D shows that h(z 0 ) = 0 q−1 and gives H (z 0 ). From (3.3), we get the approximated covariance matrix: where μ d is a mean vector depending on unknown regression parameters and auxiliary variables and V d is a covariance depending of some unknown parameters. The next section gives a flexible multivariate area-level linear mixed model for y d , d = 1, . . . , D, which allows positive and negative covariances. Depending on the employed transformation, the resulting models for z d , d = 1, . . . , D, are called alr, clr and ilr compositional mixed models. These models arise as alternative to the MLM models, for z d. or z d , described in Appendix B. For ease of exposition, we present the case q = 4 appearing in the application to real data.

Trivariate Fay-Herriot models
be the matrix of design-based covariance estimators (2.4), which can contain positive and negative covariances. Let y d = (y d1 , y d2 , y d3 ) be the alr transformation of z d .
The TFH model is defined in two stages. The sampling model is is the covariance matrix given in (3.1) and H 0 is defined in (3.4). Alternatively, y d can be defined as the clr or ilr transformation of z d . In that case, the matrix H 0 is taken from Appendix D.
Moreover, it is assumed that the μ dk 's are linearly related to r k explanatory variables associated with the k-th category in the domain d. For k = 1, 2, 3, let x dk = (x dk1 , . . . , x dkr k ) be a row vector containing the r k explanatory variables for μ dk and let X d = diag (x d1 , x d2 , x d3 ) 3×r with r = r 1 +r 2 +r 3 . Let β k be a column vector of size r k containing the regression parameters for μ dk and let β = β 1 , β 2 , β 3 r ×1 . This section introduces a TFH model by assuming (4.1) and the linking model: where the vectors u d 's are independent and independent of the vectors e d 's. Unlike the MFH models studied by Benavent and Morales (2016), the 3 × 3 covariance matrices V ud are unstructured and depend on six unknown parameters, The matrix V ud explains the covariance structure of the alr transformations of the direct estimators y d that is not taken into account by the sampling errors e d or by the auxiliary variables X d . Let I n be the n × n identity matrix and δ d be the Kronecker delta, and define where col and col are matrix operators stacking by columns and rows, respectively. In matrix form, the TFH model (4.1)+(4.2) is Under model (4.3), it holds that Further, the best linear unbiased estimator (BLUE) of β and the best linear unbiased predictors (BLUP) of u and μ are as: (4.4) The residual maximum likelihood (REML) method maximizes the joint probability density function of a vector of 3D − r independent contrasts ω = W y, where W is a 3D ×(3D −r ) matrix with linearly independent columns and such that W W = I 3D−r and W X = 0. It holds that ω is independent of the BLUEβ B given in (4.4). The joint probability density function of ω is the REML likelihood. The REML log-likelihood of model (4.3) is and P X = 0. Appendix C gives the calculation of the score vector S(θ ) = (S 1 , . . . , S 6 ) and the Fisher information matrix F(θ ) = F a,b a,b=1,...,6 , where Appendix C also gives the REML Fisher scoring algorithm for fitting the TFH model (4.3). To initiate this algorithm, a possible set of starting values isθ 4,0 =θ 5,0 =θ 6,0 = 0,θ k,0 =σ 2 uk,0 , k = 1, 2, 3, whereσ 2 uk,0 is the REML or the ML or the Prasad and Rao (1990) moment-based estimator of σ 2 uk in the k-th marginal Fay-Herriot model. They can be calculated with the R library sae. Another possibility is to take V (0) ud as the domain-level variance matrix of the alr transformations of the direct estimators of the category proportions. This is the approach described and implemented around Table 2 in the application to real data.
The output of REML Fisher scoring algorithm,θ , is the REML estimator of θ . By pluggingθ in V u , we getV u = V u (θ) andV =V u + V e . By substitutingV u in (4.4), we obtain the EBLUP of μ = X β + Zu, i.e.

Compositional predictors of proportions and counts
In practice, assumption (1) rarely holds. The direct estimator z dk is not calculated with the inverse of the inclusion probabilities, but with sampling weights (expansion factors) that are not calibrated to domain totals. Therefore, z dk is, in general, biased with respect to the sampling design distribution. Similarly,z dk. is also biased for estimatingZ dk , as it is calculated with weights equal to one. Therefore, predictors based on both models, FH and MLM, have problems to fulfil assumption (1) when dealing with real data.
Assumption (2) might be accepted if the area-level model has a "good" fit to data, which may happen if the set of auxiliary variables is highly correlated with the dependent variable. Further, by analogy to the GREG estimator, the auxiliary variables produce a calibration effect to their domain totals and a reduction of the design-based bias. In many cases, the assumptions (1) and (2) "approximately" hold and the empirical best predictor (EBP) will have some small bias, as the direct estimator has, but lower variance. So it is worthwhile to calculate EBPs based on area-level models, as they will tend to have lower design-based MSEs than direct estimators. This paper introduces a compositional model for z d that is introduced by assuming a TFH model on the alr transformed vector y d . In terms of y d , the compositional model does not fulfil assumption (1), because the direct estimator y d is not design-based unbiased for the alr transformation ofZ dk , i.e. E π [y dk ] = log(Z dk /Z d4 ), k = 1, 2, 3. For fulfilling the assumption (2), the FMEs to be predicted under the TFM model should be For the clr and ilr transformations, the alternative vectors are p clr Cd = clr −1 (μ Cd ) and p ilr Cd = ilr −1 (μ Cd ), respectively, where μ Cd = (μ Cd1 , μ Cd2 , μ Cd3 ). By predicting p Cdk instead of E C [z dk |u d ], we gain computational efficiency, but we have problems with assumption (2), asZ dk might not be considered as the realization of p Cdk ; in short, p Cdk = E π [z dk. ] ≈Z dk . However, the disadvantage of not fulfilling the assumptions (1) and (2) can be compensated, as in practice happens with the FH and MLM models, by the use of auxiliary variables correlated with the objective variable and with the modelling of the covariance structure of the categories. So that what is lost on the one hand is recovered on the other.
For the clr and ilr transformations, the EBP of p Cdk is obtained by substituting the k-th component of the alogist transformation: by the corresponding k-th component of clr −1 and ilr −1 transformations, respectively. The estimation of compositions, like domain proportions of categories of a classification variable, requires the selection of the last (q-th) category. The last category is not a control or reference category as it happens in some ANOVA-type statistical analyses. The introduced methodology is not invariant with respect to the selection of the category q, but provides together the predictors of the proportions of the q categories. In practice, the selection of the first q − 1 categories should be based on the available explanatory variables. A good approach is to select as target categories (k = 1, . . . , q − 1) those ones that can be better explained by the auxiliary variables.
For every domain d, the predictors based on multinomial or compositional models fulfil the two conditions: (1) 0 ≤ p dk ≤ 1, k = 1, . . . , q, and (2) q k=1 p dk = 1. Estimating p dk with predictors based on univariate area-level mixed models, like Fay-Herriot or binomial logit, is not a good option because the conditions (1) and (2) might not be fulfilled.

Application to Labour Force Survey data
This section gives an application of the alr compositional mixed model to the SLFS data described in Sect. 2. The TFH model is fitted to the target data and to a set of significant auxiliary aggregated variables taken from the administrative registers: PMH containing the official demographic data at municipal level, SSoc containing data from the Social Security System and SPEG containing data of employment claimants. The considered domain-level auxiliary variables are as: • SS: proportion of population registered in SSoc.
• REG: proportion of population registered as unemployed in SPEG.
The target parameters are the proportions of the four categories of the variable labour status, i.e. ≤ 15 years, employed, unemployed and inactive people per sex in counties of Galicia. We are interested in estimating domain compositions with q = 4 categories. As the explanatory variables, to15, SS and REG, are highly correlated with ≤ 15 years, employed and unemployed, respectively, for fitting compositional or multinomial mixed models to the SLFS data, we number the labour status categories from 1 to 4, so that inactive is the fourth one (q-th category). We denote the alr transformations of the direct estimators of the category proportions by y dk = log(z dk /z dq ), k = 1, 2, 3. For the sake of brevity, we do not present the data analyses with the clr or ilr transformations. Keeping these assumptions in mind, below we present the results of the application to SLFS data. Table 2 presents the domain-level correlations calculated for the sets {(y dk 1 , y dk 2 ) : d = 1, . . . , D}, k 1 , k 2 = 1, 2. Because of the alr transformation, the domain-level correlation patterns of the y-variables are different from the corresponding ones of the z-variables given in Table 1. The correlations in Table 2 are used as seeds for the matrices V ud in the Fisher scoring algorithm that calculates the REML estimators of the selected TFH model.   Table 3 presents domain-level correlations for the response variables y dk of the TFH model and the auxiliary variables x dki . For the k-th category and the i-th auxiliary variable, the correlations are calculated for the set of values {(y dk , x dki ) : d = 1, . . . , D}. Reading this table by rows, it is observed that to15 is the most positively correlated variable for y 1 and 25to54 and SS are the most positively correlated variables for y 2 . Concerning y 3 , the first five auxiliary variables have a similar positive correlation. As expected, 55to is negatively correlated with y 1 , y 2 and y 3 .
A set of appropriate auxiliary variables is selected, and the corresponding TFH model is fitted to the data (y dk , x dki ), d = 1, . . . , D, k = 1, 2, 3, i = 1, . . . , r k . For the alr transformation, Table 4 presents the estimates of the regression parameters for the TFH model and their estimated standard deviations. It also presents the p-values, defined in (4.8), for testing the hypothesis H 0 : β ki = 0, k = 1, 2, 3, i = 1, . . . , r i . Intercept parameters for y 1 , y 2 , y 3 are denoted by c 1 , c 2 , c 3 , respectively. Appendix D presents the corresponding tables for the transformations clr and ilr. Based on those tables and on further model diagnostics, we select the alr transformation.
As all the auxiliary variables are proportions, the sign and the magnitude of the regression parameters give interesting interpretations. The alr transformation of the direct estimator of the SLFS proportion of people under 16 years, y d1 , is solely explained with positive sign by the corresponding PMH proportion to15. For the proportion of employed people, y d2 tends to be greater in those domains with larger proportion of people registered in SSoc and greater proportion of working age people. The proportions of people SS and 25to54 are more relevant than 16to24 for predicting y d2 . For the proportion of unemployed people, y d3 tends to be greater in those domains with larger proportion of people registered as unemployed in SPEG and greater proportion of working age people. As expected, REG is more important for predicting y d3 than 16to24 and 25to54. We have further fitted a second TFH model with 55to in the place of 16to24 and 25to54. The second model gives similar predictions because the auxiliary variables to15, 16to24, 25to54 and 55to sum up to one, and therefore,   Table 5 gives the 95% confidence intervals (CIs) of the variance components, defined in (4.7), for θ 1 , . . . , θ 6 . We observe that the variances σ 2 u1 , σ 2 u2 , σ 2 u3 and the correlation ρ 12 are significantly greater than zero. Figure 2 plots the dispersion graphs of the target variables y d1 (left), y d2 (centre) and y d3 versus their auxiliary variables with larger regression parameter, i.e. to15, SS and REG, respectively. Accordingly with Table 3, we observe the linear patterns related to positive high linear correlations. Figures 3 and 4 plot the marginal residuals and standardized marginal residuals of the fitted TFH model. The residuals are rather symmetric around zero and do not present any relevant pattern. Further, there are few standardized residuals (7 among 306) outside the interval (−3, 3). Figure 5 plots the compositional versus the direct estimates of the proportions of people under 16 years (left), employed (centre) and unemployed (right). We observe that model-based estimators take values rather symmetrically around the direct estimates. As the direct estimators of proportions are basically design-based unbiased, Fig. 5 suggests that compositional estimators partially share this property. In addition to the compositional plug-in predictors based on the selected compositional mixed model, we also calculate the direct estimators and the multinomial predictors. Taking into account the recommendations of López-Vizcaíno et al. (2013) and the last comment of Appendix B, the MLM model is fitted to the vectors of sample counts z d. = (z d1. , z d2. , z d3. ) with sample sizes n d , d = 1, . . . , D. For the sake of comparability, we use the same auxiliary variables as the fitted compositional model. This is to say, the auxiliary variables are listed in Table 4. Figure 6 plots the direct, multinomial and compositional plug-in estimated proportions of people under 16 years (left), employed (centre) and unemployed (right) for domains sorted by sample size. For the three categories, the compositional plug-in predictors present the smoothest behaviour across domains. Figure 7 plots the design-based estimates of the root mean squared errors (RMSE) of the direct estimators (D) and the parametric bootstrap estimates of the RMSEs of the multinomial and the compositional plug-in predictors. For the sake of comparability, the RMSEs of the multinomial (MC) and the compositional (CC) predictors are calculated under the assumption that the distribution of the fitted TFH Fay-Herriot model is the true one. Therefore, we run the bootstrap procedure by generating data from the fitted TFH model. The RMSEs of the corresponding direct estimates are estimated by applying formula (2.4). Nevertheless, we know that in real life there are no true models, but useful models. Therefore, we also calculate parametric bootstrap estimates of the RMSEs of the multinomial predictors assuming that the multinomial model is valid. The new estimated RMSEs (MM) can be interpreted as a measure of the goodness of fit of the multinomial model to the data. Figure 7 shows that the compositional plug-in predictors have the best results in terms of RMSE. Figures 8 and 9 map the compositional plug-in estimated county proportions of men (left) and women (right) under 16 years old and inactive, respectively. The colours are darker in areas with higher proportions. We observe that the counties with the youngest people are in the west coast, in the provinces of Coruña and Pontevedra where it is the Atlantic motorway (from Vigo to Coruña). On the contrary, it can be observed that the counties that are in the east of Galicia (in the provinces of Lugo and Ourense), in general terms, have a lot of inactive population, more than the 50% of their population. The Costa da Morte counties (in the west) are also in this situation. Figures 10 and 11 plot the compositional plug-in estimated county proportions of men (left) and women (right) that are employed and unemployed, respectively. <=0.08 (14) >0.08 <= 0.12 (13) >0.12 <= 0.15 (13) >0.15 (13) <=0.09 (14) >0.09 <= 0.11 (13) >0.11 <= 0.13 (13) >0.13 ( <=0.05 (14) >0.05 <= 0.07 (13) >0.07 <= 0.09 (13) >0.09 (13) Proportion of unemployed people. Women <=0.04 (14) >0.04 <= 0.07 (13) >0.07 <= 0.08 (13) >0.08 (13) Fig. 11 Estimated county proportions of unemployed men (left) and women (right)  Concerning employment, we observe that the touristic counties around the Santiago trail have the largest proportions. On the contrary, the industrial areas around Vigo (south-west) and the agricultural areas of Ourense (south-east) present the largest proportions of unemployed people. Figure 9 shows that the quantiles of the distribution of the proportion of inactive women are displaced to the right with respect to that of men. Therefore, the proportions  of inactive women tend to be greater than those of men in the regions of Galicia. Figures 10 and 11 give the opposite conclusion for employed and unemployed people, respectively. Tables 6, 7 and 8 present some condensed numerical results for men and women, respectively. The tables have been constructed in two steps. The domains are sorted by province. Within each province, the domains are sorted by sample size, starting by the domain with the smallest sample size. A selection of 18 domains out of 51 is done from the positions 1, 4, 7, . . . , 49 and 51. The tables give the direct and the compositional plug-in estimates (labelled by "dir" and "mod", respectively) and the corresponding RMSE estimates. The provinces are labelled by "prov", with codes given in Sect. 2, and the sample sizes are denoted by n.

Proportion of people under 16 years. Women
Tables 6 and 7 are partitioned in three vertical sections dealing with the estimation of proportion of people under 16 years, employed and unemployed people. Table 8 contains the estimated proportions of inactive men and women and the corresponding RMSE estimates. By observing the columns of RMSEs, we conclude that the compositional plug-in predictors are preferred to the direct estimators.
We recall that the compositional predictors of the four categories are derived from the fitted TFH model by applying the alogist transformation, so they are jointly calculated. The last category (inactive) is not a reference category. Their results appear in Table 8 because of the lack of space when constructing Tables 6 and 7.

Conclusions
This paper introduces predictors of category proportions based on an area-level compositional mixed model. A TFH model is introduced for modelling the additive logratio transformations of the direct estimators of the category proportions. The additive logistic transformations of the EBLUPs under the TFH model are the proposed compositional plug-in predictors. Similarly, the centred or the isometric logratio transformation (see Appendix D) can be employed. In addition, the compositional empirical best predictors are also introduced and empirically investigated. The first predictor is easy to calculate, but the second one requires the approximation of integrals in R 3 . This numerical problem demands rather high computational time to achieve an acceptable precision. Because of this issue and the results of Simulation 2, the use of the compositional plug-in predictor is recommended.
Unlike the multinomial approach of López-Vizcaíno et al. (2013), the compositional predictors take into account the sampling design. The sampling weights are employed by means of the direct estimators of the category proportions and through the estimators of their design-based variances and covariances.
As discussed in Section 4, both models (MLM and compositional) have problems to predict the population parameterZ dk . However, these problems are compensated by the use of auxiliary variables correlated with the objective variable and with the modelling of the covariance structure of the categories. It is at this point where the compositional mixed model presents its greatest advantages compared to the MLM model. The compositional approach is more flexible and can be adapted to different correlation structures, while multinomial correlation structure is rather rigid (correlations negative) and does not fit well to both types of covariances and variances (unit level and domain level). The results of four simulation experiments and the analysis of real data confirm this assertion.
The new small area estimation methodology is applied to Labour Force Survey data from Galicia, a region in north-west of Spain, in the period January-April 2017. The selected compositional mixed model has had a better fit to the data than the corresponding multinomial mixed model, and therefore, the labour status proportions per county and sex are finally estimated by using the compositional plug-in predictors with its mean squares errors calculated by parametric bootstrap.
As for the labour market results in Galicia, we can conclude that the west coast, in general terms, is the most dynamic part with a higher proportion of employed people and people under 16 years old. There is a big problem in the south-east because there are several counties with a proportion of inactive people over the 50%. This area is essentially rural. Fixing the population in rural areas, with a decent living and income levels, is one basic requirement to ensure territorial balance and environmental sustainability. Galicia is ageing and more accentuated in rural areas, due to the emigration of young people and the fall in the birth rate. This reality, with obvious economic and sociological consequences, demands answers from social and employment policies.
Although the introduced methodology is applied to the variable "labour status" with q = 4 categories, it can be extended to categorical variables with any number q ≥ 2 of categories. However, we only recommend their use for the cases q = 2, 3, 4, 5. This paper presents the mathematical developments for the case q = 4 leading to a TFH model with six variances or covariance parameters. The cases q = 2 and q = 3 are more simpler as they require fitting univariate and bivariate Fay-Herriot models with one and three variance components, respectively. The case q = 5 yields to a MFH model with ten variance components, where the fitting algorithm might have convergence problems when the number of domains D is not big enough.
The application to real data can be carried out in the subpopulation of people aged 16 or more, where the number of categories is q = 3 (employed, unemployed and inactive). This simpler scenario requires fitting a bivariate Fay-Herriot model to the transformed survey data and obtaining the predictions of the domain proportions of the three considered categories in a similar way as it is done in the case of q = 4 categories. We have decided to follow the more complex TFH approach because of the strength of the available auxiliary variables and to guarantee the coherence of the estimates of the four category proportions. For the categories ≤ 15, employed and unemployed, we have the highly correlated variables to15, SS and REG. This last fact also motivated the selection of inactive as the fourth (reference) class.
Compositional data play an important role in public statistics. In this case, we applied the introduced methodology to the SLFS, but it is useful in other topics of the official statistics, like the classification of the population by the educational level, the income level or the type of household expenditure. In all these situations, it is necessary to take into account the simplex constraints. This work may be extended to incorporate spatial correlations between domains. In practice, it is often reasonable to assume that the effects associated with neighbouring areas are proportionally correlated with a measure of distance. Also, it is important to know the evolution of the labour market, along the quarters for the counties, in an accurate and stable form, suitable for being used in statistical offices. Therefore, extensions of the introduced methodology to models incorporating temporal correlations are also a future research task.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.