`library(BayesSampling)`

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.:

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:

where *X* is a covariate matrix of dimension *N* × *p*, with rows *X*_{i} = (*x*_{i1},...,*x*_{ip}), *i* = 1, ..., *N*;

*β* = (*β*_{1},...,*β*_{p})′ is a *p* × 1 vector of unknown parameters; and *y*, given *β*, is a random vector with mean *X**β* and known covariance matrix *V* of dimension *N* × *N*. Analogously *a* and *R* are the respective *p* × 1 prior mean vector and *p* × *p* prior covariance matrix of *β*.

Since the response vector *y* is partitioned into *y*_{s} and *y*_{s̄}, the matrix *X*, which is assumed to be known, is analogously partitioned into *X*_{s} and *X*_{s̄}, and *V* is partitioned into *V*_{s}, *V*_{s̄}, *V*_{ss̄} and *V*_{s̄s}. The first aim is to predict *y*_{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:

(…)

The *BLE_Reg()* function will apply this methodology to the given sample, calculate the Bayes Linear Estimator (and its associate variance) to the parameter *β* and for the individuals not in the sample, given the auxiliary variable values. In a simple model the auxiliary variable will have value 1 for every individual.

- Example from the help page

```
<- matrix(c(1,1,1,1,2,3,5,0),nrow=4,ncol=2)
xs <- c(12,17,28,2)
ys <- matrix(c(1,1,1,0,1,4),nrow=3,ncol=2)
x_nots <- c(1.5,6)
a <- matrix(c(10,2,2,10),nrow=2,ncol=2)
R <- diag(c(1,1,1,1))
Vs <- diag(c(1,1,1))
V_nots
<- BLE_Reg(ys,xs,a,R,Vs,x_nots,V_nots)
Estimator
Estimator#> $est.beta
#> Beta
#> 1 1.723363
#> 2 5.206675
#>
#> $Vest.beta
#> V1 V2
#> 1 0.6708234 -0.17568311
#> 2 -0.1756831 0.07225381
#>
#> $est.mean
#> y_nots
#> 1 1.723363
#> 2 6.930039
#> 3 22.550064
#>
#> $Vest.mean
#> V1 V2 V3
#> 1 1.67082340 0.49514029 -0.03190904
#> 2 0.49514029 1.39171098 0.08142307
#> 3 -0.03190904 0.08142307 1.42141940
#>
#> $est.tot
#> [1] 90.20347
#>
#> $Vest.tot
#> [1] 5.573262
```