3. Estimation of area-specific MSPE

Jiming Jiang, Thuan Nguyen and J. Sunil Rao

Previous | Next

The design-based area-specific MSPE is defined as

$$\text{MSPE}\left({\widehat{\theta }}_i\right)=\text{E}{\left({\widehat{\theta }}_i-{\theta }_i\right)}^2,(3.1)$$

where and hereafter $\text{E}$ represents the design-based expectation, and ${\widehat{\theta }}_i$ is the OBP of ${\theta }_i,$ given by (1.5) with $\psi ={\left({\beta }^{\prime }\mathrm{,}\gamma \right)}^{\prime }$ replaced by its BPE, $\widehat{\psi }={\left({\widehat{\beta }}^{\prime }\mathrm{,}\widehat{\gamma }\right)}^{\prime }.$ As noted in Jiang et al. (2011), it is difficult to obtain second-order unbiased area-specific MSPE estimator under potential model misspecification. This is because the standard asymptotic techniques used in this area, such as the Prasad-Rao linearization method (Prasad and Rao 1990), and the jackknife method (Jiang, Lahiri and Wan 2002), are no longer applicable when the underlying model is misspecified. Jiang et al. (2011) used a different technique to derive a linearization MSPE estimator which is second-order unbiased. However, the latter is not guaranteed nonnegative. Furthermore, the leading term of this MSPE estimator is an $O(1)$ function of the area-specific data, rather than all the data. More precisely, the leading term for the estimated MSPE of ${\widehat{\theta }}_i,$ where ${\widehat{\theta }}_i$ is the OBP of the $i^{\text{th}}$ small area mean, ${\theta }_i,$ is ${\left({\widehat{\theta }}_i-y_i\right)}^2+D_i\left(2{\widehat{B}}_i-1\right),$ under the Fay-Herriot model, where $y_i$ is the observation from the $i^{\text{th}}$ area (the direct estimator), $D_i$ is the (known) sampling variance, ${\widehat{B}}_i=\widehat{A}/\left(\widehat{A}+D_i\right),$ and $\widehat{A}$ is the BPE for the variance of the area-specific random effect. This is the leading term because its order is $O(1),$ while the rest of the terms in the expression of the estimated MSPE are of the order $O\left(m^{-1}\right)$ or lower. Because $y_i$ is an observation from a single small area, it has a relatively large variance, that is, the variance is $O\left(1\right),$ if $n_i$ is bounded. On the other hand, the BPE $\widehat{A}$ is obtained using data from all the small areas, and therefore has a relatively (much) smaller variance; and ${\widehat{\theta }}_i$ is a mixture of $y_i$ and the BPEs. As a result, ${\left({\widehat{\theta }}_i-y_i\right)}^2$ is the main contributor to the variance, which can be quite large due to the variation of $y_i.$ On the other hand, the term $D_i\left(2{\widehat{B}}_i-1\right)$ can be negative. Thus, as a result of the high variation of ${\left({\widehat{\theta }}_i-y_i\right)}^2,$ there is a non-vanishing probability (as $m$ increases) that the leading term, hence the estimated MSPE, is negative. If we are to take a similar linearization approach under the NER model, we can derive a second-order unbiased MSPE estimator that involves ${\overline{y}}_{i\cdot }$ in its leading term, which is based on data from a single small area. Then, once again, we run into the problem of large variation and non-vanishing probability of negative value for the MSPE estimator.

Jiang et al. (2011) also used a parametric bootstrap method to obtain an alternative MSPE estimator; however, the justification for the use of this method is questionable given the potential model misspecification. Here we propose to use the nonparametric bootstrap following Efron's original idea (Efron 1979). The method does not rely on the NER model, hence is not affected by the model misspecification. Therefore, the current method is better justified. Furthermore, the proposed MSPE estimator is guaranteed nonnegative, and positive with probability one, which is a major advantage over the linearization MSPE estimator of Jiang et al. (2011).

Suppose that the small area subpopulations, or the $N_i’\text{s,}$ are large enough, so that the sampling from the subpopulations can be treated approximately as with replacement. Let $z_{ij}={\left({x^{\prime }}_{ij}\mathrm{,}{y}_{ij}\right)}^{\prime }\mathrm{,}j=1,\dots \mathrm{,}{n}_i$ denote the (original) samples from the $i^{\text{th}}$ small area, $1\le i\le m.$ We then draw samples, $z_{ij}^{\left(a\right)}={\left[{\left\{x_{ij}^{\left(a\right)}\right\}}^{\prime },y_{ij}^{\left(a\right)}\right]}^{\prime },j=\mathrm{1,}\dots \mathrm{,}{n}_i,$ with replacement, from $\left\{z_{ij}\mathrm{,}j=\mathrm{1,}\dots \mathrm{,}{n}_i\right\},$ independently for $1\le i\le m.$ Suppose that $B$ bootstrap samples are drawn, yielding samples $z^{\left(a\right)}=\left\{z_{ij}^{\left(a\right)}\mathrm{,1}\le j\le n_i\mathrm{,1}\le i\le m\right\},1\le a\le B.$ The bootstrapped version of the BP (1.5) is

$${\tilde{\theta }}_i^{\left(a\right)}={{\overline{X}}^{\prime }}_{i\cdot }\beta +\left\{r_i+\left(1-r_i\right)\frac{n_i\gamma }{1+n_i\gamma }\right\}\left[{\overline{y}}_{i\cdot }^{\left(a\right)}-{\left\{{\overline{x}}_{i\cdot }^{\left(a\right)}\right\}}^{\prime }\beta \right]\mathrm{,}(3.2)$$

where $\beta$ and $\gamma$ are the same population parameters for $\beta$ and $\gamma ,$ respectively, as for the original population. Note that the original samples of $z_{ij}$ are assumed to satisfy the same NER model, (1.4), with $X_{ik}\left(Y_{ik}\right)$ replaced by $x_{ij}\left(y_{ij}\right).$ Because the original samples are treated as the bootstrap population, following Efron's original idea, the population parameters, $\beta ,\gamma ,$ for the bootstrap samples are the same as those for the original samples. Nevertheless, as mentioned, the proposed bootstrap procedure is nonparametric in the sense that the assumed model, (1.4), plays no role in drawing the bootstrap samples. In particular, the BPE of $\beta$ and $\gamma ,$ based on the original samples, are not used anywhere in the bootstrapping; and the population quantities of interest are ${\overline{Y}}_i,1\le i\le m,$ whose bootstrap analogies are ${\overline{y}}_{i\cdot }\mathrm{,1}\le i\le m.$ This is different from the parametric bootstrap of Jiang et al. (2011), where the BPE of the model parameters, based on the original samples, are used to draw bootstrap samples under the assumed model. Also note that, because the ${\overline{X}}_i$  are known, they are treated as known constants, and therefore do not change during the bootstrapping (it does not make sense to "estimate� something that one already knows). Other than those, the procedure follows closely the standard bootstrap idea (e.g., Efron and Tibshirani (1993); also see Chatterjee, Lahiri and Li (2008) for an application to small area estimation). The bootstrap estimator of $\text{MSPE}\left({\widehat{\theta }}_i\right)=\text{E}{\left({\widehat{\theta }}_i-{\overline{Y}}_i\right)}^2$ is

$$\widehat{\text{MSPE}}\left({\widehat{\theta }}_i\right)=\frac{1}{B}{\displaystyle \sum _{a=1}^B}{\left\{{\widehat{\theta }}_i^{\left(a\right)}-{\overline{y}}_{i\cdot }\right\}}^2\mathrm{,}(3.3)$$

where ${\widehat{\theta }}_i^{\left(a\right)}$ is (3.2) with $\beta \mathrm{,}\gamma$ replaced by their BPE based on the bootstrapped samples.

Note. One might be concerned that, because the $n_i’\text{s}$ may be small in typical SAE problems, there may not be many distinct bootstrap samples for each small area. However, the data consist of not just one, but many small areas. When all of the small areas are combined, there are, still, a lot of distinct bootstrap samples, even if the $n_i’\text{s}$ are small.

We evaluate the performance of the proposed MSPE estimator by considering the simulated example of Subection 2.1 with $b=0.5$ but under smaller sample sizes. Namely, we start with the basic sample size $m=10$ and $n_i=5,$ and then either increase $n_i,$  from 5 to 10, or increase $m,$ from 10 to 20. We first consider the design-based bias of $\widehat{\text{MSPE}}\left({\widehat{\theta }}_i\right).$ Two finite populations are generated, and then fixed, so that the finite population for $m=10$ is a subpopulation of the finite population for $m=20.$  Table 3.1 reports, for the first ten small areas (these are all the small areas that are common under different values of $m),$ the simulated true MSPE (MSPE), obtained the same way as in Section 2, the simulated mean of $\widehat{\text{MSPE}}\left({\widehat{\theta }}_i\right)\left(\widehat{\text{MSPE}}\right),$ and the percentage relative bias (%RB) defined as

$$100\times \left\{\frac{\text{E}\left(\widehat{\text{MSPE}}\right)-\text{True MSPE}}{\text{True MSPE}}\right\}\mathrm{,}$$

where the expectation is based on the simulations. Another measure of performance is the square root of the mean squared error (RMSE) over the simulations, defined as

$$\sqrt{\frac{1}{K}{\displaystyle \sum _{k=1}^K}{\left({\widehat{\text{MSPE}}}_{i\mathrm{,}k}-{\text{MSPE}}_i\right)}^2}$$

for the $i^{\text{th}}$ small area, where ${\text{MSPE}}_i$ is the true MSPE for the $i^{\text{th}}$ small area (which does not depend on $k),$ evaluated over the simulations, and ${\widehat{\text{MSPE}}}_{i\mathrm{,}k}$ is the MSPE estimate based on the $k^{\text{th}}$ simulated data set. We consider $B=100$ as the number of bootstrap samples used to evaluate the MSPE estimator, (3.3). All results are based on 1,000 simulation runs. It is seen that, overall, the results improve when either $n_i$ or $m$ increase, but, in terms of %RB, the improvement is more universal, or effective, when $n_i$ increases. This is mainly due to the fact that, as $n_i$  increases, the sample provides a better approximation to the population; hence, the bootstrap distribution better approximates to the population distribution. Also note that, depending on the area, the sign of the RB can be either positive or negative. This is mainly due to the area-to-area difference (recall that the populations are fixed) as well as the bootstrap errors. To obtain some overall measures, we report the mean and standard deviation (s.d.) of the %RBs over the ten small areas as follows: $m=\mathrm{10,}{n}_i=5:$ mean $=4.2\%,$ s.d. $=14.8\%;m=\mathrm{10,}{n}_i=10:$ mean $=1.5\%,$ s.d. $=4.2\%;m=\mathrm{20,}{n}_i=5:$ mean $=-0.6\%,$ s.d. $=8.1\%.$ The boxplots of the %RBs are presented in Figure 3.1. The plots further illustrate the pattern of improvement. On the other hand, in terms of RMSE, the improvement is much more significant when $m$ increases than when $n_i$ increases. This is because having a larger $m$ reduces the MSPEs, in general; hence, natually, the correponding MSPE estimates also drop. In other words, both the estimator and the parameter (the MSPE) decrease, which typically results in a reduction in RMSE. The summary and boxplots for RMSE are omitted.

In addition, the %RB and RMSE in Table 3.1 fluctuates quite a bit from area to area. This is mainly due to the area to area difference. Recall the small area populations are generated each with population size $N_i=1\text{,}000,$ and then fixed throughout the simulation. Although the superpopulation used to generate the small area populations, including $X$ and $Y,$ are the same, there are still some differences in the generated finite populations, especially because the population size, $N_i,$ is not very large.

Table 3.1
Empirical Performance of $\widehat{\text{MSPE}}$
Table summary
This table displays the results of Empirical Performance of $\widehat{\text{MSPE}}$. The information is grouped by $m$ (appearing as row headers), $n_i$, MSPE, $\widehat{\text{MSPE}}$ , %RB and RMSE (appearing as column headers).

$m$	$n_i$	$i$	MSPE	$\widehat{\text{MSPE}}$	%RB	RMSE	$i$	MSPE	$\widehat{\text{MSPE}}$	%RB	RMSE
10	5	1	0.041	0.042	4.5	0.103	6	0.034	0.043	26.3	0.070
10	10	1	0.036	0.036	-0.4	0.068	6	0.034	0.036	6.4	0.070
20	5	1	0.031	0.032	4.1	0.051	6	0.028	0.031	12.5	0.046
10	5	2	0.046	0.038	-16.1	0.078	7	0.032	0.040	25.4	0.078
10	10	2	0.035	0.033	-4.1	0.078	7	0.033	0.034	2.7	0.068
20	5	2	0.031	0.029	-7.2	0.050	7	0.030	0.031	3.6	0.055
10	5	3	0.038	0.042	10.2	0.121	8	0.042	0.042	-0.4	0.150
10	10	3	0.037	0.036	-1.7	0.091	8	0.033	0.035	7.5	0.067
20	5	3	0.031	0.032	4.4	0.052	8	0.030	0.031	4.1	0.058
10	5	4	0.056	0.052	-7.6	0.121	9	0.050	0.042	-15.0	0.074
10	10	4	0.037	0.040	6.3	0.072	9	0.034	0.034	-1.0	0.063
20	5	4	0.040	0.035	-11.3	0.068	9	0.034	0.030	-11.1	0.049
10	5	5	0.033	0.037	11.8	0.066	10	0.041	0.043	3.1	0.082
10	10	5	0.032	0.033	2.5	0.066	10	0.034	0.033	-2.9	0.073
20	5	5	0.024	0.025	2.9	0.052	10	0.035	0.033	-7.9	0.062

Figure 3.1 Boxplots of %RB. $1:m=\mathrm{10,}{n}_i=5;2:m=\mathrm{10,}{n}_i=10;3:m=\mathrm{20,}{n}_i=5.$

Description for Figure 3.1

We conclude this section with some comments on the theoretical side. While there have been extensive studies on MSPE estimation in SAE since Prasad and Rao's seminal paper (Prasad and Rao 1990), the vast majority of these work focus on the model-based MSPEs. See, for example, Datta, Kubokawa, Molina and Rao (2011), Lahiri (2012), and Torabi and Rao (2012) for some recent work on design-based MSPE estimation in SAE. As noted in Jiang et al. (2011), under possible model misspecification, the model-based area-specific MSPE is not consistently estimable, and this is true for the design-based, area-specific MSPE as well. The reason is that, when the model is misspecified in terms of the mean function, the MSPE is not a function of a finite number of parameters (such as $\beta ,\gamma ,$ and ${\sigma }_e^2).$ In fact, because we operate under possible model misspecification, the quantities such as ${\overline{Y}}_i^2\mathrm{,1}\le i\le m$ are involved in the expressions of the area-specific MSPEs, which should all be treated as unknown parameters. Furthermore, the effective sample size for estimating ${\overline{Y}}_i^2$ is $n_i,$ if the assumed model fails. It follows that ${\overline{Y}}_i^2$ cannot be consistently estimated using data from the area alone, if $n_i$ is bounded. Generally speaking, if the MSPE can be estimated consistently, the difference between the MSPE estimator and the MSPE is $O_{\text{P}}\left(m^{-1/2}\right);$ therefore, the bias is typically $O\left(m^{-1}\right)$ without bias correction. On the other hand, if the (area-specific) MSPE cannot be consistently estimated, the difference between the MSPE estimator and the MSPE is typically $O_{\text{P}}\left\{{\left(m\wedge n_i\right)}^{-1/2}\right\},$ where $m\wedge n=\min \left(m\mathrm{,}n\right),$ hence the bias is typically $O\left\{{\left(m\wedge n_i\right)}^{-1}\right\},$ without the bias correction. The bootstrap MSPE estimator, $\widehat{\text{MSPE}},$ has the latter property, plus that it is always nonnegative. Although it is possible to bias-correct $\widehat{\text{MSPE}}$ to reduce the order of the bias to $o\left\{{\left(m\wedge n_i\right)}^{-1}\right\}$ (e.g., Hall and Maiti 2006), the nonnegative property may be lost after the bias correction. In view of the above discussion, it seems that, under the potential model misspecification, it is reasonable to define the first and second-order unbiasedness of an area-specific MSPE estimator in terms of $O\left\{{\left(m\wedge n_i\right)}^{-1}\right\}$ and $o\left\{{\left(m\wedge n_i\right)}^{-1}\right\},$ instead of the traditional $O\left(m^{-1}\right)$ and $o\left(m^{-1}\right)$ (e.g., Rao 2003).

Previous | Next