Small area estimation methods under cut-off sampling
Section 8. Simulation experiments

• Previous
• Next

8.1  Aims and general description

In this section, we describe simulation experiments designed to compare the small sample properties of the estimators of ${\overline{Y}}_i$ discussed above in the context of cut-off sampling. Specifically, we compare the naïve direct estimator ${\widehat{\overline{Y}}}_i^{\text{HA}},$ calibration estimators ${\widehat{\overline{Y}}}_i^{\text{LCAL}}$ and ${\widehat{\overline{Y}}}_i^{\text{LCALN}},$ and the EBLUP under the nested error model ${\widehat{\overline{Y}}}_i^{\text{EBLUP}},$ under two different scenarios. In the first scenario, the values of the target variable for all the population units are generated from the same model; in the second, included and excluded units are generated from different models.

In the absence of cut-off sampling, calibration estimators are design-consistent as the domain size $n_i$ increases even if the corresponding model does not hold, but this is not the case for model-based estimators. On the other hand, under the corresponding model, the EBLUP of a linear parameter is approximately the most efficient linear and unbiased estimator, so making simulations under a model would not provide any additional knowledge. The purpose here is to see whether the model-based predictors also perform well with respect to the (cut-off sampling) design. For this reason, we run design-based simulations by generating one population vector $y\mathrm{<mo>=</mo>}{\left(y_1^{\prime }\mathrm{,}\dots \mathrm{,}{y}_m^{\prime }\right)}^{\prime }$ from the nested error model in (5.1), keeping it fixed and repeatedly drawing a new cut-off sample in each MC simulation. Allocation of units to the sets of included or excluded units is done by generating a random binary variable $c_{ij}$ for each unit $j\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}{N}_i$ and area $i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}m.$ The units $j$ with $c_{ij}\mathrm{<mo>=</mo>}1$ are assigned to $U_{iI}$ and those with $c_{ij}\mathrm{<mo>=</mo>}0$ to $U_{iE}.$ In each Monte Carlo (MC) replicate, samples are drawn, independently for each domain $i,$ from the $U_{iI}$ units, $i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}m.$

8.2  Common regression model

We consider a population of $N\mathrm{<mo>=</mo>}$ 20,000 individuals divided into $m\mathrm{<mo>=</mo>}$ 80 domains with the same size $N_i\mathrm{<mo>=</mo>}$ 250, $i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}m.$ We consider three auxiliary variables, with values generated as $x_{ij\kappa }\stackrel{\text{iid}}{{\displaystyle \sim }}N\left(3,2\right),$ $\kappa \mathrm{<mo>=</mo>}\mathrm{1,}\mathrm{2,}3.$ The binary variables $c_{ij}$ determining the allocation of units in $U_{iI}$ or $U_{iE}$ for each domain $i$ are generated independently as $c_{ij}\stackrel{\text{ind}}{{\displaystyle \sim }}\text{Bern}\left(p_{j|i}\right),$ where the probabilities $p_{j|i}\mathrm{<mo>=</mo>}\text{Pr}\left(c_{ij}\mathrm{<mo>=</mo>\; 1}\right)$ are related to the vector of auxiliary variables $x_{ij}\mathrm{<mo>=</mo>}{\left(x_{ij1}\mathrm{,}{x}_{ij2}\mathrm{,}{x}_{ij3}\right)}^{\prime }$ in the form

$$p_{j|i}\mathrm{<mo>=</mo>}\frac{\exp \left(x_{ij}^{\prime }\zeta \right)}{1+\exp \left(x_{ij}^{\prime }\zeta \right)}\mathrm{,}j\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}{N}_i\mathrm{,}i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}m\mathrm{.}$$

We take $\zeta \mathrm{<mo>=</mo>}{\left(\text{0}\text{.75}\mathrm{,}\mathrm{1,}1\right)}^{\prime }.$ Based on this value, the total number of included units (with $c_{ij}\mathrm{<mo>=</mo>}1)$ from all the domains represents roughly half of the population.

The values of the target variable $y_{ij}$ are generated from the nested error model (5.1) using $x_{ij}\mathrm{<mo>=</mo>}{\left(x_{ij1}\mathrm{,}{x}_{ij2}\mathrm{,}{x}_{ij3}\right)}^{\prime }$ and taking $\beta \mathrm{<mo>=</mo>}{\left(\mathrm{1,}\text{1}\text{.5}\mathrm{,}1\right)}^{\prime },$ ${\sigma }_u^2\mathrm{<mo>=</mo>}{\left(\text{0}\text{.75}\right)}^2$ and ${\sigma }_e^2\mathrm{<mo>=</mo>}{4}^2,$ which leads to a determination coefficient $R^2\approx \text{0}\text{.5}\text{.}$ Then, keeping the population values $\left\{\left(x_{ij}\mathrm{,}{y}_{ij}\mathrm{,}{c}_{ij}\right)\mathrm{<mo>;</mo>}j\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}{N}_i\mathrm{,}i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}m\right\}$ fixed, we draw $K\mathrm{<mo>=</mo>}$ 1,000 Monte Carlo samples $s^{\left(\text{}k\text{}\right)},$ $k\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}K.$ Each of these samples is obtained by drawing independent domain sub-samples $s_i^{\left(\text{}k\text{}\right)}$ of size $n_i$ from the units in $U_{iI}$ by SRSWOR, $i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}m.$ The domain sample sizes are taken as $n_i\in \left\{\mathrm{5,}\mathrm{10,}\mathrm{30,}50\right\},$ with each sample size repeated for 20 subsequent domains. With the data from the $k^{\text{th}}$ sample, we compute the basic direct estimator, calibration estimators at the domain level (LCAL) and at the population level (LCALN), and EBLUP. Weights, $h_{j|i}$ and $g_{j|i},$ in the calibration estimators (4.3) and (4.6) respectively are obtained using the function calib from package sampling (Tillé and Matei, 2016) of R (R Development Core Team, 2016). EBLUPs are obtained using R package sae (Molina and Marhuenda, 2015), which by default estimates the model parameters ${\sigma }_u^2,$ ${\sigma }_e^2$ and $\beta$ using restricted maximum likelihood (REML).

Let ${\widehat{\overline{Y}}}_i$ be a generic estimator of ${\overline{Y}}_i$ and ${\widehat{\overline{Y}}}_i^{\left(\text{}k\text{}\right)}$ its value obtained with $k^{\text{th}}$ sample. We evaluate the performance of estimators in terms of relative bias (RB) and relative root MSE (RRMSE) under the design, approximated empirically as

$${\text{RB}}_{\pi }\left({\widehat{\overline{Y}}}_i\right)\mathrm{<mo>=</mo>}100\frac{K^{-1}{\displaystyle {\sum }_{k\mathrm{<mo>=</mo>\; 1}}^K\left({\widehat{\overline{Y}}}_i^{\left(\text{}k\text{}\right)}-{\overline{Y}}_i\right)}}{{\overline{Y}}_i}\mathrm{,}{\text{RRMSE}}_{\pi }\left({\widehat{\overline{Y}}}_i\right)\mathrm{<mo>=</mo>}100\frac{\sqrt{K^{-1}{{\displaystyle {\sum }_{k\mathrm{<mo>=</mo>\; 1}}^K\left({\widehat{\overline{Y}}}_i^{\left(\text{}k\text{}\right)}-{\overline{Y}}_i\right)}}^2}}{{\overline{Y}}_i}\mathrm{.}$$

Averages across domains of absolute RB and of RRMSE are also calculated as

$$\overline{\text{ARB}}\mathrm{<mo>=</mo>}{m}^{-1}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^m}\left|{\text{RB}}_{\pi }\left({\widehat{\overline{Y}}}_i\right)\right|\mathrm{,}\overline{\text{RRMSE}}\mathrm{<mo>=</mo>}{m}^{-1}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^m}{\text{RRMSE}}_{\pi }\left({\widehat{\overline{Y}}}_i\right)\mathrm{.}$$

Figure 8.1 displays boxplots of percent RB for the considered estimators of the mean ${\overline{Y}}_i,$ where each boxplot is for the 20 domains in each group of sample sizes $n_i\mathrm{<mo>=</mo>}\mathrm{5,}\mathrm{10,}\mathrm{30,}50.$ We can see the large cut-off sampling bias of the basic direct estimator, with median RB exceeding 20% for all the domain sample sizes. This cut-off sampling bias is corrected by all the other estimators. Nevertheless, the LCALN estimator shows wider boxplots. This estimator gets large bias for some domains probably because its assisting model is not accounting for the domain effects. The LCAL estimator is based on a model that accounts for domain effects and performs well in terms of design bias uniformly for all the domain sample sizes, although EBLUP also performs rather well in terms of design bias.

Looking now at the RRMSE in Figure 8.2, we can see the much smaller RRMSEs of EBLUPs for all the domain sample sizes. The LCAL estimator gets closer RRMSEs as the domain sample size grows, but for $n_i\mathrm{<mo>=</mo>\; 5}$ it gets huge RRMSEs. We have seen that the LCALN can be substantially biased for some domains and it also has large RRMSEs for all the domain sample sizes. Thus, in summary, EBLUP exhibits the lowest design RRMSE and at the same time keeps the design bias under control.

Table 8.1 reports averages across all the domains of absolute RB and RRMSE, together with % share of squared bias from the total design MSE. We can see again the large cut-off sampling bias of the basic direct estimator, with a bias share of $B_{\pi }^2/{\text{MSE}}_{\pi }\approx$ 100%, in contrast to all other estimators. The LCAL estimator has the smallest average ARB, followed closely by EBLUP. LCALN performs the best in terms of bias ratio because of its large MSE. Thus, we consider that LCAL performs better. As already said, EBLUP clearly performs the best when looking at both bias and MSE.

8.3  Different regression models

In this simulation experiment, we preserve the same population values and sampling scheme as before, but the values of the target variable for the included and excluded units are generated from models with different parameter values. Of course, this is not a favorable scenario for the considered model-based estimators, but it may be realistic since, in practice, the assumed model cannot be checked for the excluded units. Thus, instead of a constant $\beta$ for all the population units, we take ${\beta }_I\mathrm{<mo>=</mo>}{\left(\mathrm{1,}\mathrm{1.5,}1\right)}^{\prime }$ for the included units and ${\beta }_E\mathrm{<mo>=</mo>}{\left(\text{0}\text{.5}\mathrm{,}\text{1}\text{.6}\mathrm{,}\text{0}\text{.5}\right)}^{\prime }$ for the excluded ones. The values of the explanatory variables and variance components ${\sigma }_u^2$ and ${\sigma }_e^2$ are taken exactly as before. Again, we draw $K\mathrm{<mo>=</mo>}$ 1,000 samples $s^{\left(\text{}k\text{}\right)}$ by independent SRSWOR within the units in domain $i$ with $c_{ij}\mathrm{<mo>=</mo>}1,$ with the same domain sample sizes $n_i$ as before. With the sample data from the $k^{\text{th}}$ sample, we compute basic direct, LCAL, LCALN and EBLUP estimates of ${\overline{Y}}_i.$

Figure 8.3 shows boxplots of the corresponding percent RBs for each domain sample size. In this case, all the estimators are biased, but the bias of the basic direct estimator becomes huge, exceeding 40% for some of the domains. The bias of LCAL and EBLUP is kept relatively small for all the domains, but that of LCALN estimator is still very large in absolute value for some of the domains. In absence of cut-off sampling, the calibration estimators are asymptotically design-unbiased as the domain sample size $n_i$ increases, even if the considered model does not hold. However, this is not true under cut-off sampling and for this reason the RBs of calibration estimators do not decrease as $n_i$ grows. Even under this unfavorable scenario of different generating models for included and excluded units, EBLUP shows a moderate bias, which is comparable to that of LCAL estimator, and performs clearly the best in terms of RRMSE.

Again, averages across all the domains of absolute RB and RRMSE are shown in Table 8.2, together with sq. bias ratio. As already noted, the basic direct estimator has a huge bias, whereas LCAL and EBLUP estimators keep an $\overline{\text{ARB}}$ below 10%. LCALN displays the lowest bias ratio because of a larger MSE. Again, EBLUP shows the best performance in terms of efficiency, with an average RRMSE also below 10%.

The simulation experiment was repeated taking a value of ${\beta }_E$ further away from ${\beta }_I,$ making the two regression models differ substantially. Results are not included due to space constraints but, as one would expect, RB and RRMSE values increase for all estimators, but conclusions are similar to the last experiment. The basic direct estimator gets the largest RB, calibration estimators and EBLUP clearly reduce the cut-off sampling bias of the basic direct estimator and EBLUP gets smaller RRMSE, specially for the domains with the smallest sample sizes.

• Previous
• Article main page
• Next