2 Bayes linear estimation for finite population

Kelly Cristina M. Gonçalves, Fernando A. S. Moura et Helio S. Migon

Previous | Next

The Bayes approach has been found to be successful in many applications, particularly when the data analysis has been improved by expert judgements. But while Bayesian models have many appealing features, their application often involves the full specification of a prior distribution for a large number of parameters. Goldstein and Wooff (2007), section 1.2, argue that as the complexity of the problem increases, our actual ability to fully specify the prior and/or the sampling model in detail is impaired. They conclude that in such situations, there is a need to develop methods based on partial belief specification.

Hartigan(1969) proposed an estimation method, termed Bayes linear estimation approach, that only requires the specification of first and second moments. The resulting estimators have the property of minimizing posterior squared error loss among all estimators that are linear in the data and can be thought of as approximations to posterior means. The Bayes linear estimation approach is fully employed in this article and is briefly described below.

2.1 Bayes linear approach

Let $y_s$ be the vector with observations and $\theta$ be the parameter to be estimated. For each value of $\theta$ and each possible estimate $d,$ belonging to the parametric space $\Theta ,$ we associate a quadratic loss function $L\left(\theta \mathrm{,}d\right)=\left(\theta -d\right)\prime \left(\theta -d\right)=tr\left(\theta -d\right)\left(\theta -d\right)\prime .$ The main interest is to find the value of $d$ that minimizes $r\left(d\right)=E\left[L\left(\theta \mathrm{,}d\right)|y_s\right],$ the conditional expected value of the quadratic loss function given the data.

Suppose that the joint distribution of $\theta$ and $y_s$ is partially specified by only their first two moments:

$$\left(\begin{array}{c}\theta \\ y_s\end{array}\right)\sim \left[\left(\begin{array}{c}a\\ f\end{array}\right)\mathrm{,}\left(\begin{array}{cc}R& AQ\\ Q{A}^{\prime }& Q\end{array}\right)\right]\mathrm{,}\text{ (2}\text{.1)}$$

where $a$ and $f,$ respectively, denote mean vectors and $R,$ $AQ$ and $Q$ the covariance matrix elements of $\theta$ and $y_s\mathrm{.}$

The Bayes linear estimator (BLE) of $\theta$ is the value of $d$ that minimizes the expected value of this quadratic loss function within the class of all linear estimates of the form $d=d\left(y_s\right)=h+H{y}_s,$ for some vector $h$ and matrix $H.$ Thus, the BLE of $\theta ,$ $\widehat{d},$ and its associated variance, $\widehat{V}\left(\widehat{d}\right),$ are respectively given by:

$$\widehat{d}=a+A\left(y_s-f\right)\text{and}\widehat{V}\left(\widehat{d}\right)=R-AQA\prime \mathrm{.}\text{ (2}\text{.2)}$$

It should be noted that the BLE depends on the specification of the first and second moments of the joint distribution partially specified in (2.1). The issue of eliciting these quantities is dealt with in sections (2.3.1) and (4.1) for some particular cases.

2.2 Bayes linear approach to finite population

Consider $U=\left\{u_1\mathrm{,}\dots \mathrm{,}{u}_N\right\}$ a finite population with $N$ units. Let $y=\left(y_1\mathrm{,}\dots \mathrm{,}{y}_N\right)\prime$ be the vector with the values of interest of the units in $U.$ The response vector $y$ is partitioned into the known observed $n­$ sample vector $y_s,$ and the non-observed vector $y_{\overline{s}}$ of dimension $N-n.$ The general problem is to predict a function of the vector $y,$ such as the total $T={\displaystyle {\sum }_{i=1}^N}{y}_i={1^{\prime }}_s{y}_s+{1^{\prime }}_{\overline{s}}{y}_{\overline{s}},$ where $1_s$ and $1_{\overline{s}}$ are the vectors of 1's of dimensions $n$ and $N-n,$ respectively. In the model-based approach, this is usually done by assuming a parametric model for the population values $y_i’s$ and then obtaining the Empirical Best Linear Unbiased Predictor (EBLUP) for the unknown vector $y_{\overline{s}}$ under this model. Usually, the mean square error of the EBLUP of $T$ is obtained by second order approximation, as well as an unbiased estimator of it. See Valliant, Dorfman and Royall (2000), chapter 2, for details.

The Bayesian approach to finite population prediction often assumes a parametric model, but it aims to find the posterior distribution of $T$ given $y_s.$ Point estimates can be obtained by setting a loss function, although in many practical problems, the posterior mean is often considered and its associated variance is given by the posterior variance, i.e.:

$$E\left(T|y_s\right)={1^{\prime }}_s{y}_s+{1^{\prime }}_{\overline{s}}E\left(y_{\overline{s}}|y_s\right)\text{and}V\left(T|y_s\right)={1^{\prime }}_{\overline{s}}V\left(y_{\overline{s}}|y_s\right)1_{\overline{s}}\mathrm{.}\text{ (2}\text{.3)}$$

It is possible to obtain an approximation to the quantities in (2.3) by using a Bayes linear estimation approach. Here, we will particularly obtain the estimators by assuming a general two-stage hierarchical model for finite population, specified only by its mean and variance-covariance matrix, presented in Bolfarine and Zacks (1992), page 76. Particular cases describing usual population structures found in practice are easily derived from (2.4). The general model can be written as:

$$y|\beta \sim \left[X\beta \text{,}V\right]\text{and}\beta \sim \left[a\mathrm{,}R\right]\mathrm{,}\text{ (2}\text{.4)}$$

where $X$ is a covariate matrix of dimension $N\times p,$ with rows $X_i=\left(x_{i1}\mathrm{,}\dots \mathrm{,}{x}_{ip}\right),$ $i=\mathrm{1,}\dots \mathrm{,}N;$ $\beta ={\left({\beta }_1\mathrm{,}\dots \mathrm{,}{\beta }_p\right)}^{\prime }$ is a $p\times 1$ vector of unknown parameters; and $y,$ given $\beta ,$ is a random vector with mean $X\beta$ and known covariance matrix $V$ of dimension $N\times N.$ Analogously $a$ and $R$ are the respective $p\times 1$ prior mean vector and $p\times p$ prior covariance matrix of $\beta .$

Since the response vector $y$ is partitioned into $y_s,$ and $y_{\overline{s}},$ the matrix $X,$ which is assumed to be known, is analogously partitioned into $X_s$ and $X_{\overline{s}},$ and $V$ is partitioned into $V_s,$ $V_{\overline{s}},$ $V_{s\overline{s}}$ and $V_{\overline{s}s}.$ The first aim is to predict $y_{\overline{s}}$ given the observed sample $y_s$ and then the total $T.$ We did this in the following steps: first, we used a joint prior distribution that is only partially specified in terms of moments, as follows:

$$\left(\begin{array}{c}{y}_{\overline{s}}\\ y_s\end{array}\right)|\beta \sim \left[\left(\begin{array}{c}{X}_{\overline{s}}\beta \\ X_s\beta \end{array}\right)\mathrm{,}\left(\begin{array}{cc}{V}_{\overline{s}}& V_{\overline{s}s}\\ V_{s\overline{s}}& V_s\end{array}\right)\right]\mathrm{.}$$

Therefore, applying the general result in equation (2.2), the BLE of $E\left(y_{\overline{s}}|y_s\mathrm{,}\beta \right)$ and the minimum expected square loss (associated variance) are given by:

$$\widehat{E}\left(y_{\overline{s}}|y_s\mathrm{,}\beta \right)=X_{\overline{s}}\beta +V_{\overline{s}s}{V}_s^{-1}\left(y_s-X_s\beta \right)\text{and}\widehat{V}\left(y_{\overline{s}}|y_s\mathrm{,}\beta \right)=V_{\overline{s}}-V_{\overline{s}s}{V}_s^{-1}{V}_{s\overline{s}}\mathrm{.}\text{ (2}\text{.5)}$$

Remark 1: It should be noted that if normality is assumed then $E\left(y_{\overline{s}}|y_s\mathrm{,}\beta \right)$ and $V\left(y_{\overline{s}}|y_s\mathrm{,}\beta \right)$ are respectively given by the right sides of (2.5). The BLE in (2.5) and its associated variance can be viewed respectively as approximations of $E\left(y_{\overline{s}}|y_s\mathrm{,}\beta \right)$ and $V\left(y_{\overline{s}}|y_s\mathrm{,}\beta \right)$ for non-normality cases.

Now, if we come back to model (2.4), we need to adapt the structure (2.1) and use the results in (2.2) to obtain the BLE of $\beta$ and its associate variance, $\widehat{V}\left(\widehat{\beta }\right),$ respectively given by:

$$\widehat{\beta }=a+R{X^{\prime }}_s{\left(X_{s}R{X^{\prime }}_s+V_s\right)}^{-1}\left(y_s-X_{s}a\right)\text{and}\widehat{V}\left(\widehat{\beta }\right)=C=R-R{X^{\prime }}_s{\left(X_{s}R{X^{\prime }}_s+V_s\right)}^{-1}{X}_{s}R\mathrm{.}\text{ (2}\text{.6)}$$

It is easy to show that the first equation in (2.6) can be rewritten as $\widehat{\beta }=C\left({X^{\prime }}_s{V}_s^{-1}{y}_s+R^{-1}a\right),$ where $C^{-1}=R^{-1}+{X^{\prime }}_s{V}_s^{-1}{X}_s.$ It should be noted that if we place a vague prior distribution on $\beta ,$ taking $R^{-1}\to 0,$ we obtain the minimum least square estimator of $\beta :$ ${\widehat{\beta }}_{LS}={\left({X^{\prime }}_s{V}_s^{-1}{X}_s\right)}^{-1}{X^{\prime }}_s{V}_s^{-1}{y}_s\mathrm{.}$

Now, applying well known properties of conditional expectations and variances, we obtain:

$$E\left[y_{\overline{s}}|y_s\right]=E\left(E\left(y_{\overline{s}}|y_s\mathrm{,}\beta \right)|y_s\right)\text{and}V\left[y_{\overline{s}}|y_s\right]=E\left(V\left(y_{\overline{s}}|y_s\mathrm{,}\beta \right)|y_s\right)+V\left(E\left(y_{\overline{s}}|y_s\mathrm{,}\beta \right)|y_s\right).\text{ (2}\text{.7)}$$

Replacing $E\left(y_{\overline{s}}|y_s\mathrm{,}\beta \right)$ and $V\left(y_{\overline{s}}|y_s\mathrm{,}\beta \right)$ in (2.7) with their respective BLE's in (2.5) and in turn, replacing $E\left(\beta |y_s\right)$ and $V\left(\beta |y_s\right)$ with $\widehat{\beta }$ and $\widehat{V}\left(\widehat{\beta }\right)$ in (2.6), we obtain the BLE of $E\left[y_{\overline{s}}|y_s\right]$ and its associated variance as:

$$\begin{array}{c}\widehat{E}\left[y_{\overline{s}}|y_s\right]=X_{\overline{s}}\widehat{\beta }+V_{\overline{s}s}{V}_s^{-1}\left(y_s-X_s\widehat{\beta }\right)\text{and}\\ \widehat{V}\left[y_{\overline{s}}|y_s\right]=V_{\overline{s}}-V_{\overline{s}s}{V}_s^{-1}{V}_{s\overline{s}}+\left(X_{\overline{s}}-V_{\overline{s}s}{V}_s^{-1}{X}_s\right)C\left(X_{\overline{s}}-V_{\overline{s}s}{V}_s^{-1}{X}_s\right)\prime \mathrm{.}\text{ (2}\text{.8)}\end{array}$$

Remark 2: Analogously to the Remark 1, when normality is assumed we have that the right sides of (2.8) are respectively the values of $E\left[y_{\overline{s}}|y_s\right]$ and $V\left[y_{\overline{s}}|y_s\right].$

The general expression of BLE for the total $T$ and its associated variance are respectively obtained by replacing $E\left(y_{\overline{s}}|y_s\right)$ and $V\left(y_{\overline{s}}|y_s\right)$ in equations in (2.3) with their respective counterparts $\widehat{E}\left[y_{\overline{s}}|y_s\right]$ and $\widehat{V}\left[y_{\overline{s}}|y_s\right]:$

$$\widehat{T}={1^{\prime }}_s{y}_s+{1^{\prime }}_{\overline{s}}\widehat{E}\left[y_{\overline{s}}|y_s\right]\text{and}\widehat{V}\left(\widehat{T}\right)={1^{\prime }}_{\overline{s}}\widehat{V}\left[y_{\overline{s}}|y_s\right]1_{\overline{s}}\mathrm{.}\text{ (2}\text{.9)}$$

It should be noted that in many applications of (2.9), the matrix $V$ is assumed diagonal, which implies $V_{\overline{s}s}=0$ and then we have:

$$\widehat{T}={1^{\prime }}_s{y}_s+{1^{\prime }}_{\overline{s}}{X}_{\overline{s}}\widehat{\beta }\text{and}\widehat{V}\left(\widehat{T}\right)={1^{\prime }}_{\overline{s}}\left[V_{\overline{s}}+X_{\overline{s}}C{X^{\prime }}_{\overline{s}}\right]1_{\overline{s}}\mathrm{.}\text{ (2}\text{.10)}$$

For the sake of illustration, we consider some examples discussed by O'Hagan (1985) and propose a new ratio estimator, which is one of the contributions of our work. All of them can be treated as special cases of the model (2.4).

2.3 Revisiting some common survey designs

2.3.1 Simple random sampling without replacement: Second order exchangeability

O'Hagan (1985) considered the simple case where the population has no relevant structure, which can be done by setting up:

$$E\left(y_i\right)=m\mathrm{,}V\left(y_i\right)=v\text{and}\text{Cov}\left(y_i\mathrm{,}{y}_j\right)=c\mathrm{,}i\mathrm{,}j=\mathrm{1,}\dots \mathrm{,}N\mathrm{,}\forall i\ne j\mathrm{.}\text{ (2}\text{.11)}$$

Remark 3: The correlation introduced in model (2.11) can be justified to mimic simple random sampling without replacement.

Applying the general result established in (2.10) to (2.11) with $\beta$ of dimension 1, $X=1_N,$ $a=m,$ $R=c$ and $V={\sigma }^{2}I,$ where ${\sigma }^2=v-c,$ we obtain the BLE of $T$ and its respective associated variance:

$${\widehat{T}}_{srs}=n{\overline{y}}_s+\left(N-n\right)\widehat{\mu }\text{and}\widehat{V}\left({\widehat{T}}_{srs}\right)=\left(N-n\right){\sigma }^2+{\left(N-n\right)}^{2}c{\sigma }^2{\left({\sigma }^2+nc\right)}^{-1}\mathrm{,}\text{ (2}\text{.12)}$$

where

$$\begin{array}{c}{\overline{y}}_s=n^{-1}{1^{\prime }}_s{y}_s\text{is}\text{the}\text{sample}\text{mean}\mathrm{,}\\ \widehat{\mu }=\omega {\overline{y}}_s+\left(1-\omega \right)m\text{is}\text{the}\text{expected}\text{value}\text{of}\text{the}\text{non-observed}\text{values}\text{of}y\text{and}\\ \omega =\frac{n{\sigma }^{-2}}{c^{-1}+n{\sigma }^{-2}}\mathrm{,}\text{where}{\sigma }^2=v-c\mathrm{.}\end{array}$$

It should be noted that $\widehat{\mu }$ is a weighted average of the prior mean $m$ and the sample mean ${\overline{y}}_s,$ where $\omega$ is the ratio between two population quantities. The mean $m$ can be viewed as the investigator's prior of the true population mean $\overline{y}.$ The uncertainty about $y_i$ is split into two components: the uncertainty about the overall level of the $y_i’\text{s}$ (between variation) and the one with respect to how much each $y_i$ may vary from that overall level (within variation). A useful measure of variability of units within the population is given by

$$S^2=\frac{1}{N-1}{\displaystyle \sum _{i=1}^N{\left(y_i-\overline{y}\right)}^2}.$$

It is not difficult to show that $E\left(S^2\right)=v-c={\sigma }^2.$ Therefore, ${\sigma }^2$ can be interpreted as a prior estimate of variability within the population. We also obtain $V\left(\overline{y}\right)=c+N^{-1}{\sigma }^2.$ In many applications, $N$ is large and thus the constant $c$ can be viewed as the between variation.

Letting $v\to \infty$ and keeping ${\sigma }^2$ fixed, that is, assuming prior ignorance, the estimates in (2.12) yield:

$$\begin{array}{c}{\widehat{T}}_{srs}=N{\overline{y}}_s\text{and}\widehat{V}\left({\widehat{T}}_{srs}\right)=N^2\left(1-\frac{n}{N}\right)\frac{{\sigma }^2}{n}\mathrm{.}\end{array}$$

These expressions are very similar to the well-known total estimate and its variance in the design-based context for the simple random sampling case. O'Hagan (1985) discussed some possibilities to avoid the difficult task of assigning a value for ${\sigma }^2.$ The most natural way to do this is to find the BLE of it, but linear in the squares and cross product variance terms. However, it requires to specify fourth order moments of the $y_i’\text{s}\text{.}$ Goldstein (1979) proposed a BLE for the variance, which uses only linear functions of data. Nevertheless, it results in a complicated expression to its associated variance of his modified BLE. O'Hagan (1985) argued that if prior information about variance components is weak, any posterior estimate is close to the standard non-Bayesian estimates using only the data, wherever such estimate is available. Therefore, he suggested, as an approximate Bayesian procedure, substituting these standard variance estimates into the BLE and its associated variance wherever appropriate. For this case, we can replace ${\sigma }^2$ with $s^2={\left(n-1\right)}^{-1}{\displaystyle {\sum }_{i=1}^n}{\left(y_i-{\overline{y}}_s\right)}^2,$ which is design-based unbiased for $S^2.$

2.3.2  Stratified simple random sampling without replacement

Denote by $y_{hi}$ the $i^{\text{th}}$ unit, $i=\mathrm{1,...,}{N}_h$ belonging to stratum $h,$ $h=\mathrm{1,..,}H.$ It is assumed that the stratum sizes, $N_h,$ are known for all strata. The second-order exchangeability within each stratum is stated in O'Hagan (1985) as:

$$E\left(y_{hi}\right)=m_h\mathrm{,}V\left(y_{hi}\right)=v_h\mathrm{,}\text{Cov}\left(y_{hi}\mathrm{,}{y}_{hj}\right)=c_h\mathrm{,}i\ne j\text{and}\text{Cov}\left(y_{hi}\mathrm{,}{y}_{lj}\right)=d_{hl}\mathrm{,}h\ne l\mathrm{.}$$

Remark 4: It is reasonable to assume that the information gained about one stratum could change the beliefs about other strata in some special applications. However, if we want to mimic the stratified simple random sampling, we should assume that observations in different strata are uncorrelated, letting $d_{hl}=0.$

The general model (2.4) can be applied to this case by setting $X=\text{diag}\left(X_1\mathrm{,}\dots \mathrm{,}{X}_H\right)$ and $V=\text{diag}\left(V_1\mathrm{,}\dots \mathrm{,}{V}_H\right)$ , with $X_h=1_{N_h}$ and $V_h={\sigma }_h^2{I}_{N_h},$ where ${\sigma }_h^2=v_h-c_h,$ $\forall$ $h=\mathrm{1,}\dots \mathrm{,}H\mathrm{,}$ $a=\left(m_1\mathrm{,}\dots \mathrm{,}{m}_H\right)\prime ,$ $R$ is an $H\times H$ matrix with $R_{hl}=c_h,$ if $h=l$ and $R_{hl}=d_{hl}$ otherwise. The BLE of $T$ and its associated variance are obtained from (2.10) and can be found in O'Hagan (1985). Models for cluster sampling can be found in Bolfarine and Zacks (1992), page 11. The BLE for cluster models can be seen in O'Hagan (1985).

Previous | Next