Local polynomial estimation for a small area mean under informative sampling
Section 3. The local polynomial estimator

• Previous
• Next

3.1  The estimation of a small area mean

The objective is to estimate the mean ${\overline{Y}}_i$ for small area $U_i$ for $i=1,\dots ,M.$ Splitting the population $U_i$ into observed units in the sample, $s_i$ of size $n_i,$ and non-observed units in the non-sampled portion, ${\overline{s}}_i=U_i/s_i$ of size $N_i-n_i,$ we can express ${\overline{Y}}_i$ as

$${\overline{Y}}_i=\frac{1}{N_i}\left({\displaystyle \sum _{j\in s_i}{y}_{ij}}+{\displaystyle \sum _{j\in {\overline{s}}_i}{y}_{ij}}\right).(3.1)$$

Given that we do not know the $y$  values for the non-observed units in sets ${\overline{s}}_i$ for $i=1,\dots ,M,$ we need to estimate them. Denoting as ${\widehat{y}}_{ij}$ the estimator of $y_{ij}$ for such units, the resulting estimator of the mean ${\overline{Y}}_i$ is

$${\widehat{\overline{Y}}}_i=\frac{1}{N_i}\left({\displaystyle \sum _{j\in s_i}{y}_{ij}}+{\displaystyle \sum _{j\in {\overline{s}}_i}{\widehat{y}}_{ij}}\right).(3.2)$$

We obtain estimators ${\widehat{y}}_{ij}$ of $y_{ij},$ for $j\in {\overline{s}}_i,$ based on an augmented model that includes an unknown smooth function of the selection probabilities $p_{j|i},$ denoted $m_0\left(p_{j|i}\right).$ The proposed augmented semi-parametric sample model is given by

$$y_{ij}={\tilde{x}}_{ij}^T{\beta }_1+m_0\left(p_{j|i}\right)+v_{1i}+e_{1ij},j=1,\dots ,n_i;i=1,\dots ,M,(3.3)$$

where $v_{1i}\stackrel{\text{iid}}{{\displaystyle \sim }}N\left(0,{\sigma }_{1v}^2\right)$ and independent of $e_{1ij}\stackrel{\text{iid}}{{\displaystyle \sim }}N\left(0,{\sigma }_{1e}^2\right).$ The vector ${\tilde{x}}_{ij}={\left(x_{ij1},\dots ,x_{ijp}\right)}^T$ in model (3.3) represents the covariates $x_{ij}$ without a constant (i.e., the intercept) and ${\beta }_1={\left({\beta }_{11},\dots ,{\beta }_{1p}\right)}^T$ a vector of fixed effects. Model (3.3) is semi-parametric as the response variable $y_{ij}$ depends linearly on the vector of auxiliary variables, ${\tilde{x}}_{ij},$ and the probability of selection $p_{j|i}$ enters non-parametrically through the smooth function $m_0(\cdot ).$

We assume that model (3.3) has a similar covariance structure with the one associated with model (1.2): the small area effects $v_{1i}$ and random errors $e_{1ij}$ are i.i.d., normally distributed and independently of one another. However, the semi-parametric model (3.3) is more flexible than the parametric model (1.2), as it does not force the function $m_0\left(p_{j|i}\right)$ to be of a specific form. There is a disadvantage to this set-up. Since model (3.3) is not a linear mixed model, the general EBLUP theory given in Section 2 cannot be applied directly to obtain estimators of $m_0\left(p_{j|i}\right),$ ${\beta }_1$ and $v_{1i}.$ Consequently, we propose to estimate (3.3) by combining the EBLUP theory for linear mixed models and the local polynomial technique (Fan and Gijbels, 1996).

We estimate (3.3) in three steps. In the first step, we obtain estimates of $m_0\left(p_{j|i}\right),$ ${\widehat{m}}_0\left(p_{j|i}\right),j=1,\dots ,N_i,i=1,\dots ,M,$ for all units in the population. These estimates are local in character as they are based on the local polynomial technique. Estimates ${\widehat{m}}_0\left(p_{j|i}\right),j\in s_i$ for the observed units are then used in the second step to obtain global estimators of ${\beta }_1$ and $v_{1i},i=1,\dots ,M.$ We denote these estimators as ${\widehat{\beta }}_{\text{glo},1}$ and ${\widehat{v}}_{\text{glo},1i},i=1,\dots ,M.$ Finally, in the third step, we use the local estimators ${\widehat{m}}_0\left(p_{j|i}\right)$ for the unobserved units, obtained in the first step, and the global estimators ${\widehat{\beta }}_{\text{glo},1}$ and ${\widehat{v}}_{\text{glo},1i}$ obtained in the second step, to estimate $y_{ij}$ for $j\in {\overline{s}}_i$ and $i=1,\dots ,M.$ The resulting estimators of $y_{ij},$ denoted as ${\widehat{y}}_{ij},$ are

$${\widehat{y}}_{ij}={\tilde{x}}_{ij}^T{\widehat{\beta }}_{\text{glo},1}+{\widehat{m}}_0\left(p_{j|i}\right)+{\widehat{v}}_{\text{glo},1i},j\in {\overline{s}}_i.(3.4)$$

The ${\widehat{y}}_{ij}’\text{s}$ are incorporated into equation (3.2) to obtain the estimator of the small area mean ${\widehat{\overline{Y}}}_i.$

We now proceed to describe the first step in more detail. Following Ruppert and Matteson (2015), we estimate the values of the unknown function $m_0\left(p_{l|k}\right)$ for all units $l\in U_k$ and small areas $k,$  with $k=1,\dots ,M,$ by using local polynomial regression. Local polynomial regression is based on the principle that a smooth function can be approximated locally by a low-degree polynomial. We approximate $m_0\left(p_{j|i}\right)$ in model (3.3) by a $q^{\text{th}}$ -degree polynomial, say $m_1\left(p_{j|i}\right),$ using a Taylor expansion around $p_{l|k}.$ The approximation is given by

$$m_1\left(p_{j|i}\right)=m_0\left(p_{l|k}\right)+{\displaystyle \sum _{a=1}^q\frac{1}{a!}{m}_0{\left(p_{l|k}\right)}^{\left(a\right)}{\left(p_{j|i}-p_{l|k}\right)}^a},j\in s_i;i=1,\dots ,M,(3.5)$$

where $m_0{\left(p_{l|k}\right)}^{\left(a\right)}$ is the $a^{\text{th}}$ derivative of $m_0\left(p_{j|i}\right)$ evaluated at $p_{l|k}.$ The function $m_1\left(p_{j|i}\right)$ depends on $l\in U_k,$ but we suppress this dependence to simplify the notation.

For each point $p_{l|k},$ $l\in U_k;k=1,\dots ,M,$ in model (3.3) we replace $m_0\left(p_{j|i}\right)$ by its approximation $m_1\left(p_{j|i}\right)$ given by (3.5). The resulting model is given by

$$y_{ij}={\tilde{x}}_{ij}^T{\beta }_1+m_0\left(p_{l|k}\right)+{\displaystyle \sum _{a=1}^q\frac{1}{a!}{m}_0{\left(p_{l|k}\right)}^{\left(a\right)}{\left(p_{j|i}-p_{l|k}\right)}^a}+v_{1i}+e_{1ij}\text{,}j\in s_i;i=1,\dots ,M.(3.6)$$

Model (3.6) is an approximate local model for (3.3) depending on the point $l\in U_k$ of the population. Estimates of ${\beta }_1$ and $v_{1i}$ based on (3.6) will be denoted by ${\widehat{\beta }}_{\text{loc},1}$ and ${\widehat{v}}_{\text{loc},1i}.$ Notice that (3.6) allows the estimation of $m_0\left(p_{l|k}\right),$ the value of the smooth function $m_0(\cdot )$ at a point $p_{l|k}.$ We express (3.6) as

$$y_{ij}={\tilde{x}}_{ij}^T{\beta }_1+u_0+{\displaystyle \sum _{a=1}^q{u}_a{\left(p_{j|i}-p_{l|k}\right)}^a}+v_{1i}+e_{1ij}\text{:}j\in s_i;i=1,\dots ,M,(3.7)$$

where $u_a=m_0{\left(p_{l|k}\right)}^{\left(a\right)}/a!$ for $a=0,\dots ,q.$ Model (3.7) is a linear mixed model with fixed parameters $\left({\beta }_1,u_0,\dots ,u_q\right)$ and random small area effects $v_{1i},i=1,\dots ,M.$

Let ${\widehat{u}}_0$ be an estimator of $u_0$ obtained by fitting model (3.7). An approximate estimator of $m_0\left(p_{l|k}\right)=u_0$ is given by ${\widehat{m}}_0\left(p_{l|k}\right)={\widehat{u}}_0.$ Since we require estimators of $m_0\left(p_{l|k}\right)$ for $l\in U_k$ and $k=1,\dots ,M,$ we use $N={\displaystyle {\sum }_{i=1}^M{N}_i}$ models (3.7). As pointed out by an Associate Editor, if $N$  is large, estimating the values of $m_0(\cdot )$ for all points in the population can be computationally intensive.

It is more convenient to work with matrix notation. To this end, we define $y_i={\left(y_{i1},\dots ,y_{i{n}_i}\right)}^T,$ ${\tilde{X}}_i={\left({\tilde{x}}_{i1}^T,\dots ,{\tilde{x}}_{i{n}_i}^T\right)}^T,$ $m_{0,i}={\left(m_0\left(p_{1|i}\right),\dots ,m_0\left(p_{n_i|i}\right)\right)}^T,$ $v_1={\left(v_{11},\dots ,v_{1M}\right)}^T$ and $e_{1i}={\left(e_{1i1},\dots ,e_{1i{n}_i}\right)}^T.$ Model (3.3) can be expressed in a matrix form by stacking the observations, and the resulting equation is

$$y=\tilde{X}{\beta }_1+m_0+Z{v}_1+e_1,(3.8)$$

where $y={\text{col}}_{\text{1}\le i\le M}\left(y_i\right),$ $\tilde{X}={\text{col}}_{1\le i\le M}\left({\tilde{X}}_i\right),$ $m_0={\text{col}}_{\text{1}\le i\le M}\left(m_{0,i}\right),$ $Z={\text{diag}}_{1\le i\le M}\left\{1_{n_i}\right\}$ and $e_1={\text{col}}_{1\le i\le M}\left(e_{1i}\right).$

For unit $l$ in small area $U_k,$ we define the $n\times \left(q+1\right)$ matrix:

$$Q=\left(\begin{array}{cccc}1& \left(p_{1|1}-p_{l|k}\right)& \cdots & {\left(p_{1|1}-p_{l|k}\right)}^q\\ ⋮& ⋮& \cdots & ⋮\\ 1& \left(p_{n_M|M}-p_{l|k}\right)& \cdots & {\left(p_{n_M|M}-p_{l|k}\right)}^q\end{array}\right),$$

where $n={\displaystyle {\sum }_{i=1}^M{n}_i}$ is the total sample size. Let $u={\left(m_0\left(p_{l|k}\right),m_0^{\left(1\right)}\left(p_{l|k}\right)/1!,\dots ,m_0^{\left(q\right)}\left(p_{l|k}\right)/q!\right)}^T$ represent the vector of derivatives of the function $m_0(\cdot )$ evaluated at $p_{l|k}.$ The terms $Q$ and $u$ depend on the unit $l\in U_k$ where the localization is realized. We omitted their dependence on the unit $l$ from small area $U_k$ in order not to burden the notation. We define vector $m_1$ obtained by stacking the $n$ values of the function $m_1(\cdot )$ defined by (3.5). That is, $m_1={\text{col}}_{\text{1}\le i\le M}\left(m_{1,i}\right)$ with $m_{1,i}={\left(m_1\left(p_{1|i}\right),\dots ,m_1\left(p_{n_i|i}\right)\right)}^T.$ This allows to approximate $m_0$ by $m_0\approx m_1.$ The vector $m_1$ is given by $m_1=Qu.$ It then follows that an approximation to (3.8) in a neighbourhood of $l\in U_k$ is

$$y=\tilde{X}{\beta }_1+Qu+Z{v}_1+e_1.(3.9)$$

Equations (3.8) and (3.9) are the matrix form equivalents of equations (3.3) and (3.7), respectively. The matrix $\tilde{X}$ in (3.9) does not include the constant term that represents the intercept, because this term is already included in $Q.$ Equation (3.9) is a standard linear mixed effects model with fixed parameters ${\beta }_{\text{fixed}}={\left({\beta }_1^T\text{},u^T\right)}^T$ and random small area effects $v_1.$ We denote by $V\left(v_1\right)=G={\sigma }_{1v}^2{I}_M,$ $V\left(e_{1i}\right)=R_i={\sigma }_{1e}^2{I}_{n_i}$ and $V\left(e_1\right)=R={\text{diag}}_{1\le i\le M}\left\{R_i\right\}$ as the respective covariance matrices of $v_1,$ $e_{1i}$ and $e_1.$ The covariance matrix of $y_i$ is given by $V\left(y_i\right)=V_i={\sigma }_{1v}^2{J}_{n_i}+{\sigma }_{1e}^2{I}_{n_i}.$ The matrices $I_M$ and $I_{n_i}$ are the identity matrices of order $M$ and $n_i$ respectively, whereas $J_{n_i}$ is the square matrix of order $n_i$ with all its elements equal to 1. It follows that $V\left(y\right)=V={\text{diag}}_{1\le i\le M}\left\{V_i\right\}.$

Assume that $V$ is known and that $v_1$ and $e_1$ are normally distributed. Using classical EBLUP theory, estimators of ${\beta }_{\text{fixed}}$ and $v_1$ can be obtained by minimizing

$$\Phi ={\left(y-\tilde{X}{\beta }_1-Qu-Z{v}_1\right)}^T{R}^{-1}\left(y-\tilde{X}{\beta }_1-Qu-Z{v}_1\right)+v_1^T{G}^{-1}{v}_1.$$

Note that all the observations that are included in $\Phi$ are equally weighted. However, we need to modify $\Phi$ to be in line with how local polynomial estimation is carried out. To this end, referring back to equation (3.7), we estimate its parameters by associating kernel weights $K\left(\left(p_{j|i}-p_{l|k}\right)/h\right)/h$ to each sampled unit $j\in s_i;i=1,\dots ,M.$ These kernel weights are chosen so as to give a larger weight to the sample points that are close to $l\in U_k,$ and a smaller weight to those that are further away. The weight $K(\cdot )$ is a probability density function and $h$ is a bandwidth controlling the size of the local neighbourhood. We explain in Section 3.2 how an optimal bandwidth can be obtained. Let $W$ be the $n\times n$ diagonal matrix of kernel weights given by

$$W=\underset{\begin{array}{l}1\le j\le n_i\\ 1\le i\le M\end{array}}{\text{diag}}\left\{\frac{1}{h}K\left(\frac{p_{j|i}-p_{l|k}}{h}\right)\right\}.$$

The matrix $W$ depends on unit $l$ from small area $U_k$ and the bandwidth $h.$ We do not include the subscripts $l\in U_k$ and $h$ in the definition of the matrix $W,$ in order not to burden notation. Following Wu and Zhang (2002), the incorporation of the kernel weights in $\Phi$ lead us to minimize ${\Phi }_W$ where

$${\Phi }_W={\left(y-\tilde{X}{\beta }_1-Qu-Z{v}_1\right)}^T{W}^{1/2}{R}^{-1}{W}^{1/2}\left(y-\tilde{X}{\beta }_1-Qu-Z{v}_1\right)+v_1^T{G}^{-1}{v}_1,$$

and $W^{1/2}$ represents the square root of the matrix $W.$

Estimating the parameters of (3.9) by minimizing ${\Phi }_W$ is equivalent to estimating those given by

$$W^{1/2}y=W^{1/2}\tilde{X}{\beta }_1+W^{1/2}Qu+W^{1/2}Z{v}_1+e_1.(3.10)$$

The weighted EBLUP based on (3.9) with the matrix of weights given by $W$ corresponds to a classical EBLUP obtained from model (3.10). Define $y_w=W^{1/2}y,$ $X_w=\left[W^{1/2}\tilde{X}\text{,}{W}^{1/2}Q\right]$ and $Z_w=W^{1/2}Z.$ Equation (3.10) can be rewritten as

$$y_w=X_w{\beta }_{\text{fixed}}+Z_w{v}_1+e_1.(3.11)$$

Let ${\widehat{\beta }}_{\text{loc,}\text{fixed}}={\left({\widehat{\beta }}_{\text{loc},1}^T,{\widehat{u}}^T\right)}^T$ and ${\widehat{v}}_{\text{loc},1}={\left({\widehat{v}}_{\text{loc},11},\dots ,{\widehat{v}}_{\text{loc},1M}\right)}^T$ be the EBLUP estimators of the fixed and random effects of (3.11). The estimators ${\widehat{\beta }}_{\text{loc,}\text{fixed}}$ and ${\widehat{v}}_{\text{loc},1}$ are based on local estimators of the variance components $\left({\sigma }_{1v}^2,{\sigma }_{1e}^2\right).$ The estimators of these components, denoted as $\left({\widehat{\sigma }}_{\text{loc},1v}^2,{\widehat{\sigma }}_{\text{loc},1e}^2\right),$ are obtained using HFC or REML methods under model (3.11). Given that $u={\left(m_0\left(p_{l|k}\right),m_0^{\left(1\right)}\left(p_{l|k}\right)/1!,\dots ,m_0^{\left(q\right)}\left(p_{l|k}\right)/q!\right)}^T,$ an estimator ${\widehat{m}}_0\left(p_{l|k}\right)$ of $m_0\left(p_{l|k}\right)$ is the first component ${\widehat{u}}_0$ of $\widehat{u}.$

Notice that ${\widehat{\beta }}_{\text{loc},1},$ ${\widehat{m}}_0\left(p_{l|k}\right)$ and ${\widehat{v}}_{\text{loc},1k}$ could be used to obtain local estimates ${\widehat{y}}_{\text{loc},kl}$ for the unknown value $y_{kl},$ where ${\widehat{y}}_{\text{loc},kl}={\tilde{x}}_{kl}^T{\widehat{\beta }}_{\text{loc},1}+{\widehat{m}}_0\left(p_{l|k}\right)+{\widehat{v}}_{\text{loc},1k}$ for $l\in {\overline{s}}_k.$ However, a referee pointed out that, in practice, this methodology would not likely to be well behaved because it requires a strong balance of the small areas across the range of the probabilities $p_{l|k}.$ If this balance is not respected, the resulting estimation would suffer severely from this localization. As a consequence, we opted for a global estimation of ${\beta }_1$ and $v_1.$

We now explain the second step of our procedure. Parameters ${\beta }_1$ and $v_1$ can be estimated globally based on the estimations ${\widehat{m}}_0\left(p_{j|i}\right)$ and the auxiliary data ${\tilde{x}}_{ij}$ associated with the sample units. For $j\in s_i$ and $i=1,\dots ,M,$ define a new variable, say $\xi ,$ as

$${\xi }_{ij}=y_{ij}-{\widehat{m}}_0\left(p_{j|i}\right),j\in s_i;i=1,\dots ,M.$$

The $n$ values ${\xi }_{ij}$ represent the differences between the observed $y_{ij}’\text{s}$ and their local estimators ${\widehat{m}}_0\left(p_{j|i}\right).$ Using model (3.3), $\xi$ satisfies the following model

$${\xi }_{ij}={\tilde{x}}_{ij}^T{\beta }_{\text{glo},1}+v_{\text{glo},1i}+e_{\text{glo},1ij},j\in s_i;i=1,\dots ,M,(3.12)$$

where $v_{\text{glo},1i}\sim N\left(0,{\sigma }_{\text{glo},1v}^2\right)$ and $e_{\text{glo},1ij}\sim N\left(0,{\sigma }_{\text{glo},1e}^2\right).$ The subscript glo indicates that (3.12) is a global model.

Given that (3.12) represents a parametric linear mixed effects model, we can use the classical (unweighted) EBLUP to estimate its parameters. Let ${\widehat{\beta }}_{\text{glo},1}$ and ${\widehat{v}}_{\text{glo},1i}$ be the respective empirical best linear unbiased estimators of ${\beta }_{\text{glo},1}$ and $v_{\text{glo},1i}.$ Let $\left({\widehat{\sigma }}_{\text{glo},1v}^2,{\widehat{\sigma }}_{\text{glo},1e}^2\right)$ be the estimators of the variance components $\left({\sigma }_{\text{glo},1v}^2,{\sigma }_{\text{glo},1e}^2\right)$ where HFC or REML can be used to estimate these parameters. We estimate $\left({\beta }_1,v_{1i},{\sigma }_{1v}^2,{\sigma }_{1e}^2\right)$ of model (3.3) by $\left({\widehat{\beta }}_{\text{glo},1},{\widehat{v}}_{\text{glo},1i},{\widehat{\sigma }}_{\text{glo},1v}^2,{\widehat{\sigma }}_{\text{glo},1e}^2\right)$ using model (3.12). The global estimators ${\widehat{\beta }}_{\text{glo},1},$ ${\widehat{v}}_{\text{glo},1i}$ and $\left({\widehat{\sigma }}_{\text{glo},1v}^2,{\widehat{\sigma }}_{\text{glo},1e}^2\right)$ are free of bias caused by informative sampling design because ${\xi }_{ij}$ is no longer related to the $p_{j|i}’\text{s}$ after conditioning on $x_{ij}.$

The third step estimates the non observed $y_{ij}$ values, for $j\in {\overline{s}}_i$ and $i=1,\dots ,M,$ by plugging into equation (3.4): i. the local estimators ${\widehat{m}}_0\left(p_{j|i}\right)$ for $j\in {\overline{s}}_i,$ obtained in the first step, and ii. the global estimators ${\widehat{\beta }}_{\text{glo},1}$ and ${\widehat{v}}_{\text{glo},1i}$ obtained in the second step. The resulting ${\widehat{y}}_{ij}’\text{s,}$ for $j\in {\overline{s}}_i,$ are inserted into (3.2) to compute the estimator ${\widehat{\overline{Y}}}_i.$ Note that ${\widehat{\overline{Y}}}_i$ requires ${\tilde{x}}_{ij}$ and $p_{j|i}$ are known for all the units of the population. A referee pointed out that, in practice, this assumption may limit the applicability of the proposed procedure. This could be remedied if National Statistical Offices provided access to the selection probabilities of all units, as they may be needed in applications such as this one.

3.2  Bandwidth selection

Local polynomials require the specification of the kernel $K(\cdot ),$ the order of the polynomial fit $q,$ as well as the bandwidth $h.$ Fan and Gijbels (1996) state that values of $q$ larger than 1 do not bring a significant improvement as compared to the linear fit $\left(q=1\right).$ Fan and Gijbels (1996) also state that the choice of $h$ is far more important than the degree of the polynomial. In what follows, we use a normal density kernel, and chose $q$  equal to one, as this leads to satisfactory results for most applications.

The optimal $h$  is determined using the cross-validation method (CV). For a given $h,$  compute the estimator of $y_{ij}$ given by (3.4) using the sample that remains after the $j^{\text{th}}$ unit has been removed from $s_i.$ Denoting the resulting estimator of $y_{ij}$ as ${\tilde{y}}_{ij},$ we follow Wu and Zhang (2002) and define the CV criterion as

$$\text{CV}\left(h\right)=\frac{1}{M}{\displaystyle \sum _{i=1}^M\frac{1}{n_i}{\displaystyle \sum _{j\in s_i}{\left(y_{ij}-{\tilde{y}}_{ij}\right)}^2.}}$$

The term $1/n_i$ takes into account the number of observations within small area $U_i.$ The optimal bandwidth $h_{\text{opt}}$ is obtained by minimizing the $\text{CV}\left(h\right).$ Given $h_{\text{opt}},$ the local polynomial estimator of the small area mean ${\overline{Y}}_i$ given by (3.2) is denoted as ${\widehat{\overline{Y}}}_i^{\text{LP}}.$

• Previous
• Article main page
• Next