Local polynomial estimation for a small area mean under informative sampling
Section 5. Simulation study

• Previous
• Next

The set-up of the simulation study follows the one used in Verret et al. (2015). We considered a population with $M=$ 15 small areas and $N_i=$ 15 units within each small area. The relatively small number of small areas and units within areas were chosen so as to alleviate the computational burden. We used a single auxiliary variable $x.$ The population $x$ -values were generated from a gamma distribution with mean 10 and variance 50. The population $y_{ij}$ -values were generated by the following model

$$y_{ij}=4+x_{ij}+v_i+e_{ij};i=1,\dots ,15;j=1,\dots ,15,(5.1)$$

where $v_i\stackrel{\text{iid}}{{\displaystyle \sim }}N\left(0,{\sigma }_v^2\right)$ and $e_{ij}\stackrel{\text{iid}}{{\displaystyle \sim }}N\left(0,{\sigma }_e^2\right)$ with ${\sigma }_v^2=$ 0.5 and ${\sigma }_e^2=$ 2.

We considered a single sample size, $n_i=$ 3, within a small area. We used Conditional Poisson Sampling (CPS) to select unequal probability samples within the small areas, with probabilities proportional to specified sizes $c_{ij}$ (see Tillé, 2006, Chapter 5). We considered two different choices of the sizes $c_{ij}$ in the simulation study. The first choice uses

$$c_{ij}=\exp \left[\frac{1}{3}\left\{-\frac{(v_i+e_{ij})}{{\sigma }_e}+\frac{{\delta }_{ij}}{5}\right\}\right],(5.2)$$

where ${\delta }_{ij}\stackrel{\text{iid}}{{\displaystyle \sim }}N\left(0,1\right).$ The size measures (5.2) are equivalent to those used by Pfeffermann and Sverchkov (2007) in their simulation study and satisfy the relationship (2.5) on the weights $w_{j|i}={\pi }_{j|i}^{-1}.$

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, $c_{ij}$ is independent of $v_i$ given $x_{ij};$ otherwise, it is called non-invariant. Invariant size measures are given by

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

Non-invariant size measures are taken as

$$c_{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},(5.4)$$

where the random pair $\left(v_i^{*},e_{ij}^{*}\right)$ is generated independently of $\left(v_i,e_{ij}\right)$ using the same distributions as $v_i$ and $e_{ij}.$ These size measures were used by Asparouhov (2006). The coefficient $\tau$ controls for the variation of the weights and the value $\alpha$ controls the level of informativeness of the sampling design. We chose $\tau =$ 0.5 and $\alpha =1,2,3$ and $\infty$ corresponding to several levels of informativeness generated by $c_{ij}$ in (5.3) and (5.4). Increasing $\alpha$ decreases informativeness, with $\alpha =\infty$ corresponding to non-informative sampling. If some of the ${\pi }_{j|i}’\text{s}$ exceeded one, they were set to one, and the probabilities were recomputed for the remaining units.

5.1  Performance of the local polynomial estimator of ${\overline{Y}}_i$

We compared the bias and mean squared error of the estimators ${\widehat{\overline{Y}}}_i^{\text{EBLUP}},$ ${\widehat{\overline{Y}}}_i^{\text{VRH}}$ and ${\widehat{\overline{Y}}}_i^{\text{LP}}.$ The EBLUP estimator ${\widehat{\overline{Y}}}_i^{\text{EBLUP}}$ based on (1.1) assumes that the sample model coincides with the population model, thereby ignoring the informativeness of the sampling design. We studied two versions of ${\widehat{\overline{Y}}}_i^{\text{VRH}}$ investigated by Verret et al. (2015) for various choices of $g(\cdot )$ that account for informativeness. They are EBLUP estimators based on the augmented sample model (1.2). They are denoted as ${\widehat{\overline{Y}}}_i^{\text{VRH1}}$ when $g\left(p_{j|i}\right)=p_{j|i}$ and ${\widehat{\overline{Y}}}_i^{\text{VRH2}}$ when $g\left(p_{j|i}\right)=\log \left(p_{j|i}\right).$ We report results only for these $g$ functions, as they outperform others given in Verret et al. (2015). Finally, ${\widehat{\overline{Y}}}_i^{\text{LP}}$ represents our new local polynomial estimator.

The bias and the mean squared error of the estimators were computed using $R=$ 1,000 simulated samples selected under a design-model approach. For each run, $r=1,\dots ,R,$ we first generated the population $y_{ij}$ -values under the population model (5.1) and computed ${\overline{Y}}_i^{\left(r\right)},$ the mean of the small area $i$ in the $r^{\text{th}}$ generated population. Samples of sizes $n_i=3$ were then selected within the small areas using CPS with probabilities proportional to specified sizes $c_{ij}^{\left(r\right)}$ given by (5.2) for the Pfeffermann and Sverchkov (2007) (PS) size measures, and (5.3) and (5.4) corresponding to the invariant and non-invariant cases in the case of the Asparouhov (2006) (AP) size measures. From each simulated sample $r$   $\left(r=1,\dots ,R\right),$ the estimates ${\widehat{\overline{Y}}}_i^{\text{EBLUP}\left(r\right)},$ ${\widehat{\overline{Y}}}_i^{\text{VRH1}\left(r\right)},$ ${\widehat{\overline{Y}}}_i^{\text{VRH2}\left(r\right)}$ and ${\widehat{\overline{Y}}}_i^{\text{LP}\left(r\right)}$ were computed for each small area $U_i.$ An optimal bandwidth $h_{\text{opt}}^{\left(r\right)}$ was found for ${\widehat{\overline{Y}}}_i^{\text{LP}\left(r\right)}$ using the cross-validation criterion. A grid of the form (0.01, 0.02, 0.03,…, 0.15) covered the possible values for $h_{\text{opt}}^{\left(r\right)}$ in populations generated by (5.1).

For a given estimator of the small area mean ${\overline{Y}}_i,$ we considered the following performance measures:

Average Absolute Bias

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

where

$${\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|.$$

Average Root Mean Squared Error

$$\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)}^{\text{2}}}}}.$$

Table 5.1 reports on the average absolute bias $\left(\overline{\mathrm{AB}}\right)$ of estimators ${\widehat{\overline{Y}}}_i^{\text{EBLUP}},$ ${\widehat{\overline{Y}}}_i^{\text{VRH1}},$ ${\widehat{\overline{Y}}}_i^{\text{VRH2}}$ and ${\widehat{\overline{Y}}}_i^{\text{LP}}$ under the PS size measures (5.2) and AP size measures (5.3 and 5.4) for $\alpha =1,2,3$ and $\infty .$

As observed in Verret et al. (2015), the $\overline{\mathrm{AB}}$ of the EBLUP estimator ${\widehat{\overline{Y}}}_i^{\text{EBLUP}}$ with just the auxiliary variable $x,$ is quite a bit larger than those based on the augmented models $(p_{j|i}$ and $\log (p_{j|i})),$ and the local polynomial method. This holds regardless of how the size measures have been generated (PS or AP). The $\overline{\mathrm{AB}}$ of ${\widehat{\overline{Y}}}_i^{\text{EBLUP}}$ attains its highest value (0.431) when the design is very informative $\left(\alpha =1\right),$ and decreases as $\alpha$ increases. This observation also holds for the estimators based on the augmented models. The inclusion of $p_{j|i}$ or $\log \left(p_{j|i}\right),$ as an augmenting variable, in the model results in small $\overline{\mathrm{AB}}’\text{s,}$ with the highest being 0.036. Comparing the $\overline{\mathrm{AB}}’\text{s}$ of the local polynomial estimator ${\widehat{\overline{Y}}}_i^{\text{LP}}$ to those associated with the VRH augmented models, we observe that they are comparable for $\alpha =1$ and $\alpha =2,$ and slightly larger for $\alpha \ge 3.$

Table 5.2 reports the simulation results on the average root mean squared error $\left(\overline{\text{RMSE}}\right)$ of the estimators for both the PS size measures (5.2) and the AP size measures (5.3 and 5.4) for $\alpha =1,2,3$ and $\infty .$ The EBLUP, ${\widehat{\overline{Y}}}_i^{\text{EBLUP}},$ based on model (1.1) without the augmenting variable $g\left(p_{j|i}\right),$ has the largest $\overline{\text{RMSE}}’\text{s}$ (0.740 for I and 0.752 for NI) for the AP size measures corresponding to $\alpha =1,$ and 0.685 for the PS size measure. The $\overline{\text{RMSE}}$ decreases as $\alpha$ increases: 0.608 for I and 0.610 for NI in the case of non-informative sampling $\left(\alpha =\infty \right).$ The $\overline{\text{RMSE}}’\text{s}$ for ${\widehat{\overline{Y}}}_i^{\text{VRH1}},$ ${\widehat{\overline{Y}}}_i^{\text{VRH2}}$ and ${\widehat{\overline{Y}}}_i^{\text{LP}}$ are significantly smaller than those associated with ${\widehat{\overline{Y}}}_i^{\text{EBLUP}}$ when sampling is very informative $\left(\alpha =1\right)$ and for the PS size measure. There are small differences in terms of $\overline{\text{RMSE}}$ between our non-parametric approach and the parametric approach in Verret et al. (2015).

When the sampling is less informative $\left(\alpha =3\right),$ the local linear estimator ${\widehat{\overline{Y}}}_i^{\text{LP}}$ is better than ${\widehat{\overline{Y}}}_i^{\text{EBLUP}},$ but its $\overline{\text{RMSE}}$ is slightly larger than those associated with the parametric estimators ${\widehat{\overline{Y}}}_i^{\text{VRH1}}$ and ${\widehat{\overline{Y}}}_i^{\text{VRH2}}.$ In this case, we observe that the estimated function $m_0\left(p_{j|i}\right)$ is close to a flat line, and this implies that the local linear approximation is not as appropriate. This explains why ${\widehat{\overline{Y}}}_i^{\text{LP}}$ is slightly worse than ${\widehat{\overline{Y}}}_i^{\text{VRH1}}$ and ${\widehat{\overline{Y}}}_i^{\text{VRH2}}$ when the level of informativeness of the sampling is low. A local polynomial estimator performs well when the function $m_0(\cdot )$ is meaningfully non-constant.

When the sample is non-informative $\left(\alpha =\infty \right),$ ${\widehat{\overline{Y}}}_i^{\text{EBLUP}}$ is better than ${\widehat{\overline{Y}}}_i^{\text{VRH1}},$ ${\widehat{\overline{Y}}}_i^{\text{VRH2}}$ and ${\widehat{\overline{Y}}}_i^{\text{LP}}$ in both invariant and non-invariant case. This conclusion is somewhat different from that of Verret et al. (2015) where for $\alpha =\infty$ their estimators ${\widehat{\overline{Y}}}_i^{\text{EBLUP}},$ ${\widehat{\overline{Y}}}_i^{\text{VRH1}}$ and ${\widehat{\overline{Y}}}_i^{\text{VRH2}}$ have equal $\overline{\mathrm{AB}}$ and $\overline{\text{RMSE}}$ values. Verret et al. (2015) used both larger populations and samples, and this may explain why their augmented models produced estimators as good as the population model under non-informative sampling designs. Under our simulation set-up, we found that the $\overline{\mathrm{AB}}$ and $\overline{\text{RMSE}}$ of the EBLUP are small for $\alpha$ values larger than 6: this corresponds to a sample design that is almost non-informative. In this case, we recommend using EBLUP.

5.2  Performance of the MSE estimators

We now turn to the performance of the bootstrap procedures for estimating the MSEs of the EBLUP, VRH and local polynomial estimators. Let ${\widehat{\overline{Y}}}_i$ be an estimator of ${\overline{Y}}_i$ and ${\text{mse}}_{\text{boot}}\left({\widehat{\overline{Y}}}_i\right)$ be the bootstrap estimator of $\text{MSE}\left({\widehat{\overline{Y}}}_i\right).$ From $R=$ 1,000 simulated populations and samples, we first computed measures of MSE values as

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

where ${\overline{Y}}_i^{\left(r\right)}$ is the true mean, and ${\widehat{\overline{Y}}}_i^{\left(r\right)}$ is the value of the estimator for the $r^{\text{th}}$ population. Let ${\text{mse}}_{\text{boot}}\left({\widehat{\overline{Y}}}_i\right)$ be the bootstrap estimator of $\text{MSE}\left({\widehat{\overline{Y}}}_i\right).$ It is denoted as ${\text{mse}}_{\text{boot}}\left({\widehat{\overline{Y}}}_i^{\text{EBLUP}}\right)$ for the EBLUP estimator ${\widehat{\overline{Y}}}_i^{\text{EBLUP}},$ and corresponds to the parametric (unconditional) bootstrap method given by equation (4.2). For our local polynomial estimator ${\widehat{\overline{Y}}}_i^{\text{LP}}$ and the Verret et al. (2015) estimators, ${\widehat{\overline{Y}}}_i^{\text{VRH1}}$ and ${\widehat{\overline{Y}}}_i^{\text{VRH2}},$ the mse values, denoted as ${\text{mse}}_{\text{boot}}\left({\widehat{\overline{Y}}}_i^{\text{LP}}\right)$ and ${\text{mse}}_{\text{boot}}\left({\widehat{\overline{Y}}}_i^{\text{VRH }j}\right),$ for $j=1$ and $j=2$ respectively, are computed using the conditional parametric bootstrap method of Section 4. For each selected sample in the $r^{\text{th}}$  simulated population $\left(r=1,\dots ,R\right),$ we used $B=$ 400 bootstraps to compute the $r^{\text{th}}$ value of ${\text{mse}}_{\text{boot}}\left({\widehat{\overline{Y}}}_i\right),$ that we denote as ${\text{mse}}_{\text{boot}}^{\left(r\right)}\left({\widehat{\overline{Y}}}_i\right).$ We considered two measures to evaluate the performance of ${\text{mse}}_{\text{boot}}\left({\widehat{\overline{Y}}}_i\right)\text{:}$ average absolute relative bias and average confidence interval. These measures are defined as follows:

Average Absolute Relative Bias:

$$\overline{\text{ARB}}\text{=}\frac{1}{M}{\displaystyle \sum _{i=1}^M\left|\frac{E\left({\text{mse}}_{\text{boot}}\left({\widehat{\overline{Y}}}_i\right)\right)}{\text{MSE}\left({\widehat{\overline{Y}}}_i\right)}-1\right|},$$

where

$$E\left({\text{mse}}_{\text{boot}}\left({\widehat{\overline{Y}}}_i\right)\right)\text{=}\frac{1}{R}{\displaystyle \sum _{r=1}^R{\text{mse}}_{\text{boot}}^{\left(r\right)}\left(\widehat{\overline{Y}}\right).}$$

Average Confidence Level:

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

where ${\text{CL}}_i\text{=}{R}^{-1}{\displaystyle {\sum }_{r=1}^R\text{I}\left({\overline{Y}}_i^{\left(r\right)}\in {\text{IC}}^{\left(r\right)}\right)}$ and ${\text{IC}}^{\left(r\right)}=\left[{\widehat{\overline{Y}}}_i^{\left(r\right)}\pm 1.96\sqrt{{\text{mse}}_{\text{boot}}^{\left(r\right)}\left({\widehat{\overline{Y}}}_i\right)}\right].$

Table 5.3 reports simulation results on the average relative bias $\left(\overline{\text{ARB}}\right)$ of the MSE estimators for both the PS size measures (5.2) and Asparouhov size measures (5.3 and 5.4) for $\alpha =1,2,3$ and $\infty .$

The $\overline{\mathrm{ARB}}$ of ${\widehat{\overline{Y}}}_i^{\text{EBLUP}},$ based on the model without the augmenting variable $g\left(p_{j|i}\right),$ is very large when the sampling is very informative $\left(\alpha =1\right)\text{:}$ 39.9% for I and 46.6% for NI. The $\overline{\text{ARB}}$ gradually decreases to around 5% under non-informative sampling $\left(\alpha =\infty \right).$ The $\overline{\mathrm{ARB}}’\text{s}$ of both the parametric and non-parametric estimators are smaller in general than 10%, with the exception of 14.4% for the ${\widehat{\overline{Y}}}_i^{\text{VRH2}}$ estimator that uses $\log \left(p_{j|i}\right)$ as an augmenting variable.

Table 5.4 reports simulation results on the average confidence level $\left(\overline{\text{CL}}\right)$ associated with the MSE estimators for both the PS size measures (5.2) and the AP size measures (5.3 and 5.4) for $\alpha =1,2,3$ and $\infty$ and nominal level of 0.95.

The EBLUP estimator ${\widehat{\overline{Y}}}_i^{\text{EBLUP}}$ has the worst coverage when the sample design is very informative. The coverage improves as the design becomes less informative. The coverage of the other estimators is between 93% and 95%, with the exception of ${\widehat{\overline{Y}}}_i^{\text{VRH2}}$ (the one that includes $\log (p_{j|i}))$ with coverage slightly lower.

5.3  Inclusion of an augmenting variable

The local polynomial approach results in an automatic way of obtaining a reasonable augmented model that is a function of the selection probabilities $p_{j|i}.$ However, given that one does not know whether the design is informative or not, should we always include an augmenting variable in the model? If the sample design is not informative it is reasonable to use model (1.1). Note that in this case, including the augmenting variables, $p_{j|i}$ or $\log \left(p_{j|i}\right),$ has a very small impact either on the absolute relative bias of the estimator and absolute relative bias of the estimated MSE. A similar conclusion was obtained in Verret et al. (2015) who used a larger population and sample size.

The same question arises with respect to the use of the local polynomial procedure. In this case, the conclusions are not quite as clear. If the design is very informative, the local polynomial approach gains in terms of absolute bias and mean squared error when $\alpha =1$ or $\alpha =2.$ When the sampling design is less informative $\left(\alpha =3\right)$ the parametric approach in Verret et al. (2015) is the better choice, but by a very small margin.

In a practical situation, the value of $\alpha$ is not known and the decision to use the augmenting variable in a parametric or nonparametric model should be taken. To this end, we follow the suggested procedure in Verret et al. (2015) to provide some guidelines on how to decide on this choice for an arbitrary data set. Define $u_{ij}=v_i+e_{ij},$ and fit the following model $y_{ij}={\beta }_0+{\beta }_1{x}_{ij}+u_{ij}$ to the sample data by ordinary least squares (OLS). The residuals are ${\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 5.1 displays residual plots of $\left({\widehat{m}}_0\left(p_{j|i}\right),{\tilde{u}}_{ij}\right),j=1,\dots ,N_i;i=1,\dots ,M$ for the AP measures $\alpha =1,2,3$ and $\infty$ in the invariant case. For $\alpha =1,$ the relationship between ${\tilde{u}}_{ij}$ and ${\widehat{m}}_0\left(p_{j|i}\right)$ is clearly linear, suggesting that the design is informative. As $\alpha$ increases, the design is less informative. Note that $m_0\left(p_{j|i}\right)$ is constant when $\alpha =\infty .$ Similar observations hold for the non-invariant case. For the PS size measures the graph resembled the one given in Figure 5.1 when $\alpha =1.$

Table 5.5 provides the estimated correlation coefficients, $\widehat{\rho }\text{=}\text{cor}\left({\tilde{u}}_{ij},{\widehat{m}}_0\left(p_{j|i}\right)\right),$ for PS and AP size measures for $\alpha =1,2,3$ and $\infty .$

In terms of $\overline{\text{RMSE}},$ we noticed in Section 5.1 that ${\widehat{\overline{Y}}}_i^{\text{EBLUP}}$ is better than the estimators based on augmented models for $\alpha \ge 6.$ Results not presented in Table 5.5 show that for $\alpha \ge 6,$ the absolute value of the correlation coefficient is less than 0.1. On the basis of this limited simulation, a user could decide on the choice of the estimator to use for a real data set as follows: i. If $\left|\widehat{\rho }\right|$ is larger than 0.5, use ${\widehat{\overline{Y}}}_i^{\text{LP}}\text{};$ ii. If $\left|\widehat{\rho }\right|$ is less than 0.1, use ${\widehat{\overline{Y}}}_i^{\text{EBLUP}}\text{};$ iii. otherwise use ${\widehat{\overline{Y}}}_i^{\text{VRH1}}$ or ${\widehat{\overline{Y}}}_i^{\text{VRH2}}.$

• Previous
• Article main page
• Next