3. Parameter estimation

Jae-kwang Kim, Seunghwan Park and Seo-young Kim

Previous | Next

Now, we discuss estimation of the model parameters in (2.3). The GLS estimator of $\beta =\left({\beta }_0\mathrm{,}{\beta }_1\right)$ can be obtained by minimizing

$$Q^{*}\left({\beta }_0\mathrm{,}{\beta }_1\right)={\displaystyle \sum _{h=1}^H}\frac{{\left({\overline{y}}_{1h}-{\beta }_0-{\beta }_1{\overline{x}}_h\right)}^2}{V\left({\overline{y}}_{1h}-{\beta }_0-{\beta }_1{\overline{x}}_h\right)}\mathrm{.}(3.1)$$

Since

$$V\left({\overline{y}}_{1h}-{\beta }_0-{\overline{x}}_h{\beta }_1\right)={\sigma }_{e\mathrm{,}h}^2+\left(-{\beta }_1\mathrm{,1}\right){\Sigma }_h{\left(-{\beta }_1\mathrm{,1}\right)}^{\prime }\mathrm{,}(3.2)$$

where ${\sigma }_{e\mathrm{,}h}^2=V\left({\overline{e}}_{1h}\right)$ and ${\Sigma }_h=V\left\{{\left(a_h\mathrm{,}{b}_h\right)}^{\prime }\right\},$ we can express

$$Q^{*}\left({\beta }_0\mathrm{,}{\beta }_1\right)={\displaystyle \sum _{h=1}^H}{w}_h\left({\beta }_1\right){\left({\overline{y}}_{1h}-{\beta }_0-{\beta }_1{\overline{x}}_h\right)}^2\mathrm{,}(3.3)$$

where $w_h\left({\beta }_1\right)={\left\{{\sigma }_{e\mathrm{,}h}^2+\left(-{\beta }_1\mathrm{,1}\right){\Sigma }_h{\left(-{\beta }_1\mathrm{,1}\right)}^{\prime }\right\}}^{-1}\mathrm{.}$ Now, by solving $\partial Q^{*}/\partial \beta =0,$ we have

$${\widehat{\beta }}_0={\overline{y}}_w-{\widehat{\beta }}_1{\overline{x}}_w(3.4)$$

and

$${\widehat{\beta }}_1=\frac{{\displaystyle \sum _{h=1}^H}{w}_h\left({\widehat{\beta }}_1\right)\left\{\left({\overline{x}}_h-{\overline{x}}_w\right)\left({\overline{y}}_{1h}-{\overline{y}}_{1w}\right)-C\left(a_h\mathrm{,}{b}_h\right)\right\}}{{\displaystyle \sum _{h=1}^H}{w}_h\left({\widehat{\beta }}_1\right)\left\{{\left({\overline{x}}_h-{\overline{x}}_w\right)}^2-V\left(a_h\right)\right\}}\mathrm{,}(3.5)$$

where

$$\left({\overline{x}}_w\mathrm{,}{\overline{y}}_w\right)={\left\{{\displaystyle \sum _{h=1}^H}{w}_h\left({\widehat{\beta }}_1\right)\right\}}^{-1}{\displaystyle \sum _{h=1}^H}{w}_h\left({\widehat{\beta }}_1\right)\left({\overline{x}}_h\mathrm{,}{\overline{y}}_h\right)\mathrm{.}$$

Note that the weight $w_h\left({\beta }_1\right)$  depends on ${\beta }_1.$  Thus, the solution (3.5) can be obtained by an iterative algorithm. Once ${\widehat{\beta }}_1$  is computed by (3.5), then ${\widehat{\beta }}_0$  is obtained by (3.4).

Now, we discuss the estimation of model variance ${\sigma }_{e\mathrm{,}h}^2.$  The simplest method is the Method of Moments (MOM). That is, we can use

$$E\left\{{\left({\overline{y}}_{1h}-{\beta }_0-{\overline{x}}_h{\beta }_1\right)}^2-{\beta }_1^{2}V\left(a_h\right)+2{\beta }_{1}C\left(a_h\mathrm{,}{b}_h\right)-V\left(b_h\right)\right\}={\sigma }_{e\mathrm{,}h}^2(3.6)$$

to obtain an unbiased estimator of ${\sigma }_{e\mathrm{,}h}^2.$  Under the nested error model in (2.4), we have ${\sigma }_{e\mathrm{,}h}^2={\sigma }_e^2$ and

$$E\left\{{\left({\overline{y}}_{1h}-{\beta }_0-{\overline{x}}_h{\beta }_1\right)}^2-{\beta }_1^{2}V\left(a_h\right)+2{\beta }_{1}C\left(a_h\mathrm{,}{b}_h\right)-V\left(b_h\right)\right\}={\sigma }_e^2\mathrm{.}(3.7)$$

Thus, similarly to Fuller (2009), the MOM estimator of ${\sigma }_e^2$  can be obtained by

$${\widehat{\sigma }}_e^2={\displaystyle \sum _{h=1}^H}{\kappa }_h\left\{{\left({\overline{y}}_{1h}-{\widehat{\beta }}_0-{\overline{x}}_h{\widehat{\beta }}_1\right)}^2-\left(-{\widehat{\beta }}_1\mathrm{,1}\right){\Sigma }_h\left(-{\widehat{\beta }}_1\mathrm{,1}\right)\right\}(3.8)$$

where

$${\kappa }_h\propto {\left\{{\widehat{\sigma }}_e^2+\left(-{\widehat{\beta }}_1\mathrm{,1}\right){\Sigma }_h\left(-{\widehat{\beta }}_1\mathrm{,1}\right)\right\}}^{-1}$$

and ${\displaystyle {\sum }_{h=1}^H}{\kappa }_h=1.$ Because ${\kappa }_h$  depends on ${\widehat{\sigma }}_e^2,$  the solution (3.8) can be obtained iteratively, using ${\widehat{\sigma }}_e^2=0$ as an initial value. Fay and Herriot (1979) used an alternative method which is based on the iterative solution to nonlinear equation:

$${\displaystyle \sum _{h=1}^H}\frac{{\left({\overline{y}}_{1h}-{\widehat{\beta }}_0-{\widehat{\beta }}_1{\overline{x}}_h\right)}^2}{{\sigma }_e^2+\left(-{\widehat{\beta }}_1\mathrm{,1}\right){\Sigma }_h{\left(-{\widehat{\beta }}_1\mathrm{,1}\right)}^{\prime }}=H-2.$$

Writing the above equation as $g\left({\sigma }_e^2\right)=H-2,$ a Newton-type method for $g\left(\theta \right)=0$ with $\theta ={\sigma }_e^2$ can be obtained by

$${\theta }^{\left(t+1\right)}={\theta }^{\left(t\right)}+\frac{1}{g^{\prime }\left({\theta }^{\left(t\right)}\right)}\left(H-2-g\left({\theta }^{\left(t\right)}\right)\right)(3.9)$$

where

$$g^{\prime }\left(\theta \right)=-{\displaystyle \sum _{h=1}^H}\frac{{\left({\overline{y}}_{1h}-{\widehat{\beta }}_0-{\widehat{\beta }}_1{\overline{x}}_h\right)}^2}{{\left\{\theta +\left(-{\widehat{\beta }}_1\mathrm{,1}\right){\Sigma }_h{\left(-{\widehat{\beta }}_1\mathrm{,1}\right)}^{\prime }\right\}}^2}\mathrm{.}$$

Assuming ${\sigma }_{e\mathrm{,}h}^2\equiv {\sigma }_e^2,$  we now describe the whole parameter estimation procedure as follows:

• Step 1 Compute the initial estimator of $\left({\beta }_0\mathrm{,}{\beta }_1\right)$  by setting ${\widehat{\sigma }}_e^2=0$ in (3.4) and (3.5).
• Step 2 Based on the current value of $\left({\widehat{\beta }}_0\mathrm{,}{\widehat{\beta }}_1\right),$  compute ${\widehat{\sigma }}_e^2$  using the iterative algorithm in (3.9).
• Step 3 Use the current value of ${\widehat{\sigma }}_e^2,$  compute the updated estimator of $\left({\beta }_0\mathrm{,}{\beta }_1\right)$  by (3.4) and (3.5).
• Step 4 Repeat [Step 2]-[Step 3] until convergence.

The proposed parameter estimation method estimates $\beta =({\beta }_0\mathrm{,}{\beta }_1)$ by the GLS and estimates ${\sigma }_e^2$  by the MOM iteratively. Note that the estimation of $\beta$  is based on data from all areas. If separate regression models are used, then the proposed parameter estimation method can be applied to the groups of areas. Instead of this separate iterative estimation method, we can also consider another method based on maximum likelihood estimation (MLE) under parametric distributional assumptions. See Carroll, Rupert, and Stefanski (1995) and Schafer (2001) for further discussion of MLE for parameters in the measurement error models.

Remark 2 If ${\sigma }_{e\mathrm{,}h}^2={\sigma }_e^2$ is not true, we can consider some alternative model such as

$${\overline{e}}_h\sim \left(\mathrm{0,}{\overline{X}}_h{\sigma }_e^2\right)\mathrm{.}(3.10)$$

To check whether model (3.10) holds, one can compute

$${\nu }_h={\left({\overline{y}}_{1h}-{\widehat{\beta }}_0-{\overline{x}}_h{\widehat{\beta }}_1\right)}^2-{\widehat{\beta }}_1^{2}V\left(a_h\right)+2{\widehat{\beta }}_1\widehat{C}\left(a_h\mathrm{,}{b}_h\right)-V\left(b_h\right)(3.11)$$

and plot ${\nu }_h$  on ${\overline{x}}_h.$  If the plot shows a linear relationship, then (3.10) can be treated as a reasonable model. Under model (3.10), we can obtain   ${\sigma }_e^2$  by a ratio method:

$${\widehat{\sigma }}_e^2=\frac{{\displaystyle \sum _{h=1}^H}{\kappa }_h{\nu }_h}{{\displaystyle \sum _{h=1}^H}{\kappa }_h{\widehat{\overline{X}}}_h}(3.12)$$

where

$${\kappa }_h\propto {\left\{{\widehat{\overline{X}}}_h{\widehat{\sigma }}_e^2+\left(-{\widehat{\beta }}_1\mathrm{,1}\right){\Sigma }_h\left(-{\widehat{\beta }}_1\mathrm{,1}\right)\right\}}^{-1}$$

with ${\displaystyle {\sum }_{h=1}^H}{\kappa }_h=1,$  ${\widehat{\overline{X}}}_h$  is defined in (2.9), and   ${\nu }_h$  is defined in (3.11). Because   ${\kappa }_h$  also depends on   ${\sigma }_e^2,$  the solution (3.12) can be obtained iteratively.

Remark 3 We can also consider a transformation   ${\overline{x}}_h^{*}=T\left({\overline{x}}_h\right)$ and   ${\overline{y}}_{1h}^{*}=T\left({\overline{y}}_{1h}\right)$ to improve the approximation to asymptotic normality. To check the departure from normality, plot   $n_{ha}\overline{V}\left({\overline{x}}_h\right)$  on   ${\overline{x}}_h.$  If the plot shows some structural relationship of   ${\overline{x}}_h$  then the normality assumption can be doubted. Now, consider the following transformation

$$T\left(x\right)=\log \left(x\right)\mathrm{.}(3.13)$$

Note that the asymptotic variance of   ${\overline{x}}_h^{*}=T\left({\overline{x}}_h\right)$ is equal to

$$V\left({\overline{x}}_h^{*}\right)\doteq \frac{1}{{\left({\overline{x}}_h\right)}^2}V\left({\overline{x}}_h\right)\mathrm{.}$$

Such transformation is a variance stabilizing transformation and is useful when we want to improve the approximation to normality.

Once the GLS estimator   ${\widehat{\overline{X}}}_h^{*}$ of   ${\overline{X}}_h^{*}$ is obtained, then we need to apply the inverse transformation to obtain the best estimator of   ${\overline{X}}_h=T^{-1}\left({\overline{X}}_h^{*}\right):=Q\left({\overline{X}}_h^{*}\right).$ Simply applying the inverse transformation will lead to biased estimation. To correct for the bias, we can use a second-order Taylor linearization. Using a Taylor expansion, we have

$$Q\left({\widehat{\overline{X}}}_h^{*}\right)\doteq Q\left({\overline{X}}_h^{*}\right)+Q^{\prime }\left({\overline{X}}_h^{*}\right)\left({\widehat{\overline{X}}}_h^{*}-{\overline{X}}_h\right)+\frac{1}{2}{Q}^{\prime }{}^{\prime }\left({\overline{X}}_h^{*}\right){\left({\widehat{\overline{X}}}_h^{*}-{\overline{X}}_h\right)}^2$$

and so, if we use   $Q\left({\widehat{\overline{X}}}_h^{*}\right)$ as an estimator for   ${\overline{X}}_h=Q\left({\overline{X}}_h^{*}\right),$ we have, ignoring the smaller order terms,

$$E\left\{Q\left({\widehat{\overline{X}}}_h^{*}\right)\right\}={\overline{X}}_h+\frac{1}{2}{Q}^{\prime }{}^{\prime }\left({\overline{X}}_h^{*}\right)V\left({\widehat{\overline{X}}}_h^{*}\right)\mathrm{.}$$

For the transformation in (3.13), we have   $Q\left({\overline{X}}_h^{*}\right)=\exp \left({\overline{X}}_h^{*}\right)$ and so   $Q^{\prime }{}^{\prime }\left({\overline{X}}_h^{*}\right)={\overline{X}}_h.$ Thus,   ${\widehat{\overline{X}}}_h=Q\left({\widehat{\overline{X}}}_h^{*}\right),$ we have

$$E\left({\widehat{\overline{X}}}_h\right)\cong {\overline{X}}_h+\frac{1}{2}{\overline{X}}_{h}V\left({\widehat{\overline{X}}}_h^{*}\right)$$

and the bias-corrected estimator of   ${\overline{X}}_h$  is

$${\widehat{\overline{X}}}_{h\mathrm{,}bc}=\frac{{\widehat{\overline{X}}}_h}{1+0.5V\left({\widehat{\overline{X}}}_h^{*}\right)}\mathrm{,}(3.14)$$

where   $V\left({\widehat{\overline{X}}}_h^{*}\right)$ is computed by the MSE estimation method which will be discussed in Section 4.

Previous | Next