Small area estimation methods under cut-off sampling
Section 5. EBLUP under the nested error model

• Previous
• Next

Estimators described so far use only the outcome information coming from the domain. This means that, when the domain sample size $n_i$ is small, these estimators might be inefficient even in the absence of cut-off sampling. Small area (or indirect) estimation methods are designed to reduce the variance by increasing the effective sample size; see Rao and Molina (2015) for a comprehensive account of small area estimation methods. In this section, we focus on model-based methods, which provide estimators with good properties under the distribution induced by the model. Since the model-based properties are well known, we wish to analyze whether the estimators have good properties under the sampling-replication mechanism, which does not assume that the model actually holds.

We consider a very popular unit level model introduced by Battese, Harter and Fuller (1988) and often called nested error model. Similarly as for model $m_2$ in (4.16), this model assumes a constant linear regression for all the population units, but allows for unexplained heterogeneity between the domains by including random domain effects $u_i$ apart from model errors $e_{ij}.$ This model, denoted model $m_3,$ assumes

$$\begin{array}{ll}{y}_{ij}\hfill & =x_{ij}^{\prime }\beta +u_i+e_{ij}\mathrm{,}{u}_i\stackrel{\text{iid}}{{\displaystyle \sim }}N\left(\mathrm{0,}{\sigma }_u^2\right)\mathrm{,}\hfill \\ e_{ij}\hfill & \stackrel{\text{iid}}{{\displaystyle \sim }}N\left(\mathrm{0,}{\sigma }_e^2\right)\mathrm{,}j\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}{N}_i\mathrm{,}i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}m\mathrm{,}(5.1)\hfill \end{array}$$

where area effects $u_i$ and errors $e_{ij}$ are all mutually independent. The vectors $\beta$ and $\theta \mathrm{<mo>=</mo>}{\left({\sigma }_u^2\mathrm{,}{\sigma }_e^2\right)}^{\prime }$ are unknown. Setting ${\sigma }_u^2\mathrm{<mo>=</mo>}0$ in (5.1), we obtain model $m_2$ given in (4.16). If $y_i\mathrm{<mo>=</mo>}{\left(y_{i1}\mathrm{,}\dots \mathrm{,}{y}_{i{N}_i}\right)}^{\prime }$ denotes the vector of outcomes for domain $i$ and $X_i\mathrm{<mo>=</mo>}{\left(x_{i1}\mathrm{,}\dots \mathrm{,}{x}_{i{N}_i}\right)}^{\prime }$ the corresponding design matrix, the model in matrix notation reads

$$y_i\stackrel{\text{ind}}{{\displaystyle \sim }}N\left(X_i\beta \mathrm{,}{V}_i\right)\mathrm{,}{V}_i\mathrm{<mo>=</mo>}{\sigma }_u^2{1}_{N_i}{1}_{N_i}^{\prime }+{\sigma }_e^2{I}_{N_i}\text{}\mathrm{,}i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}m\mathrm{,}(5.2)$$

where $1_k$ denotes a vector of ones of size $k$ and $I_k$ is the $k\times k$ identity matrix.

We consider linear domain parameters defined as $H_i\mathrm{<mo>=</mo>}{b}_i^{\prime }{y}_i,$ where $b_i$ is a non-stochastic vector of known elements. The domain mean $H_i\mathrm{<mo>=</mo>}{\overline{Y}}_i\mathrm{<mo>=</mo>}{N}_i^{-1}{\displaystyle {\sum }_{j\mathrm{<mo>=</mo>\; 1}}^{N_i}}{y}_{ij}$ is obtained with $b_i\mathrm{<mo>=</mo>}{N}_i^{-1}{1}_{N_i}.$

A sample $s_i$ is supposed to be drawn from the set of included units in domain $i,$ that is, $s_i\subset U_{iI}.$ We denote by $r_i\mathrm{<mo>=</mo>}\left(U_{iI}-s_i\right)\cup U_{iE}$ the set of non-sampled units from domain $U_i,$ which includes those non-sampled units from $U_{iI}$ and all the units in $U_{iE}.$ Note that $U_i\mathrm{<mo>=</mo>}{s}_i\cup r_i\mathrm{<mo>=</mo>}{U}_{iI}\cup U_{iE}.$ Then, the overall sample $s$ is composed of the samples $s_i$ drawn from the sets of included units in each area $U_{iI},$ $i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}m,$ that is, $s\mathrm{<mo>=</mo>}{s}_1\cup \cdots \cup s_m.$

We decompose the domain vector $y_i$ and the design and covariance matrices $X_i$ and $V_i$ into the corresponding subvectors and submatrices for sample and out-of-sample units, indicated with subscripts $s$ and $r$ respectively, as follows

$$y\mathrm{<mo>=</mo>}\left(\begin{array}{c}{y}_{is}\\ y_{ir}\end{array}\right)\mathrm{,}X\mathrm{<mo>=</mo>}\left(\begin{array}{c}{X}_{is}\\ X_{ir}\end{array}\right)\mathrm{,}{V}_i\mathrm{<mo>=</mo>}\left(\begin{array}{cc}{V}_{is}& V_{isr}\\ V_{irs}& V_{ir}\end{array}\right)\mathrm{.}$$

The linear parameter $H_i\mathrm{<mo>=</mo>}{b}_i^{\prime }{y}_i$ can then be expressed as $H_i\mathrm{<mo>=</mo>}{b}_{is}^{\prime }{y}_{is}+b_{ir}^{\prime }{y}_{ir}.$ Under model (5.1), the best linear unbiased predictor (BLUP) of $H$ is the model-unbiased linear function of the sample data ${\widehat{H}}_i\mathrm{<mo>=</mo>}{\alpha }_{is}^{\prime }{y}_{is},$ which minimizes the model mean squared error (MSE), ${\text{MSE}}_{m_3}\left({\widehat{H}}_i\right)\mathrm{<mo>=</mo>}{E}_{m_3}{\left({\widehat{H}}_i-H_i\right)}^2.$ The BLUP of $H_i\mathrm{<mo>=</mo>}{b}_{is}^{\prime }{y}_{is}+b_{ir}^{\prime }{y}_{ir}$ is then

$${\widehat{H}}_i^{\text{BLUP}}\left(\text{}\theta \text{}\right)\mathrm{<mo>=</mo>}{b}_{is}^{\prime }{y}_{is}+b_{ir}^{\prime }\left[X_{ir}{\tilde{\beta }}_s+V_{irs}{V}_{is}^{-1}\left(y_{is}-X_{is}^{\prime }{\tilde{\beta }}_s\right)\right]\mathrm{,}(5.3)$$

where ${\tilde{\beta }}_s$ is the weighted least squares estimator of $\beta ,$ given by

$${\tilde{\beta }}_s\mathrm{<mo>=</mo>}{\tilde{\beta }}_s\left(\text{}\theta \text{}\right)\mathrm{<mo>=</mo>}{\left({\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^m}{X}_{is}^{\prime }{V}_{is}^{-1}{X}_{is}\right)}^{-1}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^m}{X}_{is}^{\prime }{V}_{is}^{-1}{y}_{is}\mathrm{.}(5.4)$$

The BLUP of $H_i$ given in (5.3) depends on the true values of the variance components $\theta \mathrm{<mo>=</mo>}{\left({\sigma }_u^2\mathrm{,}{\sigma }_e^2\right)}^{\prime },$ which are typically unknown. Replacing them by corresponding model-consistent estimators $\widehat{\theta }\mathrm{<mo>=</mo>}{\left({\widehat{\sigma }}_u^2\mathrm{,}{\widehat{\sigma }}_e^2\right)}^{\prime },$ we obtain the so-called empirical BLUP (EBLUP), denoted ${\widehat{H}}_i^{\text{EBLUP}}\mathrm{<mo>=</mo>}{\widehat{H}}_i^{\text{BLUP}}\mathrm{<mo\; stretchy="false">(</mo>}\widehat{\theta }\mathrm{<mo\; stretchy="false">)</mo>\; .}$

If the domain sampling fraction, $n_i/N_i,$ is negligible, the BLUP of ${\overline{Y}}_i$ may be expressed as the weighted average

$${\widehat{\overline{Y}}}_i^{\text{BLUP}}\cong {\gamma }_{is}\left[{\overline{y}}_{is}+{\left({\overline{X}}_i-{\overline{x}}_{is}\right)}^{\prime }{\tilde{\beta }}_s\right]+\left(1-{\gamma }_{is}\right){\overline{X}}_i^{\prime }{\tilde{\beta }}_s\mathrm{,}(5.5)$$

where ${\gamma }_{is}\mathrm{<mo>=</mo>}{\sigma }_u^2/\mathrm{<mo\; stretchy="false">(</mo>}{\sigma }_u^2+{\sigma }_e^2/n_i\mathrm{<mo\; stretchy="false">)</mo>}$ is in the $\left(\mathrm{0,}1\right)$ interval and tends to 1 as $n_i\to \infty$ (Rao and Molina, 2015). Thus, for domains with large sample size $n_i,$ ${\widehat{\overline{Y}}}_i^{\text{BLUP}}$ approaches the survey regression estimator ${\overline{y}}_{is}+{\left({\overline{X}}_i-{\overline{x}}_{is}\right)}^{\prime }{\tilde{\beta }}_s,$ whereas for domains with small sample size $n_i,$ ${\widehat{\overline{Y}}}_i^{\text{BLUP}}$ borrows strength from the other domains by approaching the regression-synthetic estimator ${\overline{X}}_i^{\prime }{\tilde{\beta }}_s.$ Replacing the variance components in $\theta \mathrm{<mo>=</mo>}{\left({\sigma }_u^2\mathrm{,}{\sigma }_e^2\right)}^{\prime }$ by consistent estimators $\widehat{\theta }\mathrm{<mo>=</mo>}{\left({\widehat{\sigma }}_u^2\mathrm{,}{\widehat{\sigma }}_e^2\right)}^{\prime }$ in the BLUP, denoting ${\widehat{\gamma }}_{is}\mathrm{<mo>=</mo>}{\widehat{\sigma }}_u^2/\left({\widehat{\sigma }}_u^2+{\widehat{\sigma }}_e^2/n_i\right)$ and ${\widehat{\beta }}_s\mathrm{<mo>=</mo>}{\tilde{\beta }}_s\left(\text{}\widehat{\theta }\text{}\right),$ we obtain the EBLUP of ${\overline{Y}}_i,$ given by

$${\widehat{\overline{Y}}}_i^{\text{EBLUP}}\cong {\widehat{\gamma }}_{is}\left[{\overline{y}}_{is}+{\left({\overline{X}}_i-{\overline{x}}_{is}\right)}^{\prime }{\widehat{\beta }}_s\right]+\left(1-{\widehat{\gamma }}_{is}\right){\overline{X}}_i^{\prime }{\widehat{\beta }}_s\mathrm{.}(5.6)$$

The BLUP is unbiased and optimal under model $m_3$ in the sense of minimizing the MSE under that model. We now study its design properties, which do not assume that the model is correct and hence account for bias under model departures. To that end, we consider the census regression parameter for the included units, defined as $B_I\mathrm{<mo>=</mo>}{\left({\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^m}{X}_{iI}^{\prime }{V}_{iI}^{-1}{X}_{iI}\right)}^{-1}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^m}{X}_{iI}^{\prime }{V}_{iI}^{-1}{y}_{iI},$ where $y_{iI},$ $X_{iI}$ and $V_{iI}$ are the corresponding sub-vector and sub-matrices of $y_i,$ $X_i$ and $V_i,$ for the included units $\left(j\in U_{iI}\right).$ Again, we consider the theoretical version of the BLUP defined in terms of $B_I,$

$${\tilde{\overline{Y}}}_i^{\text{BLUP}}\mathrm{<mo>=</mo>}{\gamma }_{is}\left[{\overline{y}}_{is}+{\left({\overline{X}}_i-{\overline{x}}_{is}\right)}^{\prime }\text{}{B}_I\right]+\left(1-{\gamma }_{is}\right){\overline{X}}_i^{\prime }{B}_I\mathrm{.}(5.7)$$

If each sample $s_i$ is drawn from the corresponding domain $U_{iI}$ by simple random sampling without replacement (SRSWOR), then $E_{\pi }\left({\overline{y}}_{is}\right)\mathrm{<mo>=</mo>}{\overline{Y}}_{iI}$ and $E_{\pi }\left({\overline{x}}_{is}\right)\mathrm{<mo>=</mo>}{\overline{X}}_{iI}.$ Using these facts, it is easy to calculate the design-bias of ${\tilde{\overline{Y}}}_i^{\text{BLUP}}$ under SRSWOR, which is given by

$$B_{\pi }\left({\tilde{\overline{Y}}}_i^{\text{BLUP}}\right)\mathrm{<mo>=</mo>}{\gamma }_{is}\frac{N_{iE}}{N_{iI}}\left[\left({\overline{Y}}_i-{\overline{X}}_i^{\prime }{B}_I\right)-\left({\overline{Y}}_{iE}-{\overline{X}}_{iE}^{\prime }{B}_I\right)\right]+\left(1-{\gamma }_{is}\right)\left({\overline{X}}_i^{\prime }{B}_I-{\overline{Y}}_i\right)\mathrm{.}(5.8)$$

This bias will be small if (5.1) holds for the whole population, in which case $E_{m_3}\text{}\left({\overline{Y}}_i\right)\mathrm{<mo>=</mo>}{\overline{X}}_i\beta$ and $E_{m_3}\text{}\left({\overline{Y}}_{iE}\right)\mathrm{<mo>=</mo>}{\overline{X}}_{iE}\beta .$ Using these results when taking expectation under model $m_3$ in (5.8), we get $B_{m_3\text{}\mathrm{,}\pi }\left({\tilde{\overline{Y}}}_i^{\text{BLUP}}\right)\mathrm{<mo>=</mo>}0.$ In fact, the same result also holds under model $m_2.$

Concerning variance, if $s_i$ is obtained by SRSWOR within $U_{iI},$ the design-variance of the theoretical BLUP estimator is given by

$$V_{\pi }\left({\tilde{\overline{Y}}}_i^{\text{BLUP}}\right)\mathrm{<mo>=</mo>}{\gamma }_{is}^2{V}_{\pi }\left({\overline{y}}_{is}-{\overline{x}}_{is}{B}_I\right)\mathrm{<mo>=</mo>}\frac{{\gamma }_{is}^2}{N_i^2}{V}_{\pi }\left({\widehat{Y}}_i-{\widehat{X}}_i^{\prime }{B}_I\right)\mathrm{.}$$

Hence, if the census least squared (LS) regression lines for the domains from model (4.10) are similar to the national census weighted least squared (WLS) regression line from model (5.1), that is, if $B_I\approx B_{iI},$ then the variance of the BLUP for ${\overline{Y}}_i$ reduces to that of the LCAL estimator of ${\overline{Y}}_i$ obtained from (4.17), multiplied by the factor ${\gamma }_{is}^2\in \left(0,1\right).$

Under more general sampling designs within $U_{iI},$ we consider the pseudo-EBLUP of ${\overline{Y}}_i$ proposed by You and Rao (2002) instead of the EBLUP. Defining the analogous theoretical estimator that uses the weighted sample means ${\overline{y}}_{iw}\mathrm{<mo>=</mo>}{\left({\displaystyle {\sum }_{j\in s_i}}{w}_{j|i}\right)}^{-1}{\displaystyle {\sum }_{j\in s_i}}{w}_{j|i}{y}_{ij}$ and ${\overline{x}}_{iw}\mathrm{<mo>=</mo>}{\left({\displaystyle {\sum }_{j\in s_i}}{w}_{j|i}\right)}^{-1}{\displaystyle {\sum }_{j\in s_i}}{w}_{j|i}{x}_{ij}$ instead or the unweighted ones ${\overline{y}}_{is}$ and ${\overline{x}}_{is}$ in (5.7), we obtain the same expressions for the design bias and variance, with ${\gamma }_{is}$ changed to ${\gamma }_{iw}\mathrm{<mo>=</mo>}{\sigma }_u^2/\left({\sigma }_u^2+{\sigma }_e^2{\delta }_{iw}\right),$ for ${\delta }_{iw}\mathrm{<mo>=</mo>}{\left({\displaystyle {\sum }_{j\in s_i}}{w}_{j|i}\right)}^{-2}{\displaystyle {\sum }_{j\in s_i}}{w}_{j|i}^2.$

• Previous
• Article main page
• Next