Comparison of some positive variance estimators for the Fay-Herriot small area model 3. Review of REML and adjusted maximum likelihood methodsComparison of some positive variance estimators for the Fay-Herriot small area model 3. Review of REML and adjusted maximum likelihood methods

• Previous
• Next

3.1 REML method

We consider the combined Fay-Herriot model (2.3) with ${\sigma }_v^2> 0.$ The REML variance estimator of ${\sigma }_v^2$ is obtained by maximizing the residual likelihood function with respect to ${\sigma }_v^2:$

$$L_{\text{REML}}\left({\sigma }_v^2\right)\propto {\left|\left[{\displaystyle \sum _{i=1}^m{z}_i{z}_i^{\prime }}/\left({\sigma }_v^2+{\psi }_i\right)\right]\right|}^{-1/2}{{\displaystyle {\prod }_{i=1}^m\left({\sigma }_v^2+{\psi }_i\right)}}^{-1/2}\exp \left\{-\frac{1}{2}{y}^{\prime }Py\right\}$$

where $y=\left(y_1,\dots ,y_m\right)\prime ,P=V^{-1}-V^{-1}Z{\left(Z^{\prime }{V}^{-1}Z\right)}^{-1}{Z}^{\prime }{V}^{-1},$ $V=\text{Var}\left(y\right),$ and $Z={\left(z_1,\dots ,z_m\right)}^{\prime }.$ (Cressie 1992, Datta and Lahiri 2000 and Rao 2003, chapter 6). The REML variance estimator is given by:

$${\widehat{\sigma }}_{v\text{REML}}^2=\max \left({\tilde{\sigma }}_{v\text{REML}}^2,0\right),(3.1)$$

where ${\tilde{\sigma }}_{v\text{REML}}^2$ is the converging value of the REML algorithm. The asymptotic bias and variance of the REML estimator, up to the second order, are respectively given by:

$$\text{Bias}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)=o\left(\frac{1}{m}\right)\text{ and }V\left({\widehat{\sigma }}_{v\text{REML}}^2\right)=\frac{2}{\text{tr}\left(V^{-2}\right)}+o\left(\frac{1}{m}\right)\text{.}(3.2)$$

A second order unbiased estimator of the MSE of the EBLUP under REML variance estimation is given by (Datta and Lahiri 2000 and Chen and Lahiri 2008, 2011):

$$\text{mse}\left\{{\widehat{\theta }}_i\left({\widehat{\sigma }}_{v\text{REML}}^2\right)\right\}=\{\begin{array}{ll}{g}_{1i}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)+g_{2i}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)+2g_{3i}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)\hfill & \text{if }{\widehat{\sigma }}_{v\text{REML}}^2>0\hfill \\ g_{2i}\left(0\right)\hfill & \text{if }{\widehat{\sigma }}_{v\text{REML}}^2\text{= }0.\hfill \end{array}(3.3)$$

Remark 3.1. When ${\widehat{\sigma }}_v^2=0,$ the EBLUP reduces to the synthetic estimator. However, note that when

$${\widehat{\sigma }}_v^2=0,g_{1i}\left({\widehat{\sigma }}_v^2\right)=0,g_{2i}\left({\widehat{\sigma }}_v^2\right)=z_i^{\prime }{\left[{\displaystyle \sum _{i=1}^m{z}_i{z}_i^{\prime }}/{\psi }_i\right]}^{-1}{z}_i,$$

and $g_{3i}\left({\widehat{\sigma }}_v^2\right)=\overline{V}\left({\widehat{\sigma }}_v^2\right)/{\psi }_i>0,$ i.e., $\text{mse}\left\{{\widehat{\theta }}_i\left({\widehat{\sigma }}_v^2\right)\right\}$ is not a continuous function of ${\widehat{\sigma }}_v^2.$ We will see in the empirical study that when conditioning on $\left\{{\widehat{\sigma }}_v^2=0\right\},$ the MSE estimator in (3.3) has significant negative bias, unless the underlying signal to noise ratio ${\sigma }_v^2/{\psi }_i$ is negligible.

3.2 Adjusted maximum likelihood methods

The adjusted maximum likelihood variance estimators are derived from optimizing either the profile (AM) or the residual (AR) likelihood adjusted with the factor $h\left({\sigma }_v^2\right).$ As noted in the introduction, the AM.LL and AR.LL estimators use the adjustment factor $h_{\text{LL}}\left({\sigma }_v^2\right)={\sigma }_v^2,$ and the AM.YL and AR.YL estimators use the adjustment factor

$$h_{\text{YL}}\left({\sigma }_v^2\right)={\left\{\text{arc}\tan \left[{\displaystyle \sum _{i=1}^m{\sigma }_v^2}/\left({\sigma }_v^2+{\psi }_i\right)\right]\right\}}^{1/m}.$$

We denote by ${\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2$ and ${\widehat{\sigma }}_{v\text{AM}\text{.YL}}^2$ the variance estimators obtained by maximizing the adjusted profile likelihood functions, with respect to ${\sigma }_v^2:$

$$L_{\text{AM}\text{.*}}\left({\sigma }_v^2\right)\propto h\left({\sigma }_v^2\right)\cdot {{\displaystyle {\prod }_{i=1}^m\left({\sigma }_v^2+{\psi }_i\right)}}^{-1/2}\exp \left\{-\frac{1}{2}{y}^{\prime }Py\right\},(3.4)$$

where $h\left({\sigma }_v^2\right)=h_{\text{LL}}\left({\sigma }_v^2\right)$ and $h\left({\sigma }_v^2\right)=h_{\text{YL}}\left({\sigma }_v^2\right)$ for AM.LL and AM.YL respectively. The matrix $P$ is as in (3.1). The bias of the AM estimators up to the second order (denoted by $\approx )$ is:

$$B\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2\right)\approx \frac{\text{tr}\left\{P-V^{-1}\right\}+2/{\sigma }_v^2}{\text{tr}\left(V^{-2}\right)}=O\left(\frac{1}{m}\right)\text{and}B\left({\widehat{\sigma }}_{v\text{AM}\text{.YL}}^2\right)\approx \frac{\text{tr}\left\{P-V^{-1}\right\}}{\text{tr}\left(V^{-2}\right)}=O\left(\frac{1}{m}\right),(3.5)$$

(Li and Lahiri 2011 and Yoshimori and Lahiri 2014). The AR.LL and AR.YL variance estimators, denoted by ${\widehat{\sigma }}_{v\text{AR}\text{.LL}}^2$ and ${\widehat{\sigma }}_{v\text{AR}\text{.YL}}^2,$ are obtained by maximizing the adjusted residual (AR) likelihood functions with respect to ${\sigma }_v^2:$

$$L_{\text{AR}\text{.*}}\left({\sigma }_v^2\right)\propto h\left({\sigma }_v^2\right)\cdot {\left|{\displaystyle \sum _{i=1}^m{z}_i{z^{\prime }}_i}/\left({\sigma }_v^2+{\psi }_i\right)\right|}^{-1/2}{{\displaystyle {\prod }_{i=1}^m\left({\sigma }_v^2+{\psi }_i\right)}}^{-1/2}\exp \left\{-\frac{1}{2}{y}^{\prime }Py\right\}(3.6)$$

where $h\left({\sigma }_v^2\right)=h_{\text{LL}}\left({\sigma }_v^2\right)$ and $h\left({\sigma }_v^2\right)=h_{\text{YL}}\left({\sigma }_v^2\right)$ for AR.LL and AR.YL respectively and $P$ is as in (3.1), The asymptotic bias of the AR estimators are given, respectively by:

$$B\left({\widehat{\sigma }}_{v\text{AR}\text{.LL}}^2\right)\approx \frac{2/{\sigma }_v^2}{\text{tr}\left(V^{-2}\right)}=O\left(\frac{1}{m}\right)\text{ and }B\left({\widehat{\sigma }}_{v\text{AR}\text{.YL}}^2\right)=o\left(\frac{1}{m}\right).(3.7)$$

Under the regularity conditions given in Section 2 and ${\sigma }_v^2>0,$ the two LL and the two YL variance estimators exist and are $\sqrt{m}-$ consistent (Li and Lahiri 2011 and Yoshimori and Lahiri 2014). Lahiri and co-authors proposed the following MSE estimators:

$$\text{mse}\left\{{\widehat{\theta }}_i(\cdot )\right\}=g_{1i}(\cdot )+g_{2i}(\cdot )+2g_{3i}(\cdot )-{\psi }_i^2\cdot B(\cdot )/{\left(\cdot +{\psi }_i\right)}^2(3.8)$$

where the argument in $(\cdot )$ above is either ${\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2,$ ${\widehat{\sigma }}_{v\text{AR}\text{.LL}}^2$ or ${\widehat{\sigma }}_{v\text{AM}\text{.YL}}^2$ under AM.LL, AR.LL and AM.YL variance estimation respectively, and under ${\widehat{\sigma }}_{v\text{AR}\text{.YL}}^2:$

$$\text{mse}\left\{{\widehat{\theta }}_i\left({\widehat{\sigma }}_{v\text{AR}\text{.YL}}^2\right)\right\}=g_{1i}\left({\widehat{\sigma }}_{v\text{AR}\text{.YL}}^2\right)+g_{2i}\left({\widehat{\sigma }}_{v\text{AR}\text{.YL}}^2\right)+2g_{3i}\left({\widehat{\sigma }}_{v\text{AR}\text{.YL}}^2\right).(3.9)$$

Estimators (3.8) and (3.9) are unbiased up to the second order.

Remark 3.2. The sampling errors do not need to be normally distributed for the consistency and asymptotic normality of the LL and YL estimators (see, for example, Rubin-Bleuer et al. 2011).

3.3 Optimization algorithms

Given the data, the REML likelihood function may attain its maximum value at ${\sigma }_v^2=0,$ even when the true underlying value of ${\sigma }_v^2$ is positive. On the other hand, the LL and YL likelihoods always attain their maximum value at ${\sigma }_v^2>0.$ Yet, the YL residual likelihood is very close to the REML likelihood. Empirical studies show that the scoring algorithm under AR.YL yields ${\widehat{\sigma }}_{v\text{AR}\text{.YL}}^2=0$  in almost as large a percentage as under REML for data sets following a Fay-Herriot model with a small but non-zero true underlying variance. This happens when the scoring algorithm misses the positive maximum value of the AR.YL likelihood and outputs a zero value (see Appendix B for details). To avoid this problem, we use a grid method for optimization (Estevao 2014). In our study, we set the upper boundary of the search interval as $\text{1,000}\times {\sigma }_v^2,$ since we know ${\sigma }_v^2$ a priori. For applications with real data we suggest to obtain an initial estimate ${\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2$ by the method of scoring and set $\text{1,000}\times {\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2$ as the upper boundary. Then keep increasing the boundary until the variance estimate lies within the search interval.

• Previous
• Next