Model-based small area estimation under informative sampling 4. Simulation studyModel-based small area estimation under informative sampling 4. Simulation study

• Previous
• Next

4.1 Implementation

A design-model (pm) approach was used for the simulation study by generating data for the $N={\displaystyle {\sum }_{i=1}^M{N}_i}$ population units according to a specified model, and then selecting a sample according to a specified design. The process of generating population data and then selecting a sample is repeated $R$ times. We next describe the steps to implement the process. The population data, $y_{ij},$ for $M=99$ areas and $N_i=100$ units within each area $i$ were generated from the simple nested error linear regression model

$$y_{ij}={\beta }_0+{\beta }_1{x}_{ij}+{\nu }_i+e_{ij};i=1,\dots ,99;j=1,\dots ,100,(4.1)$$

where ${\beta }_0=1,{\beta }_1=1,v_i\stackrel{\text{iid}}{\sim }N\left(0,{\sigma }_v^2=0.5\right)$ and independent of $e_{ij}\stackrel{\text{iid}}{\sim }N\left(0,{\sigma }_e^2=2\right).$ The population $x_{ij}-$ values were generated from a gamma distribution with mean 10 and variance 50, and held fixed over the simulation of population $y_{ij}-$ values from (4.1).

We considered different sample sizes within areas by fixing $n_i=5$ for the first 33 areas, $n_i=7$ for the next 33 areas and $n_i=9$ for the last 33 areas. This was done to study the effect of unequal sample sizes on the choice of the augmenting variable $g_{ij}=g\left(p_{j|i}\right).$ Samples of specified sizes, $n_i,$ were selected within the areas with probabilities proportional to specified sizes, $b_{ij},$ using the Rao-Sampford (Rao 1965 and Sampford 1967) method of sampling with unequal probabilities and without replacement. The latter method ensures that the inclusion probabilities ${\pi }_{j|i}$ are proportional to the sizes $b_{ij},\text{i}\text{.e}.,{\pi }_{j|i}=n_i{b}_{ij}/B_i=n_i{p}_{j|i},j=1,\dots ,N_i,$ where $B_i$ is the total of the $b_{ij}$ in area $i.$

We considered two different choices of the sizes $b_{ij}$ in the simulation study. The first choice uses

$$\begin{array}{lll}{b}_{ij}\hfill & =\hfill & \exp \left[\left\{-\left(y_{ij}-{\beta }_0-{\beta }_1{x}_{ij}\right)/{\sigma }_e+{\delta }_{ij}/5\right\}/3\right]\hfill \\ \hfill & =\hfill & \exp \left[\left\{-\left(v_i+e_{ij}\right)/{\sigma }_e+{\delta }_{ij}/5\right\}/3\right],\hfill \end{array}(4.2)$$

where ${\delta }_{ij}\stackrel{\text{iid}}{\sim }N\left(0,1\right).$ The size measures (4.2) are equivalent to those used by Pfeffermann and Sverchkov (2007) in their simulation study and satisfy the relationship (1.2) on the weights $w_{j|i}={\pi }_{j|i}^{-1}.$ Following PS, the area effects $v_i$ and the unit errors $e_{ij}$ were truncated to $\pm 2.5{\sigma }_v$ and $\pm 2.5{\sigma }_e$ to avoid extreme selection probabilities.

The second choice of size measures, following Asparouhov (2006), involves two different types of size measures: invariant (I) and non-invariant (NI). For the invariant case, $b_{ij}$ is independent of $v_i$ given $x_{ij};$ otherwise, it is called non-invariant. Invariant size measures are given by

$$b_{ij}={\left[1+\exp \left\{-\tau \left(\frac{1}{\alpha }{e}_{ij}+\sqrt{1-\frac{1}{{\alpha }^2}}{e}_{ij}^{*}\right)\right\}\right]}^{-1}.(4.3)$$

Non-invariant size measures are taken as

$$b_{ij}={\left[1+\exp \left\{-\tau \left(\frac{1}{\alpha }\left(v_i+e_{ij}\right)+\sqrt{1-\frac{1}{{\alpha }^2}}\left(v_i^{*}+e_{ij}^{*}\right)\right)\right\}\right]}^{-1}.(4.4)$$

The coefficient $\tau$ in (4.3) and (4.4), chosen as 0.5, ensures that the variation of the weights $w_{j|i}$ would not be too large within a simulation run. The random pair $\left(v_i^{*},e_{ij}^{*}\right)$ was generated independently of $\left(v_i,e_{ij}\right)$ from the same distributions as $v_i$ and $e_{ij}$ to ensure that the weight variation would be comparable between various levels of $\alpha .$ If some of the ${\pi }_{j|i}$ exceeded one, they were set to one, and the probabilities were recomputed for the remaining units. The $\alpha -$ values in (4.3) and (4.4), chosen as 1, 2, 3 or $\infty ,$ control the level of informativeness. Increasing $\alpha$ decreases informativeness, with $\alpha =\infty$ corresponding to non-informative sampling. Various dependencies in the simulations were introduced as follows, in order to increase the precision of comparisons between different estimators: All the four error components $\left(v_i,e_{ij},v_i^{*},e_{ij}^{*}\right)$ were first generated. Population $y-$ values, as well as invariant and non-invariant probabilities of selection, were then generated from those errors. For a given generated population, eight samples were selected: an invariant sample and a non-invariant sample for each value of $\alpha$ considered.

It may be noted that the weights $w_{j|i}$ obtained from the size measures (4.3) and (4.4) may not satisfy condition (1.2) of PS. We nevertheless fitted (1.2) to those weights to compute $\widehat{b}$ needed in the bias-adjusted estimator ${\widehat{\overline{Y}}}_i^{\text{PS}}.$

Using the design-model (pm) approach, $R=\text{1,000}$ samples were generated under the size measures (4.2) and the size measures (4.3) and (4.4). From each simulated sample $r\left(r=1,\dots ,R\right),$ the estimates ${\widehat{\overline{Y}}}_i^{H\left(r\right)},{\widehat{\overline{Y}}}_{i\left(a\right)}^{H\left(r\right)}$ and ${\widehat{\overline{Y}}}_i^{\text{PS}\left(r\right)}$ were computed for each small area $i;$ for the YR method only ${\widehat{\mu }}_i^{\text{YR}\left(r\right)}$ and ${\widehat{\mu }}_{i\left(a\right)}^{\text{YR}\left(r\right)}$ were computed. Also, the MSE estimates, $\text{mse}{\left({\widehat{\mu }}_i^H\right)}^{\left(r\right)},\text{mse}{\left({\widehat{\mu }}_{i\left(a\right)}^H\right)}^{\left(r\right)},\text{mse}{\left({\widehat{\mu }}_i^{\text{YR}}\right)}^{\left(r\right)}$ and $\text{mse}{\left({\widehat{\mu }}_{i\left(a\right)}^{\text{YR}}\right)}^{\left(r\right)},$ associated with ${\widehat{\mu }}_i^H,{\widehat{\mu }}_{i\left(a\right)}^H,{\widehat{\mu }}_i^{\text{YR}}$ and ${\widehat{\mu }}_{i\left(a\right)}^{\text{YR}},$ were computed. As noted earlier, we did not include the bootstrap MSE estimator of ${\widehat{\overline{Y}}}_i^{\text{PS}},$ proposed by Pfeffermann and Sverchkov (2007), in the simulation study. Also, for simplicity, we did not include the MSE estimators of ${\widehat{\overline{Y}}}_i^H$ and ${\widehat{\overline{Y}}}_{i\left(a\right)}^H$ because the latter estimators performed similarly to ${\widehat{\mu }}_i^H$ and ${\widehat{\mu }}_{i\left(a\right)}^H$ in terms of MSE.

We considered the following performance measures for a given estimator, say of the small area mean ${\overline{Y}}_i.$ Average absolute bias $\left(\overline{\text{AB}}\right)$ is measured by

$$\overline{\text{AB}}=\frac{1}{M}{\displaystyle \sum _{i=1}^M{\text{AB}}_i}$$

with

$${\text{AB}}_i=\left|\frac{1}{R}{\displaystyle \sum _{r=1}^R\left({\widehat{\overline{Y}}}_i^{\left(r\right)}-{\overline{Y}}_i^{\left(r\right)}\right)}\right|$$

where ${\widehat{\overline{Y}}}_i^{\left(r\right)}$ and ${\overline{Y}}_i^{\left(r\right)}$ are the values of ${\widehat{\overline{Y}}}_i$ and ${\overline{Y}}_i$ for the $r^{\text{th}}$ simulated sample and population. Efficiency of an estimator ${\widehat{\overline{Y}}}_i$ is measured by the average root MSE

$$\overline{\text{RMSE}}=\frac{1}{M}{\displaystyle \sum _{i=1}^M\sqrt{\frac{1}{R}{\displaystyle \sum _{r=1}^R{\left({\widehat{\overline{Y}}}_i^{\left(r\right)}-{\overline{Y}}_i^{\left(r\right)}\right)}^2}}}.$$

Turning to the performance of MSE estimators $\text{mse}\left({\widehat{\mu }}_i^H\right),\text{mse}\left({\widehat{\mu }}_{i\left(a\right)}^H\right),\text{mse}\left({\widehat{\mu }}_i^{\text{YR}}\right)$ and $\text{mse}\left({\widehat{\mu }}_{i\left(a\right)}^{\text{YR}}\right)$ in estimating MSEs, we first calculated reliable measures of MSEs by increasing $R=\text{1,000}$ to $T=\text{10,000}$ simulated samples. The MSE of an estimator ${\widehat{\mu }}_i$ is then calculated as

$$\text{MSE}\left({\widehat{\mu }}_i\right)=\frac{1}{T}{\displaystyle \sum _{t=1}^T{\left({\widehat{\mu }}_i^{\left(t\right)}-{\overline{Y}}_i^{\left(t\right)}\right)}^2},$$

where ${\widehat{\mu }}_i^{\left(t\right)}$ and ${\overline{Y}}_i^{\left(t\right)}$ denote the values of ${\widehat{\mu }}_i$ and ${\overline{Y}}_i$ for the $t^{\text{th}}$ simulated sample and population. For MSE estimation, we retained the original $R$ simulated samples and calculated the expected values $E\left[\text{mse}\left({\widehat{\mu }}_i\right)\right]=R^{-1}{\displaystyle {\sum }_{r=1}^R\text{mse}{\left({\widehat{\mu }}_i\right)}^{\left(r\right)}},$ where $\text{mse}{\left({\widehat{\mu }}_i\right)}^{\left(r\right)}$ denotes the value of the MSE estimate for the $r^{\text{th}}$ simulated sample. The average absolute relative bias $\left(\overline{\text{ARB}}\right)$ of a MSE estimator $\text{mse}\left({\widehat{\mu }}_i\right)$ is then calculated as

$$\overline{\text{ARB}}\left[\text{mse}\left({\widehat{\mu }}_i\right)\right]=M^{-1}{\displaystyle \sum _{i=1}^M\left|\frac{E\left[\text{mse}\left({\widehat{\mu }}_i\right)\right]}{\text{MSE}\left({\widehat{\mu }}_i\right)}-1\right|}.$$

4.2 Results under the Pfeffermann and Sverchkov size measures

Table 4.1 reports the simulation results on the average absolute bias $\left(\overline{\text{AB}}\right)$ and the average root mean square error $\left(\overline{\text{RMSE}}\right)$ of the estimators ${\widehat{\overline{Y}}}_i^H,{\widehat{\overline{Y}}}_{i\left(a\right)}^H,{\widehat{\mu }}_i^{\text{YR}},{\widehat{\mu }}_{i\left(a\right)}^{\text{YR}}$ and ${\widehat{\overline{Y}}}_i^{\text{PS}}$ under the PS size measures (4.2). The average absolute RB $\left(\overline{\text{ARB}}\right)$ of the MSE estimators $\text{mse}\left({\widehat{\mu }}_i^H\right),\text{mse}\left({\widehat{\mu }}_{i\left(a\right)}^H\right),\text{mse}\left({\widehat{\mu }}_i^{\text{YR}}\right)$ and $\text{mse}\left({\widehat{\mu }}_{i\left(a\right)}^{\text{YR}}\right)$ are also reported. Four different choices of the augmenting variable $g_{ij}$ were studied: $p_{j|i},w_{j|i},n_i{w}_{j|i}=p_{j|i}^{-1}$ and $\log p_{j|i}.$ Bootstrap estimator of $\text{MSE}\left({\widehat{\overline{Y}}}_i^{\text{PS}}\right),$ proposed by Pfeffermann and Sverchkov (2007), is not included in our study because the bootstrap simulation is very computer intensive.

Table 4.1 shows that the $\overline{\text{AB}}$ of the EBLUP estimator ${\widehat{\overline{Y}}}_i^H$ is large (= 0.456) relative to the corresponding augmented model EBLUP, ${\widehat{\overline{Y}}}_{i\left(a\right)}^H,$ for the four choices of $g_{ij}.$ Also, the choice $g_{ij}=w_{j|i}$ leads to larger $\overline{\text{AB}}$ compared to the other three choices (0.131 compared to 0.042 or less). The customary pseudo-EBLUP, ${\widehat{\mu }}_i^{\text{YR}},$ surprisingly performed well $\left(\overline{\text{AB}}=0.044\right)$ even though it was obtained under the assumption of noninformative sampling. This good performance is perhaps due to the use of weights in ${\widehat{\mu }}_i^{\text{YR}}.$ Augmented pseudo-EBLUP, ${\widehat{\mu }}_{i\left(a\right)}^{\text{YR}},$ leads to further reduction in $\overline{\text{AB}}.$ The PS estimator, ${\widehat{\overline{Y}}}_i^{\text{PS}},$ performs well relative to ${\widehat{\overline{Y}}}_{i\left(a\right)}^H:\overline{\text{AB}}=\mathrm{0.033.}$

Turning to $\overline{\text{RMSE}},$ Table 4.1 shows that ${\widehat{\overline{Y}}}_i^H$ has the largest value (= 0.617) due to large $\overline{\text{AB}},$ followed by ${\widehat{\mu }}_i^{\text{YR}}$ and ${\widehat{\overline{Y}}}_i^{\text{PS}}$ with values 0.442 and 0.416 respectively. On the other hand, the augmented model estimators performed significantly better relative to ${\widehat{\overline{Y}}}_i^{\text{PS}}$ and ${\widehat{\mu }}_i^{\text{YR}}.$ For example, the choice $g_{ij}=p_{j|i}$ gives $\overline{\text{RMSE}}=\mathrm{0.151.}$ Among the four choices of $g_{ij},$ the choice $w_{j|i}$ gives the largest $\overline{\text{RMSE}}\left(=0.242\right).$ We also calculated $\overline{\text{AB}}$ and $\overline{\text{RMSE}}$ of the approximate EBLUP estimators ${\widehat{\mu }}_i^H$ and ${\widehat{\mu }}_{i\left(a\right)}^H.$ We found that the values are practically the same as the corresponding values for ${\widehat{\overline{Y}}}_i^H$ and ${\widehat{\overline{Y}}}_{i\left(a\right)}^H.$

Finally, with respect to MSE estimation, $\text{mse}\left({\widehat{\mu }}_i^H\right)$ exhibits largest $\overline{\text{ARB}}:53.1\%$ compared to 3.8% for ${\widehat{\mu }}_i^{\text{YR}},$ although $\overline{\text{RMSE}}$ for ${\widehat{\mu }}_i^{\text{YR}}$ is larger compared to ${\widehat{\mu }}_{i\left(a\right)}^H$ based on $p_{j|i}$ or $n_i{w}_{j|i}.$ The MSE estimators $\text{mse}\left({\widehat{\mu }}_{i\left(a\right)}^H\right)$ and $\text{mse}\left({\widehat{\mu }}_{i\left(a\right)}^{\text{YR}}\right)$ lead to small $\overline{\text{ARB}}\left(<7\%\right)$ except for the choice $w_{j|i}$ which leads to $\overline{\text{ARB}}=62.6\%$ for ${\widehat{\mu }}_{i\left(a\right)}^H$ and $\overline{\text{ARB}}=39.6\%$ for ${\widehat{\mu }}_{i\left(a\right)}^{\text{YR}}.$

4.3 Selection of the augmenting variable

In this section we illustrate the selection of the augmenting variable by generating data for the $N$ population units from model (4.1) and then selecting a sample from the population data according to the Rao-Sampford method using size measures (4.2). Letting $u_{ij}=v_i+e_{ij},$ we fitted the model $y_{ij}={\beta }_0+{\beta }_1{x}_{ij}+u_{ij}$ to the sample data by ordinary least squares (OLS) and obtained the residuals ${\tilde{u}}_{ij}=y_{ij}-{\tilde{\beta }}_0-{\tilde{\beta }}_1{x}_{ij},$ where ${\tilde{\beta }}_0$ and ${\tilde{\beta }}_1$ are the OLS estimators of ${\beta }_0$ and ${\beta }_1$ respectively.

Figure 4.1 gives residual plots of $\left({\tilde{u}}_{ij},p_{j|i}\right),\left({\tilde{u}}_{ij},\log p_{j|i}\right),\left({\tilde{u}}_{ij},n_i{w}_{j|i}\right)$ and $\left({\tilde{u}}_{ij},w_{j|i}\right).$ All four plots clearly indicate informative sampling. Linear relationships between $u_{ij}$ and the two choices $p_{j|i}$ and $\log p_{j|i}$ suggest that either of them should work well. The choice $w_{j|i}$ indicates some non-linearity and wider scatter in the residual plot, and this choice led to the largest $\overline{\text{RMSE}}$ among the four choices, as shown in Table 4.1. The choice $n_i{w}_{j|i}$ also indicates some non-linearity but less scatter in the residual plot.

We have also fitted the augmented model (1.4) with $g\left(p_{j|i}\right)=p_{j|i}$ and calculated the OLS residuals ${\tilde{u}}_{0ij}=y_{ij}-{\tilde{\beta }}_{00}-{\tilde{\beta }}_{01}{x}_{ij}-{\tilde{\delta }}_0{p}_{j|i}.$ All the residuals ${\tilde{u}}_{0ij}$ are less than 2.0 in absolute value, suggesting adequacy of the augmented model.

4.4 Results under Asparouhov size measures

Table 4.2 reports the simulation results on $\overline{\text{AB}}$ under the Asparouhov size measures (4.3) and (4.4). It shows, as in Table 4.1 for the PS size measures, that $\overline{\text{AB}}$ of the EBLUP is large (0.437 for the invariant size measures (I) and 0.440 for non-invariant size measures (NI)) when the augmenting variable, $g_{ij},$ is not included in the model and sampling is very informative $\left(\alpha =1\right).$ Also, $\overline{\text{AB}}$ decreases as $\alpha$ increases. On the other hand, under the same model $\overline{\text{AB}}$ associated with pseudo-EBLUP is much lower: 0.048 for I and 0.047 for NI when $\alpha =1,$ and $\overline{\text{AB}}$ decreases as $\alpha$ increases. The PS estimator under the same model also exhibits lower $\overline{\text{AB}}$ (about 0.01) regardless of the choice of the value of $\alpha .$ Inclusion of $p_{j|i}$ or $n_i{w}_{j|i}$ or $\log p_{j|i}$ as augmenting variable in the model also leads to small $\overline{\text{AB}}$ for the EBLUP (0.02 or less) regardless of the value of $\alpha .$ On the other hand, the choice $w_{j|i}$ as the augmenting variable leads to larger $\overline{\text{AB}}$ (0.14 for $\alpha =1$ and 2), except for non-informative sampling $\left(\alpha =\infty \right).$ This poor performance of the choice $w_{j|i}$ is probably due to the fact that $w_{j|i}={\left(n_i{p}_{j|i}\right)}^{-1}$ depends on $n_i$ when the area sample sizes, $n_i,$ are not equal, unlike the other choices of $g_{ij}.$ Pseudo-EBLUP performed similarly to EBLUP under the augmented model in terms of $\overline{\text{AB}}.$

Table 4.3 reports the simulation results on the average root mean squared error $\left(\overline{\text{RMSE}}\right)$ using the Asparouhov size measures (4.3) and (4.4). It shows that the EBLUP, based on model (1.4) without the augmenting variable $g_{ij},$ has the largest $\overline{\text{RMSE}}$ (0.596 for I and 0.619 for NI) when the sampling is very informative $\left(\alpha =1\right).$ The $\overline{\text{RMSE}}$ gradually decreases to around 0.42 as the sampling becomes non-informative $\left(\alpha =\infty \right).$ On the other hand, $\overline{\text{RMSE}}$ of both the pseudo-EBLUP (without the $g_{ij}-$ term in the model) and PS do not depend on $\alpha ,$ and lead to significant reduction: $\overline{\text{RMSE}}$ of the pseudo-EBLUP is around 0.44 and $\overline{\text{RMSE}}$ of PS is slightly smaller, around 0.42. Increase in $\overline{\text{RMSE}}$ of the pseudo-EBLUP and PS over the EBLUP under non-informative sampling $\left(\alpha =\infty \right)$ is also small. On the other hand, EBLUP and pseudo-EBLUP under the augmented model lead to large reduction in MSE when the sampling is very informative $\left(\alpha =1\right),$ particularly for the choices $p_{j|i}$ and $\log p_{j|i}:\overline{\text{RMSE}}$ less than 0.15. The choice of $w_{j|i}$ leads to larger $\overline{\text{RMSE}}$ (around 0.29) when $\alpha =1$ but it is still much smaller than the $\overline{\text{RMSE}}$ for the pseudo-EBLUP without the $g_{ij}-$ term and the PS. As $\alpha$ increases, $\overline{\text{RMSE}}$ is roughly the same for EBLUP (under the augmented model), pseudo-EBLUP and PS.

Table 4.4 reports the simulation result on the average absolute relative bias $\left(\overline{\text{ARB}}\right)$ of MSE estimators under the Asparouhov size measures (4.3) and (4.4). It shows that $\overline{\text{ARB}}$ of the MSE estimator of the EBLUP, based on the model without the augmenting variable $g_{ij},$ is very large when the sampling is very informative $\left(\alpha =1\right):$ 52.8% for I and 59.1% for NI. $\overline{\text{ARB}}$ gradually decreases to around 5% under non-informative sampling $\left(\alpha =\infty \right).$ The use of $\log p_{j|i}$ as an augmenting variable leads to large reduction in $\overline{\text{ARB}}$ $\left(<9\%\right)$ and the choices $p_{j|i}$ and $n_i{w}_{j|i}$ also perform well in terms of $\overline{\text{ARB}}$ except for the case of NI and $\alpha =1$ which leads to 18.5% and 12.9% respectively. Again, $w_{j|i}$ is not a good choice because it leads to $\overline{\text{ARB}}$ as large as 40% when $\alpha =1.$ The MSE estimator associated with the pseudo-EBLUP (without $g_{ij})$ also performs well, except for NI and $\alpha =1,$ leading to $\overline{\text{ARB}}$ of 19.5%. Use of $\log p_{j|i}$ as auxiliary variable leads to $\overline{\text{ARB}}$ less than 8% for the MSE estimator associated with the pseudo-EBLUP. We have not included the PS bootstrap MSE estimator in our study.

Overall, our simulation study indicates that the use of augmented models under informative sampling leads to EBLUPs that perform well in terms of $\overline{\text{AB}}$ and $\overline{\text{RMSE}}$ of the estimators, and $\overline{\text{ARB}}$ of MSE estimators, provided that the augmenting variable is chosen properly. The bias-adjusted estimators of PS also perform well, even though they led to larger $\overline{\text{RMSE}}$ under the PS size measures (4.2). Pseudo-EBLUP estimators (without the augmenting variable) also perform well and further improvement may be achieved under augmented models.

• Previous
• Next