Sample allocation for efficient model-based small area estimation
Section 2. Allocations which utilize the model

• Previous
• Next

2.1 Choosing the model

Pfeffermann (2013) presents a wide variety of models and methods for small area estimation. Our model is one of this assortment, a unit-level mixed model

$$y_{dk}=x_{dk}^{{}^{\prime }}\beta +v_d+e_{dk};k=1,\dots ,N_d;d=1,\dots ,D,(2.1)$$

where $v_d’\text{s}$ are random area effects with mean zero and variance ${\sigma }_v^2$ and $e_{dk}’\text{s}$ are random effects with mean zero and variance ${\sigma }_e^2.$ Furthermore, $E\left(y_{dk}\right)=x_{dk}^{{}^{\prime }}\beta$ and $V\left(y_{dk}\right)={\sigma }_v^2+{\sigma }_e^2$ (total variance). Matrix $V$ is the variance-covariance matrix of the study variable $y.$ This model can be used when unit-level values are available for the auxiliary variables $x.$ We use one auxiliary variable in our study.

Two important measures are needed in developing one of these types of allocations. The first one is a common intra-area correlation $\rho$ and the second one is the ratio $\delta$ between variance components. They are defined as follows:

$$\rho ={\sigma }_v^2/\left({\sigma }_v^2+{\sigma }_e^2\right)\text{and}\delta ={\sigma }_e^2/{\sigma }_v^2=1/\rho -1.(2.2)$$

Before estimating area parameters, the variance components, regression coefficients and area effects must be estimated from the sample data. The BLUE estimator (Best Linear Unbiased Estimator) of $\beta ,$ noted $\tilde{\beta },$ is obtained according to the theory of the general linear model, and it is replaced with its EBLUP estimate $\widehat{\beta }.$

The EBLUP estimate (predicted value) for the area total $Y_d$ of the study variable is the sum of the observed $y-$ values and predicted $y-$ values for units outside the sample:

$${\widehat{Y}}_{d,\text{Eblup}}={\displaystyle \sum _{k\in s_d}{y}_{dk}}+{\displaystyle \sum _{k\in {\overline{s}}_d}{\widehat{y}}_{dk}=}{\displaystyle \sum _{k\in s_d}{y}_{dk}}+{\displaystyle \sum _{k\in {\overline{s}}_d}{x}_{dk}^{{}^{\prime }}\widehat{\beta }}+\left(N_d-n_d\right){\widehat{v}}_d.(2.3)$$

We use the Prasad-Rao approximation (See Rao 2003) of MSE (Mean Squared Error) for finite populations:

$$\text{mse}\left({\widehat{Y}}_{d,\text{Eblup}}\right)=g_{1d}\left({\widehat{\sigma }}_v^2,{\widehat{\sigma }}_e^2\right)+g_{2d}\left({\widehat{\sigma }}_v^2,{\widehat{\sigma }}_e^2\right)+2g_{3d}\left({\widehat{\sigma }}_v^2,{\widehat{\sigma }}_e^2\right)+g_{4d}\left({\widehat{\sigma }}_v^2,{\widehat{\sigma }}_e^2\right),(2.4)$$

where the four components $g_{1d},$ $g_{2d},$ $g_{3d}$ and $g_{4d}$ are defined as follows:

$$\begin{array}{ll}{g}_{1d}\left({\widehat{\sigma }}_v^2,{\widehat{\sigma }}_e^2\right)\hfill & ={\left(N_d-n_d^{*}\right)}^2\left(1-{\widehat{\gamma }}_d\right){\widehat{\sigma }}_v^2,\hfill \\ g_{2d}\left({\widehat{\sigma }}_v^2,{\widehat{\sigma }}_e^2\right)\hfill & ={\left(N_d-n_d^{*}\right)}^2{\left({\overline{x}}_d^{*}-{\widehat{\gamma }}_d{\overline{x}}_d\right)}^{\prime }{\left(X^{\prime }{V}^{-1}X\right)}^{-1}\left({\overline{x}}_d^{*}-{\widehat{\gamma }}_d{\overline{x}}_d\right),\hfill \\ g_{3d}\left({\widehat{\sigma }}_v^2,{\widehat{\sigma }}_e^2\right)\hfill & ={\left(N_d-n_d^{*}\right)}^2{\left(n_d^{*}\right)}^{-2}{\left({\widehat{\sigma }}_v^2+{\widehat{\sigma }}_e^2{\left(n_d^{*}\right)}^{-1}\right)}^{-3}[{\widehat{\sigma }}_e^{4}V\left({\widehat{\sigma }}_v^2\right)\hfill \\ \hfill & +{\widehat{\sigma }}_v^{4}V\left({\widehat{\sigma }}_e^2\right)-2{\widehat{\sigma }}_e^2{\widehat{\sigma }}_v^2\text{Cov}\left({\widehat{\sigma }}_e^2,{\widehat{\sigma }}_v^2\right)],\hfill \\ g_{4d}\left({\widehat{\sigma }}_v^2,{\widehat{\sigma }}_e^2\right)\hfill & =\left(N_d-n_d^{*}\right){\widehat{\sigma }}_e^2.(2.5)\hfill \end{array}$$

The area sample sizes $n_d^{*}$ depend on the sample and are not fixed. The component $g_{1d}$ contains the area-specific ratio ${\widehat{\gamma }}_d={\widehat{\sigma }}_v^2/\left({\widehat{\sigma }}_v^2+{\widehat{\sigma }}_e^2/n_d^{*}\right).$ According to Nissinen (2009, page 53), the $g_{1d}$ component (later simply $g\text{1})$ contributes generally over 90% of the estimated MSE. This component represents uncertainty as regards the variation between the areas. Of course this variation must be strong enough so that such a high proportion for $g\text{1}$ exists.

Unfortunately, the idea of an analytical solution, which means minimizing the sum of MSE’s over areas subject to $n={\displaystyle {\sum }_{d=1}^D{n}_d},$ is difficult and laborious to accomplish because components of the MSE approximation (2.5) include sample information which is unknown, and some components contain complex matrix and variance-covariance operations. We have examined this allocation problem for the first time in an experimental study (Keto and Pahkinen 2009). Now we have developed an allocation based only on the component $g\text{1}$ and auxiliary variable $x.$ The reasoning for this solution is that because $x$ and $y$ are correlated, the between-area variation in $x$ is transferred to $y.$

2.2 Model-based g1 – allocation

The $g1-$ allocation utilizes the auxiliary variable $x$ and the adjusted homogeneity coefficient (Keto and Pahkinen 2014). This coefficient is an approximation of an intra-class correlation (ICC) known of cluster sampling. We regard one area as one cluster. First, simple ANOVA between areas is carried out, and then the adjusted homogeneity measure of variation between the areas can be computed:

$$R_{ax}^2=1-R^2\left(x\right)=1-\text{MSW}/S_x^2,(2.6)$$

where $R^2\left(x\right)$ is the coefficient of determination from regression analysis, MSW (Mean Square Within) is the mean SS (Sum of Squares) of areas and $S_x^2$ is the variance of the auxiliary variable $x.$

Because MSE of the area total is complex, we use only the component $g\text{1,}$ which appears in (2.4) and (2.5), for the reason we have given in Section 2.1. We search for the minimum for the sum of $g\text{1}’s$ over areas:

$${\displaystyle \sum _{d=1}^D{g}_{1d}\left({\sigma }_v^2,{\sigma }_e^2\right)}={\displaystyle \sum _{d=1}^D{\left(N_d-n_d\right)}^2{\left(n_d/{\sigma }_e^2+1/{\sigma }_v^2\right)}^{-1}}(2.7)$$

subject to $n={\displaystyle {\sum }_{d=1}^D{n}_d}.$

We use Lagrange’s multiplier method to find the solution. Therefore, we define the function $F$ of sample sizes $n^{\prime }=\left(n_1,n_2,\dots ,n_D\right)$ and $\lambda :$

$$F\left(n,\lambda \right)={\displaystyle \sum _{d=1}^D{g}_{1d}\left({\sigma }_v^2,{\sigma }_e^2\right)=}{\displaystyle \sum _{d=1}^D{\left(N_d-n_d\right)}^2{\left(n_d/{\sigma }_e^2+1/{\sigma }_v^2\right)}^{-1}}+\lambda \left({\displaystyle \sum _{d=1}^D{n}_d-n}\right).(2.8)$$

We set the derivative of $F$ with respect to the area sample size $n_d$ to zero and solve for $n_d.$ The expression for area sample size $n_d^{g1}$ is as follows:

$$n_d^{g1}=\frac{\left(N_d+\delta \right)\left(n+\delta D\right)}{N+\delta D}-\delta =\frac{N_{d}n-\left(N-N_{d}D-n\right)\left(1/\rho -1\right)}{N+D\left(1/\rho -1\right)},(2.9)$$

where the ratio $\delta$ and the intra-area correlation $\rho$ are defined in (2.2). The only unknown member in (2.9) is the intra-area correlation $\rho .$ Therefore we substitute the known homogeneity measure (2.6) of the auxiliary variable $x$ for $\rho .$ Thus the final expression for computing area sample sizes is

$$n_d^{g1}=\frac{N_{d}n-\left(N-N_{d}D-n\right)\left(1/R_{ax}^2-1\right)}{N+D\left(1/R_{ax}^2-1\right)}.(2.10)$$

It is easy to prove that ${\displaystyle {\sum }_{d=1}^D{n}_d^{g1}=n}.$ The computed sample sizes are rounded to the nearest integer. Sometimes compromises must be made. It can be concluded by the examination of (2.10) that the sample size increases when the size of area $N_d$ increases, but not proportionally. Under certain circumstances, such as low homogeneity coefficient, low overall sample size $n$ or small size of area, $N_d$ can lead to negative area sample size $n_d^{g1}.$ In this situation the negative value is changed to zero. A special case occurs if the total variation is only between areas causing value one to the measure of homogeneity (2.6), and (2.10) is reduced to proportional allocation.

2.3 Model-assisted MC-allocation

Molefe and Clark (2015) have used the following composite estimator for estimating the mean of the study variable $y$ for area $d:$

$${\tilde{y}}_d^C=\left(1-{\phi }_d\right){\overline{y}}_{dr}+{\phi }_d{\widehat{\beta }}^{\prime }{\overline{X}}_d.(2.11)$$

This estimator is a combination of two estimators: the synthetic estimator ${\widehat{\overline{Y}}}_{d\left(\text{syn}\right)}={\widehat{\beta }}^{\prime }{\overline{X}}_d,$ where $\widehat{\beta }$ is the estimated regression coefficient and ${\overline{X}}_d$ is the area population means of auxiliary variables $x,$ and a direct estimator ${\overline{y}}_{dr}={\overline{y}}_d+{\widehat{\beta }}^{\prime }\left({\overline{x}}_d-{\overline{X}}_d\right),$ where ${\overline{y}}_d$ and ${\overline{x}}_d$ are the area $d$ sample means of $y$ and $x.$ We use one auxiliary variable in our study. The coefficients ${\phi }_d$ are set with the intent to minimize the MSE of the estimator (2.11). The approximated design-based MSE of the estimator under certain conditions and assumptions is given by the expression

$${\text{MSE}}_p\left({\tilde{y}}_d^C;{\overline{Y}}_d\right)\approx {\left(1-{\phi }_d\right)}^2{v}_{d\left(\text{syn}\right)}+{\phi }_d^2{B}_d^2,(2.12)$$

where $v_{d\left(\text{syn}\right)}$ is the sampling variance of the synthetic estimator ${\widehat{\overline{Y}}}_{d\left(\text{syn}\right)}$ and $B_d={\beta }_U^{{}^{\prime }}{\overline{X}}_d-{\overline{Y}}_d$ is the bias when ${\widehat{\overline{Y}}}_{d\left(\text{syn}\right)}$ is used to estimate ${\overline{Y}}_d,$ with ${\beta }_U$ denoting the approximate design-based expectation of $\widehat{\beta }.$

The population contains $N$ units and $D$ strata defined by areas, and stratified sampling is used. A random sample SRSWOR (Simple Random Sampling without Replacement) of $n_d$ units is selected from stratum $d\left(d=1,\dots ,D\right)$ containing $N_d$ units. The relative size of area $d$ is $P_d=N_d/N.$

A two-level linear model $\xi$ conditional on the values of $x$ is assumed, with uncorrelated stratum random effects $u_d$ and random effects ${\epsilon }_i:$

$$\begin{array}{ll}{y}_i\hfill & ={\beta }^{\prime }{x}_i+u_d+{\epsilon }_i\hfill \\ E_{\xi }\left(u_d\right)\hfill & =E_{\xi }\left({\epsilon }_i\right)=0\hfill \\ V_{\xi }\left(u_d\right)\hfill & ={\sigma }_{ud}^2\hfill \\ V_{\xi }\left({\epsilon }_i\right)\hfill & ={\sigma }_{ed}^2\hfill \end{array}\},(2.13)$$

where $i$ refers to all units in stratum $d.$ This model implies that $V_{\xi }\left(y_i\right)={\sigma }_{ud}^2+{\sigma }_{ed}^2$ for all population units and ${\mathrm{cov}}_{\xi }\left(y_i,y_j\right)$ equals ${\rho }_d{\sigma }_d^2$ for units $i\ne j$ in the same stratum and zero for units from different strata, where ${\rho }_d={\sigma }_{ud}^2/\left({\sigma }_{ud}^2+{\sigma }_{ed}^2\right).$ A simplifying assumption that ${\rho }_d=\rho$ are equal for all strata is defined.

After making some other simplifying assumptions and solving the optimal weight ${\phi }_d$ in (2.12), the final approximate optimum anticipated MSE or approximate model assisted mean squared error is obtained of (2.12):

$${\text{AMSE}}_d=E_{\xi }{\text{MSE}}_p\left({\tilde{y}}_d^C\left[{\phi }_{d\left(\text{opt}\right)}\right];{\overline{Y}}_d\right)\approx {\sigma }_d^2\rho \left(1-\rho \right){\left[1+\left(n_d-1\right)\rho \right]}^{-1}.(2.14)$$

Next the criterion $F$ using anticipated MSE’s of the small area mean and overall mean estimators for model-assisted allocation is defined and developed into the final approximative form:

$$\begin{array}{ll}F\hfill & ={\displaystyle \sum _{d=1}^D{N}_d^q}{\text{AMSE}}_d+G{N}_{+}^{\left(q\right)}{E}_{\xi }{\mathrm{var}}_p\left({\widehat{\overline{Y}}}_r\right)\hfill \\ \hfill & \approx {\displaystyle \sum _{d=1}^D{N}_d^q{\sigma }_d^2\rho \left(1-\rho \right){\left[1+\left(n_d-1\right)\rho \right]}^{-1}}+G{N}_{+}^{\left(q\right)}{\displaystyle \sum _{d=1}^D{\sigma }_d^2{P}_d^2{n}_d^{-1}\left(1-\rho \right).}(2.15)\hfill \end{array}$$

Optimal sample sizes for the areas are obtained by minimizing (2.15) subject to ${\displaystyle {\sum }_d{n}_d=n}.$ Expression (2.15) follows the idea of Longford (2006). The weight $N_d^q$ reflects the inferential priority (importance) for area $d,$ with $0\le q\le 2,$ and $N_{+}^{\left(q\right)}={\displaystyle {\sum }_{d=1}^D{N}_d^q}.$ The quantity $G$ is a relative priority coefficient on the population level. Ignoring the goal of estimating the population mean corresponds to $G=0,$ and the attention is then only focused on area level estimation. On the other hand, the larger the value of $G,$ the more the second component in (2.15) dominates and the more the area level estimation is ignored.

We assume first that the population estimation has no priority $\left(G=0\right)$ and the unit survey cost are fixed. In this case minimization of (2.15) with respect of $n_d$ has a unique solution

$$n_{d\text{},\text{opt}}=\frac{n\sqrt{{\sigma }_d^2{N}_d^q}}{{\displaystyle {\sum }_{d=1}^D\sqrt{{\sigma }_d^2{N}_d^q}}}+\frac{1-\rho }{\rho }\left(\frac{\sqrt{{\sigma }_d^2{N}_d^q}}{D^{-1}{\displaystyle {\sum }_{d=1}^D\sqrt{{\sigma }_d^2{N}_d^q}}}-1\right).(2.16)$$

The formula (2.16) contains two unknown parameters, the intra-class correlation $\rho$ and the area-specific variance ${\sigma }_d^2.$ We replace $\rho$ with an adjusted homogeneity coefficient of the auxiliary variable $x.$ This coefficient is an approximation of the ICC (Intra-Class Correlation) (Section 2.2). Parameter ${\sigma }_d^2$ is replaced with the variance of $x$ in area $d.$ The reason for both replacements is that $y$ is correlated with $x.$ If also the population estimation has a priority $\left(G>0\right)$ then (2.16) does not apply and $F$ must be minimized numerically by using, for example, the NLP method, as we have done (Excel Solver, NLP option).

• Previous
• Article main page
• Next