Comparison of some positive variance estimators for the Fay-Herriot small area model 6. Simulation results and analysisComparison of some positive variance estimators for the Fay-Herriot small area model 6. Simulation results and analysis

• Previous
• Next

6.1 Monte Carlo Distribution of the variance estimators

Table 6.1 shows that the REML variance estimator has the lowest bias $\left({\sigma }_v^2=1\right)$ and the highest variance. The lower efficiency of REML may be due to it not being a smooth function of the data caused by its split definition (3.1). The MIX estimator inherits some of this low efficiency. The other variance estimators have lower variability, higher positive bias but the conditional expectation of AM.YL and AR.YL given ${\widehat{\sigma }}_{\text{vREML}}^{\text{2}}=0$ is close to zero. The unconditional bias of AM.LL is higher than the unconditional bias of the MIX. By definition of the MIX estimator, the conditional bias of the MIX and AM.LL estimators coincide. The MIX estimator also converges faster than the other estimators. For example, given the probability distribution over the 10,000 variance estimates with $m=45,$ we calculated the probability of estimates lying within an interval containing ${\sigma }_v^2=1.$ The probability that the MIX estimates lie between 0.6 and 1.4 is 0.47 whereas the probability that AM.YL estimates lie between 0.6 and 1.4 is 0.16. Furthermore, the probability that MIX estimates are smaller than 0.2 is 0.05 whereas the probability that AM.YL estimates are smaller than 0.2 is 0.53.

6.2 True MSE of the EBLUP, average relative bias and average root relative MSE of the MSE estimators

All variance estimators are consistent and asymptotically normal with variance converging at the same rate. They differ in their bias: REML, AR.YL and MIX have bias of the order of $o\left(1/m\right)$ whereas AM.LL and AM.YL have bias of the order $O\left(1/m\right).$ The bias inherent in the last three methods impacts the estimation of the MSE of the EBLUP even for a moderate number of areas.

For $m=100,$ Tables 6.2a and 6.2b show that as ${\sigma }_v^2/{\psi }_i$ increases, the MSE of the EBLUP decreases and this relationship holds irrespective of the number of areas. We observe that the MSE of ${\widehat{\theta }}_i$ under the REML and the MIX variance estimators are slightly higher than the rest of the MSEs, due to the higher variability inherent in these variance estimators. Table 6.2a presents results for the Taylor linearization MSE estimator and the two parametric MSE estimators under REML, AM.LL, AR.YL and AM.YL variance estimation. Table 6.2b presents results for the following MSE estimators under the MIX variance estimation: RB_Y1 defined in (4.3), RB_Y2 defined in (4.2), M_et_al, defined in (4.5), PB MSE and naïve PB MSE estimators. Among the Taylor MSE estimators, RB_Y1 and M_et_al under MIX exhibit the lowest bias. Among the bootstrap MSE estimators PB under MIX and Naive PB under AR.YL exhibit the lowest bias. Turning to the RRMSE of the MSE estimators, it decreases as ${\sigma }_v^2/{\psi }_i$ increases. Differences between the RB_Y2 MSE estimator under the MIX and the Taylor MSE estimator under the AM.YL seem small but consistent. While ARB is lower for the RB_Y1, the M_et_al and the Naive MSE estimators under the MIX method than for the RB_Y2 under the MIX, and also lower for the Taylor and the Naive PB under the AR.YL method than for the RB_Y2 under the MIX, the opposite happens in terms of RRMSE. This can be explained in part due to the extreme negative conditional bias exhibited by these MSE estimators (i.e., the RB_Y1 and the M_et_al under the MIX and the Taylor and the Naive PB under the AR.YL method) as shown in Table 6.3. Even for $m=100$ there is a relatively high proportion (16%) of populations that yield ${\widehat{\sigma }}_{v\text{REML}}^2=0$ and in these populations, estimates from most variance methods and most MSE estimators are farthest below the true value. That is, for these MSE estimators, the conditional MSE estimators do not fare well. The PB MSE estimator seems to adjust well for bias, but it is more variable than the Naive PB MSE. When we also include the ARB, the RRMSE and the $AR{B}_{ℂ}$ in the evaluation, the RB_Y2 under the MIX method, followed closely by Naïve PB under the MIX seems to perform the best. This may suggest the superiority of RB_Y2 and Naive under MIX for $m=100,$ which is a moderate number of areas for this data.

Tables 6.4a and b below display results for $m=45$ with 9 areas per ${\sigma }_v^2/{\psi }_i.$ The AM.YL yields MSEs smaller than the MIX, with differences in MSEs of at most 2%. As the number of areas decreases, the bias of the variance estimators increase and the MSE estimators are affected by this. Indeed, the ARB of all MSE estimators have increased. In particular, the ARB of the Taylor MSE estimators under YL and LL variance estimation and the ARB of RB_Y2, have increased by 100% over the ARB with 100 areas. In terms of RRMSE, the Taylor MSE under the AM.YL has slightly lower RRMSE than the RB_Y2 under the MIX method for very small ${\sigma }_v^2/{\psi }_i.$ In general, the variability (in RRMSE) of the RB_Y2 is lower than that of the Taylor under LL and YL estimation and than that of the RB_Y1 and the M_et_al. This may be due in part to the underestimation of the MSEs for the populations with zero REML estimates, which, for $m=45,$ range around 30% of all populations. Table 6.5 illustrates this better: given ${\widehat{\sigma }}_{vR\text{EML}}^2=0,$ there is serious underestimation in RB_Y1 and M_et_al.

Taking into account the ARB, the RRMSE and the $AR{B}_{ℂ}$ of the MSE estimators, the Naive PB MSE estimator under the MIX performs the best for larger ${\sigma }_v^2/{\psi }_i.$ Table 6.6 displays performance measures, averaged over the five sampling variance groups, for the three Taylor MSE estimators under the MIX with data from the same model described in 5.1 but with three different values of ${\sigma }_v^2.$ The RB_Y2 performs better when ${\sigma }_v^2=1,$ but as ${\sigma }_v^2$ becomes smaller, the M_et_al MSE estimator has an advantage, precisely because it was constructed under the premise that ${\sigma }_v^2$ is approximately zero.

Tables 6.7a and 6.7b below show the outcomes for $m=15$ areas with 3 areas per ${\sigma }_v^2/{\psi }_i.$ Differences in MSEs per variance estimation method are at most 5%.

There is no monotone relationship between ARB or RRMSE and ${\sigma }_v^2/{\psi }_i,$ which could be an indication that the second order approximation to estimating the MSE is poor under every method of variance estimation. The ARB of all Taylor MSE estimators under the LL and the YL methods of variance estimation are unacceptably high and the same is true for the RRMSE. The RB_Y2 under the MIX does not fare well either. The reason for this last outcome is clear: the high % of zero REML estimates (43%) implies the MIX coincides with AM.LL for the zero REML populations. Thus, the MIX has a positive bias for $m=15,$ and the RB_Y2 does not account for this bias. The RB_Y1 accounts for the bias in the MIX, but the bias estimator is not very precise for $m=15.$ The M_et_al MSE estimator almost coincides with the ARB and RRMSE of the Taylor MSE estimator under the REML variance estimation, because by definition they are equal when ${\widehat{\sigma }}_{v\text{REML}}^2=0.$ The ${\text{ARB}}_{ℂ}$ of the three Taylor MSE estimators under the MIX is poor. Taking into account all performance measures, the bootstrap MSE estimators perform better than the Taylor MSE estimators. For $m=15$ areas with 3 areas per ${\sigma }_v^2/{\psi }_i,$ PB under MIX performs the best, followed by the Naive under AR.YL and AM.YL.

Summarizing, under the Fay-Herriot model with positive ${\sigma }_v^2,$ and among the positive variance estimators under study, the MIX and the AR.YL variance estimators are the only ones with negligible asymptotic bias. The AM.YL and the LL variance estimators have a larger asymptotic bias. On the other hand, our simulation showed that for a moderate number of areas and for populations that yield zero REML estimates, both YL variance estimators were negatively biased, and produced EBLUPs that were close to the synthetic estimator of the mean. In contrast, the MIX, built as the combination of the AM.LL and the REML, was only mildly negatively biased in these populations. Moreover, the unconditional distribution of the MIX approached normality much faster than those of the other variance estimators.

In terms of MSE of the EBLUP, there were considerable gains in precision over the direct estimator, under all methods of variance estimation considered here, even for a small number of areas. The AM.LL and both the AM.YL and the AR.YL variance estimators carried lower variability than the REML and the MIX. It impacted only minimally the MSE of the EBLUP: differences among MSEs for the same signal to noise ratio were small. These differences widened as either the number of areas or the signal to noise ratio decreased. Thus, it may possible that for an extremely low signal to noise ratio, the MSE under MIX would be somewhat larger than under the AM.YL variance estimator.

Under the MIX method of variance estimation, we compared three different Taylor-type MSE estimators and two bootstrap MSE estimators. All three Taylor estimators of the MSE under MIX (RB_Y1, RB_Y2 and M_et_al) are unbiased up to the second order. Also the Taylor-type estimators of the MSE under the LL and the YL are unbiased up to the second order. RB_Y1, AM.LL and AM.YL may yield negative MSE estimates.

The Taylor MSE under the REML method of variance estimation and the M_et_al under the MIX coincide by definition, hence their performance measures have negligible differences (their true MSEs are different, however in our study, for $m=100,$ the MIX coincided with the REML 84% of the time). For a moderate number of areas, which for this data could be $m=45\text{ or }100,$ and for populations that yield zero REML estimates, both the Taylor MSE estimators under the REML and the M_et_al MSE estimators do not account for the variation due to the estimation of ${\sigma }_v^2$ and this is reflected in their very negative ${\text{ARB}}_{ℂ},$ which is below -60% for the smaller signal to noise ratios. On the other hand, the RB_Y1 does account for the variation due to the estimation of ${\sigma }_v^2,$ but its ${\text{ARB}}_{ℂ}$ is also very negative: the RB_Y1 is a split MSE estimator that for populations with ${\widehat{\sigma }}_{v\text{REML}}^2=0,$ it subtracts a factor of the unconditional bias of the AM.LL, which is always positive, whereas a better formula for a split MSE estimator would be to use an estimator of the conditional bias $E\left({\widehat{\sigma }}_v^2/{\widehat{\sigma }}_{v\text{REML}}^2=0\right).$ Indeed, even for a moderate number of areas $\left(m=100\right),$ Table 6.1 shows that the unconditional bias of the MIX is 49% whereas the conditional bias of the MIX is -37%.

The PB MSE estimator under the AR.YL and the MIX methods adjusted well for the bias but paid in terms of variance. Among all the MSE estimators it appears that the Naive Bootstrap MSE estimator performed best, and even better under the MIX variance estimation, when taking into account the three measures ARB, ${\text{ARB}}_{ℂ}$ and RRMSE together. We found that for a moderate number of areas, the RB_Y2 had the lowest RRMSE among the Taylor estimators under the MIX method. On the other hand, M_et_al is most reliable when the true underlying variance ${\sigma }_v^2$ is very small: in this case M_et_al is effectively the MSE estimator of the synthetic estimator of the small area mean. We do not recommend relying on the second order approximation to the MSE when $m$ is small: the approximation (2.6) to the MSE does not necessarily hold, the performance measures obtained from our study are very unstable and they may vary from data set to data set.

In conclusion, under the hypothesis of ${\sigma }_v^2>0,$ the relative performances of competing positive variance estimators depend on the size of ${\sigma }_v^2,$ the signal to noise ratio, the number of areas and the objective function. For a moderate number of areas, the MIX variance estimator appeared to perform better than the LL and the YL estimators in this study; under the MIX method, the Naive PB MSE estimator had the lowest ${\text{ARB}}_{ℂ}$ and RRMSE combined; the M_et_al MSE estimator under the MIX variance estimator performed marginally better than the RB_Y1 when the underlying ${\sigma }_v^2$ was very small. However, the percentage of REML zeros yielded under the simulation model shows that an outcome of ${\widehat{\sigma }}_{v\text{REML}}^2=0$ and/or negative tests of hypothesis do not necessarily mean that ${\sigma }_v^2$ is sufficiently small to rely on M_et_al. In the absence of other information, the Naive PB estimator under the MIX appears to perform better.

Acknowledgements

The authors would like to thank Professor J.N.K. Rao from Carleton University for his useful comments and to Victor Estevao from Statistics Canada for developing the grid maximization especially for this project. We also would like to thank the reviewers for their careful appraisal of our paper and for their suggestions to improve this paper.

Appendix A

Proof of Theorem 4.1

The asymptotic variance of ${\widehat{\sigma }}_{v\text{MIX}}^2$ is given by: $\overline{V}\left({\widehat{\sigma }}_{v\text{MIX}}^2\right)=\underset{m\to \infty }{\lim }E{\left({\widehat{\sigma }}_{v\text{MIX}}^2-{\sigma }_v^2\right)}^2$

We show that $E{\left({\widehat{\sigma }}_{v\text{MIX}}^2-{\sigma }_v^2\right)}^2\le E{\left({\widehat{\sigma }}_{v\text{REML}}^2-{\sigma }_v^2\right)}^2+o\left(1/m\right)\text{ as }m\to \infty .$

$$\begin{array}{ll}E{\left({\widehat{\sigma }}_{v\text{MIX}}^2-{\sigma }_v^2\right)}^2\hfill & ={\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2>0\right\}}{\int }{\left({\widehat{\sigma }}_{v\text{REML}}^2-{\sigma }_v^2\right)}^{2}dP}+{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}{\int }{\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\sigma }_v^2\right)}^{2}dP}\hfill \\ \hfill & \begin{array}{l}\le {\displaystyle \underset{\Omega }{\int }{\left({\widehat{\sigma }}_{v\text{REML}}^2-{\sigma }_v^2\right)}^{2}dP}+{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}{\int }{\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\sigma }_v^2\right)}^{2}dP}=E{\left({\widehat{\sigma }}_{v\text{REML}}^2-{\sigma }_v^2\right)}^2\text{(A}\text{.1)}\\ +o\left(\frac{1}{m}\right).\end{array}\hfill \end{array}$$

Indeed, by the Holder and Minkowski inequalities, with any $1<p<\infty ,1/p+1/q=1,$ and setting $X\equiv {\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\sigma }_v^2\right)}^2=O_p\left(1/m\right)$ and the indicator $I\left({\widehat{\sigma }}_{\text{vREML}}^{\text{2}}=0\right)$ of populations with ${\widehat{\sigma }}_{v\text{REML}}^2=0,$ we have:

$$\begin{array}{ll}{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2<0\right\}}{\int }{\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\sigma }_v^2\right)}^{2}dP}\hfill & \le {\left({\displaystyle \underset{\Omega }{\int }{\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\sigma }_v^2\right)}^{2p}dP}\right)}^{1/p}.{\left(P\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}\right)}^{1/q}\hfill \\ \hfill & ={\left(O\left(\frac{1}{m^p}\right)\right)}^{1/p}.{\left(o\left(1\right)\right)}^{1/q}=o\left(\frac{1}{m}\right),\hfill \end{array}(\text{A}.2)$$

since ${\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\sigma }_v^2\right)}^2$ is uniformly bounded and ${\widehat{\sigma }}_{v\text{REML}}^2\stackrel{P}{{\displaystyle \to }}{\sigma }_v^2>0.$ Note that the AM.LL and REML estimators of ${\sigma }_v^2$ are uniformly bounded as a consequence of their almost sure convergence to ${\sigma }_v^2$ (see, for example, Yuan and Jennrich 1998).

Proof of Theorem 4.2

We denote by ${\widehat{\sigma }}_{v\text{ML}}^2$ the maximum likelihood variance estimator.

We show first that ${\widehat{\sigma }}_{v\text{REML}}^2-{\widehat{\sigma }}_{v\text{ML}}^2=O_p\left(1/m\right).$ Let $G_{*}\left({\sigma }_v^2\right)=\partial \log \left(L_{*}\right)/\partial {\sigma }_v^2=0$ be the estimating equation that yields the variance estimator *. Equation (3.4) implies:

$$G_{\text{AM}\text{.LL}}\left({\sigma }_v^2\right)-G_{\text{ML}}\left({\sigma }_v^2\right)=\partial \log {\sigma }_v^2/\partial {\sigma }_v^2=\frac{1}{m{\sigma }_v^2}=O\left(\frac{1}{m}\right).(\text{A}.3)$$

With $G_{\text{ML}}^{\prime }(\cdot )\triangleq \left(\partial G_{\text{ML}}/\partial {\sigma }_v^2\right)(\cdot )$ and $G_{\text{ML}}^{\prime \text{}\prime }(\cdot )\triangleq \left(\partial G_{\text{ML}}^{\prime }/\partial {\sigma }_v^2\right)(\cdot ),$ equation (A.3) implies:

$$G_{\text{ML}}^{\prime }\left({\sigma }_v^2\right)-G_{\text{AM}\text{.LL}}^{\prime }\left({\sigma }_v^2\right)=O\left(\frac{1}{m}\right).(\text{A}.4)$$

Now, using equation (A.4), the $\sqrt{m}-$ consistency of the ML and AM.LL estimators of ${\sigma }_v^2,$ the two-term Taylor expansion of $G_{\text{ML}}(\cdot )\text{ and }{G}_{\text{AM}\text{.LL}}(\cdot )$ at ${\sigma }_v^2$ and $G_{\text{ML}}^{\prime }\left({\sigma }_v^2\right)=O\left(1\right)$ as $m\to \infty ,$ the left- hand side in (A.3) is equal to:

$$\begin{array}{l}=G_{\text{ML}}^{\prime }\left({\sigma }_v^2\right)\left({\widehat{\sigma }}_{v\text{ML}}^2-{\sigma }_v^2\right)-G_{\text{AM}\text{.LL}}^{\prime }\left({\sigma }_v^2\right)\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\sigma }_v^2\right)+O_p\left(\frac{1}{m}\right)\\ =G_{\text{ML}}^{\prime }\left({\sigma }_v^2\right)\left({\widehat{\sigma }}_{v\text{ML}}^2-{\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2\right)+\left(G_{\text{ML}}^{\prime }\left({\sigma }_v^2\right)-G_{\text{AM}\text{.LL}}^{\prime }\left({\sigma }_v^2\right)\right)\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\sigma }_v^2\right)+O_p\left(\frac{1}{m}\right)\\ =G_{\text{ML}}^{\prime }\left({\sigma }_v^2\right)\left({\widehat{\sigma }}_{v\text{ML}}^2-{\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2\right)+O_p\left(\frac{1}{m^{3/2}}\right)+O_p\left(\frac{1}{m}\right).\end{array}$$

The last equality above implies

$${\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\widehat{\sigma }}_{v\text{ML}}^2=O_p\left(\frac{1}{m}\right)\text{ as }m\to \infty .(\text{A}.5)$$

Similarly, we establish a relationship between $G_{\text{REML}}\left({\sigma }_v^2\right)$ and $G_{\text{ML}}\left({\sigma }_v^2\right):$ given that $\text{tr}\left(V^{-1}Z{\left(Z^{\prime }{V}^{-1}Z\right)}^{-1}{Z}^{\prime }{V}^{-1}\right)=O\left(1\right)$ follows from conditions 1 through 3 in Section 3 and equation (3.1), we have:

$$G_{\text{REML}}\left({\sigma }_v^2\right)-G_{\text{ML}}\left({\sigma }_v^2\right)=\frac{1}{m}\text{tr}\left(V^{-1}Z{\left(Z^{\prime }{V}^{-1}Z\right)}^{-1}{Z}^{\prime }{V}^{-1}\right)=O\left(\frac{1}{m}\right)\text{ as }m\to \infty ,(\text{A}.6)$$

Equation (A.6) and the same argument as with the AM.LL estimator, imply:

$${\widehat{\sigma }}_{v\text{REML}}^2-{\widehat{\sigma }}_{v\text{ML}}^2=O_p\left(\frac{1}{m}\right)\text{ as }m\to \infty .(\text{A}.7)$$

Equations (A.5) and (A.7) combined, yield:

$$\left({\widehat{\sigma }}_{v\text{REML}}^2-{\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2\right)=O_P\left(\frac{1}{m}\right).(\text{A}.8)$$

Now we express the bias of the MIX estimator by:

$$B_{\text{MIX}}\left({\widehat{\sigma }}_{v\text{MIX}}^2\right)={\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2>0\right\}}{\int }\left({\widehat{\sigma }}_{v\text{REML}}^2-{\sigma }_v^2\right)dP}+{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}{\int }\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\sigma }_v^2\right)dP}.$$

We add and subtract ${\displaystyle {\int }_{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}\left({\widehat{\sigma }}_{v\text{REML}}^2-{\sigma }_v^2\right)dP}$ from the right-hand side of the equation above to obtain:

$$\begin{array}{ll}{B}_{\text{MIX}}\left({\widehat{\sigma }}_{v\text{MIX}}^2\right)\hfill & ={\displaystyle \underset{\Omega }{\int }\left({\widehat{\sigma }}_{v\text{REML}}^2-{\sigma }_v^2\right)dP}+{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}{\int }\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\widehat{\sigma }}_{v\text{REML}}^2\right)dP}\hfill \\ \hfill & =\text{Bias}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)+{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}{\int }\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\widehat{\sigma }}_{v\text{REML}}^2\right)dP}.\hfill \end{array}(\text{A}.9)$$

Now, since ${\widehat{\sigma }}_{v\text{AM}\text{.LL}}^{\text{2}}-{\widehat{\sigma }}_{v\text{REML}}^2$ is uniformly bounded, we apply the Holder and Minkowski inequality with $p=q=2$ and equation (A.8) to the last term in (A.9) to obtain:

$$\begin{array}{ll}{B}_{\text{MIX}}\left({\widehat{\sigma }}_{v\text{MIX}}^2\right)\hfill & =\text{Bias}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)+{\left({\displaystyle \underset{\Omega }{\int }{\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\widehat{\sigma }}_{v\text{REML}}^2\right)}^{2}dP}\right)}^{1/2}\cdot P{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}^{1/2}\hfill \\ \hfill & =\text{Bias}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)+O\left(\frac{1}{m}\right)\cdot o\left(1\right)=\text{Bias}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)+o\left(\frac{1}{m}\right).\hfill \end{array}(\text{A}.10)$$

Proof of Remark 4.2: ${\text{mse}}_0$ is unbiased up to the second order

$$\begin{array}{ll}E\left({\text{mse}}_0\right)-\text{MSE}\left({\widehat{\theta }}_i\right)\hfill & ={\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2>0\right\}}{\int }\left(g_{1i}+g_{2i}+2g_{3i}\right)\left({\widehat{\sigma }}_{v\text{REML}}^2\right)dP}+{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}{\int }{g}_{2i}}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)dP-\text{MSE}\hfill \\ \hfill & \begin{array}{l}=\left[{\displaystyle \underset{\Omega }{\int }\left(g_{1i}+g_{2i}+2g_{3i}\right)\left({\widehat{\sigma }}_{v\text{REML}}^2\right)dP}-\text{MSE}\right]\\ +{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}{\int }{g}_{2i}}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)dP-{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}{\int }\left(g_{1i}+g_{2i}+2g_{3i}\right)\left({\widehat{\sigma }}_{v\text{REML}}^2\right)dP}\end{array}\hfill \\ \hfill & =\left[o\left(\frac{1}{m}\right)\right]-{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}{\int }2g_{3i}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)dP},\hfill \end{array}(\text{A}.11)$$

since $g_{1i}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)=g_{1i}\left(0\right)=0$ in $\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}$ and $g_{2i}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)$ cancels out in (A.11). But

$$g_{3i}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)=g_{3i}\left(0\right)=\frac{\overline{V}\left(0\right)}{{\psi }_i}=O_p\left(\frac{1}{m}\right)$$

and is uniformly bounded under the regularity conditions given in Section 2, hence the last term in (A.11) is also an $o\left(1/m\right),$ which renders ${\text{mse}}_0$ unbiased up to the second order.

Appendix B

B.1 Comparison between REML and AR.YL using the scoring algorithm

The scoring algorithm could sometimes yield zero estimates for the likelihood of the AR.YL. Indeed, for data sets simulated under the model given in Section 5, with $m=45$ and ${\sigma }_v^2=1,$ the REML and AR.YL scoring algorithms yielded 28% and 26% zeros respectively. Figures B.1 to B.3 illustrate the why: the likelihoods correspond to a single population generated under the model with ${\sigma }_v^2=1$ for which ${\widehat{\sigma }}_{v\text{REML}}^2=0.$

Figure B.2 shows that the maximum value of the AR.YL likelihood is very near the border. The scoring algorithm may often miss the maximum and yield a zero value. Figure B.3 shows that the AM.LL likelihood has a maximum value that differentiates better from the border.

B.2 Treatment of zeros in the parametric bootstrap

For each estimate ${\widehat{\sigma }}_v^2={\widehat{\sigma }}_v^2\left(y^{\left(r\right)}\right),r=1,\dots 10K,$ and each method of variance estimation:

1. Generate a large number B of random area effects $v_i^{\left(b\right)}\stackrel{\text{i}\text{.i}\text{.d}\text{.}}{\sim }N\left(0,{\widehat{\sigma }}_v^2\right),b=1,\dots ,B,$ and generate, independently of $v_i^{\left(b\right)},$ sampling errors $e_i^{\left(b\right)}\stackrel{\text{i}\text{.i}\text{.d}\text{.}}{\sim }N\left(0,{\psi }_i\right),$ $i=1,\dots ,m,b=1,\dots ,B.$ Generate bootstrap data $y_i^{\left(b\right)}={\theta }_i^{\left(b\right)}+e_i^{\left(b\right)},{\theta }_i^{\left(b\right)}=x_i^{\prime }\widehat{\beta }+v_i^{\left(b\right)},i=1,\dots ,m.$ If ${\widehat{\sigma }}_{v\text{REML}}^2\left(y^{\left(r\right)}\right)=0,$ then generate $\left(y_i^{\left(b\right)},{\theta }_i^{\left(b\right)}\right),$ $b=1,\dots ,B,$ from the synthetic model (see also Rao and Molina 2015).
2. Fit the model to the bootstrap data and obtain ${\widehat{\sigma }}_v^{2\left(b\right)};$ for the MIX estimator calculate ${\widehat{\sigma }}_{v\text{MIX}}^{2\left(b\right)}={\widehat{\sigma }}_{v\text{REML}}^{2\left(b\right)}\text{ if }{\widehat{\sigma }}_v^{2\left(b\right)}$ is positive and ${\widehat{\sigma }}_{\text{vMIX}}^{\text{2}\left(b\right)}={\widehat{\sigma }}_{v\text{AM}}^{2\left(b\right)}$ otherwise.
3. Now obtain ${\widehat{\beta }}^{\left(b\right)},$ the corresponding EBLUP ${\widehat{\theta }}_i^{\left(b\right)},$ the bootstrap components $g_{1i}^{\left(b\right)}=g_{1i}\left({\widehat{\sigma }}_v^2\left(y^{\left(b\right)}\right)\right),g_{2i}^{\left(b\right)}=g_{2i}\left({\widehat{\sigma }}_v^2\left(y^{\left(b\right)}\right)\right)$ and ${\overline{g}}_{ji}^{\text{PB}}=B^{-1}{\displaystyle {\sum }_b{g}_{ji}^{\left(b\right)}},j=1,2.$
4. The Naive MSE bootstrap estimator is ${\text{mse}}_{\text{naive}}=B^{-1}{\displaystyle {\sum }_{b=1}^B{\left({\widehat{\theta }}_i^{\left(b\right)}-{\theta }_i^{\left(b\right)}\right)}^2}.$
5. The PB MSE estimator (which is adjusted for bias (Pfeffermann and Glickman 2004) is: ${\text{mse}}_{\text{PB}}\left({\widehat{\theta }}_i\right)=g_{1i}\left({\widehat{\sigma }}_v^2\right)+g_{2i}\left({\widehat{\sigma }}_v^2\right)-{\overline{g}}_{1i}^{\text{PB}}-{\overline{g}}_{2i}^{\text{PB}}+{\text{mse}}_{\text{naive}}.$
6. To calculate ${\text{ARB}}_{ℂ},$ average $\left({\text{mse}}_{\text{PB}}^{\left(r\right)}\left({\widehat{\theta }}_i\right)-MSE\left({\widehat{\theta }}_i\right)\right)/MSE\left({\widehat{\theta }}_i\right)$ over the populations with $\left(r\right)/{\widehat{\sigma }}_{v\text{REML}}^2\left(y^{\left(r\right)}\right)=0$ and do similarly with ${\text{ARB}}_{ℂ}$ of ${\text{mse}}_{\text{naive}}.$

References

Chen, S., and Lahiri, P. (2008). On mean squared prediction error estimation in small area estimation problems. Communications in Statistics-Theory and Methods, 37, 1792-1798.

Chen, S., and Lahiri, P. (2011). On the estimation of Mean Squared Prediction Error in small area estimation. Calcutta Statistical Association Bulletin, 63, (Special 7^th Triennial Proceedings Volume), Nos. 249-252.

Cressie, N. (1992). REML estimation in empirical Bayes smoothing of census undercount. Survey Methodology, 18, 1, 75-94.

Das, K., Jiang, J. and Rao, J.N.K. (2004). Mean squared error of empirical predictor. The Annals of Statistics, 32, 2, 818-840.

Datta, G., and Lahiri, P. (2000). A unified measure of uncertainty of estimated best linear unbiased predictors in small area estimation problems. Statistica Sinica, 10, 613-627.

Estevao, V. (2014). Grid optimization algorithm for maximum likelihood. Internal report, Statistical Research and Innovation Division (SRID), Statistics Canada.

Fay, R.E., and Herriot, R.A. (1979). Estimation of income from small places: An application of James-Stein Procedures to census data. Journal of the American Statistical Association, 74, 269-277.

Lahiri, P., and Li, H. (2009). Generalized maximum likelihood method in linear mixed models with an application in small area estimation. In Proceedings of the Federal Committee on Statistical Methodology Research Conference, available at http//www.fcsm.gov/events/papers2009.html.

Lahiri, P., and Pramanik, S. (2011). Discussion of “Estimating random effects via adjustment for density maximization” by C. Morris and R. Tang. Statistical Science, 26, 2, 291-295.

Li, H., and Lahiri, P. (2011). An adjusted maximum likelihood method for solving small area estimation problems. Journal of Multivariate Analysis, 101, 882-892.

Molina, I., Rao, J.N.K. and Datta, G.S. (2015). Small area estimation under a Fay-Herriot model with preliminary testing for the presence of random effects. Survey Methodology, 41, 1, 1-19.

Morris, C.N. (2006). Mixed model prediction and small area estimation (with discussions). Test, 15, 72-76.

Pfeffermann, D., and Glickman, H. (2004). Mean squarred error approximation in small area estimation by use of parametric and non-parametric bootstrap. Proceedings of the American Statistical Association, Section on Survey Research Methods, Alexandria, VA. 4167-78.

Rao, J.N.K. (2003). Small Area Estimation. New York: John Wiley & Sons, Inc.

Rao, J.N.K., and Molina, I. (2015). Small Area Estimation, second edition. New York: John Wiley & Sons, Inc.

Rubin-Bleuer, S., and Schiopu-Kratina, I. (2005). On the two-phase framework for joint model and design-based inference. The Annals of Statistics, 33, 6, 2789‑2810.

Rubin-Bleuer, S., and You, Y. (2012). A positive variance estimator for the Fay-Herriot small area model. SRID-2012-009E, Statistical Research and Innovation Division (SRID), Statistics Canada.

Rubin-Bleuer, S., Yung, W. and Landry, S. (2010). Adjusted maximum likelihood method for a small area model accounting for time and area effects. SRID-2010-006E, Statistical Research and Innovation Division (SRID), Statistics Canada.

Rubin-Bleuer, S., Yung, W. and Landry, S. (2011). Adjusted maximum likelihood method for a small area model accounting for time and area effects. Long abstract, Small Area Estimation, (SAE 20122) in Trier, Germany, International Statistical Institute Satellite Conference.

Rubin-Bleuer, S., Yung, W. and Landry, S. (2012). Variance Component Estimation through the Adjusted Maximum Likelihood Approach. Presentation at the Conference in Honour of the 75^th birthday of J.N.K. Rao Carleton University, May 2012, Ottawa.

Yoshimori, M., and Lahiri, P. (2014). A new adjusted maximum likelihood method for the Fay-Herriot small area model. Journal of Multivariate Analysis, 124, 281‑294.

Yuan, P. (2009). Comparison of SAE methods of variance estimation. Internal document, Statistical Research and Innovation Division (SRID), Statistics Canada.

Yuan, K.H., and Jennrich, R. (1998). Asymptotics of estimating equations under natural conditions. Journal of Multivariate Analysis, 65, 2, 245-260.

• Previous
• Next