Development of a small area estimation system at Statistics Canada

Section 5. Hierarchical Bayes (HB) method

• Previous
• Next

The basic Fay-Herriot area level model includes a linear sampling model for direct survey estimates and a linear linking model for the parameters of interest. Such models are matched because ${\theta }_i$ appears as a linear function in both the sampling and linking models. There are instances when these equations are not matched such as when a function, $h\left({\theta }_i\right),$ is modelled as a linear function of explanatory variables instead of ${\theta }_i.$ The sampling model and linking model pair is

$${\widehat{\theta }}_i={\theta }_i+e_i(5.1)$$

and

$$h\left({\theta }_i\right)=z_i^T\beta +b_i{v}_i(5.2)$$

where $e_i\stackrel{\text{ind}}{{\displaystyle \sim }}N\left(0,{\psi }_i\right)$ and $v_i\stackrel{\text{ind}}{{\displaystyle \sim }}N\left(0,{\sigma }_v^2\right).$

The model pair given by (5.1) and (5.2) is referred to as an unmatched model. Nonlinear linking models are often needed in practice to provide a better model fit to the data. For example, if the parameter of interest is a probability or a rate within the range of 0 and 1, a linear linking model with normal random effects may not be appropriate. A linking model, in this case, could be a logistic or log-linear model. Such a model was used to adjust counts for detailed levels for the 2011 Census of Canada. A good description of what is involved to carry out such an adjustment can be found in Dick (1995) and You, Rao and Dick (2004).

The production system includes the following choices of $h\left({\theta }_i\right)$

$$h\left({\theta }_i\right)=\{\begin{array}{ll}{\theta }_i\hfill & :\text{ Matched Fay-Herriot }\left(\text{FH}\right)\text{ model}\hfill \\ \log \left({\theta }_i\right)\hfill & :\text{ Unmatched log-linear model}\hfill \\ \log \left({\theta }_i/\left({\theta }_i+C_i\right)\right)\hfill & :\text{ Unmatched log census undercount model}\text{.}\hfill \end{array}(5.3)$$

The inclusion of $h\left({\theta }_i\right)={\theta }_i$ corresponds to the matched model represented by equations (3.1) and (3.2). An advantage of choosing the Hierarchical Bayes method is that the estimated ${\sigma }_v^2$ cannot be negative. The function $\log \left({\theta }_i\right)\text{,}$ where ${\theta }_i$ is equal to the population mean ${\overline{Y}}_i,$ was used in Fay and Herriot (1979). Their context was to estimate per capita income (PCI) for small places in the United States with a population less than 1,000. The function $h\left({\theta }_i\right),$ $\log \left({\theta }_i/\left({\theta }_i+C_i\right)\right),$ was included to support the methodology to estimate the net undercoverage in Canadian Censuses. In this model, ${\theta }_i$ represents the number of individuals not counted in the census, while $C_i$ is the known census count. As a result, ${\theta }_i/\left({\theta }_i+C_i\right)$ is the proportion of individuals undercounted by the Census.

The sampling variances, ${\psi }_i,$ are assumed known for all the linking models represented by (5.2). The variances are assumed to be estimated for the first two functions (the matched Fay-Herriot and unmatched log-linear model) given in (5.3). If the sampling variances, ${\psi }_i,$ are assumed known, then the unknown parameters in the sampling model (5.1) and the linking model (5.2) can be presented in a hierarchical Bayes (HB) framework as follows

$$\left[{\widehat{\theta }}_i|{\theta }_i\right]\sim N\left({\theta }_i,{\psi }_i\right),i=1,\dots ,m$$

and

$$\left[h\left({\theta }_i\right)|\beta ,{\sigma }_v^2\right]\sim N\left(z_i^T\beta ,b_i^2{\sigma }_v^2\right).$$

If the sampling variances are unknown, they are estimated by adding

$$\left[d_i{\stackrel{⌢}{\psi }}_i|{\psi }_i\right]\sim {\psi }_i{\chi }_{d_i}^2$$

where ${\chi }_{d_i}^2$ follows a chi-square distribution with $d_i=\left(n_i-1\right)$ degrees of freedom.

The model parameters $\beta ,$ ${\sigma }_v^2$ and ${\psi }_i$ (when it is unknown) are assumed to obey prior distributions. The distributions used in the production system for $\beta$ and ${\sigma }_v^2$ are the flat prior, $\pi \left(\beta \right)\propto 1,$ and $\pi \left({\sigma }_v^2\right)\propto {\left({\sigma }_v^2\right)}^{-1/2}.$ If ${\psi }_i$ is estimated, the prior $\pi \left({\psi }_i\right)\propto {\left({\psi }_i\right)}^{-1/2}$ is added to the Bayesian model. These prior distributions are multiplied by the density functions of the distributions associated with the sampling and linking models. This yields a joint likelihood function in terms of the model parameters. This function is used to obtain a full conditional (posterior) distribution for each of the unknown parameters. For some of these, the resulting distribution has a tractable or well-known form. For others, the resulting distribution is a product of density functions with no known form. All HB methods involve estimation of the model parameters through repeated sampling of their respective full conditional distributions.

Markov Chain Monte Carlo (MCMC) methods are used to obtain estimates from the full conditional distribution of each parameter. Gibbs sampling is used repeatedly to sample from the full conditional distributions. The Gibbs sampling method (Gelfand and Smith, 1990) with the Metropolis-Hastings algorithm (Chib and Greenberg, 1995) are used to find the posterior means and posterior variances; see Estevao et al. (2015) for details. The various estimators of ${\theta }_i$ resulting from (5.3) are denoted as ${\widehat{\theta }}_i^{\text{HB}}.$

5.1  Benchmarked HB estimator

Benchmarking of the estimators uses the difference adjustment method described in Section 3.2. That is, the benchmarked estimators ${\widehat{\theta }}_i^{\text{HB}}$ are computed as

$${\widehat{\theta }}_i^{\text{HB},b}=\{\begin{array}{ll}{\widehat{\theta }}_i^{\text{HB}}+{\alpha }_i\left({\widehat{\theta }}^{*}-{\displaystyle {\sum }_{d\in A}{\omega }_d{\widehat{\theta }}_d^{\text{HB}}}\right)\hfill & \text{for}i\in A\hfill \\ z_i^T\widehat{\beta }\hfill & \text{for}i\in \overline{A}\hfill \end{array}$$

where ${\alpha }_i={\left({\displaystyle {\sum }_{i\in A}{\omega }_i^2\left({\ddot{\psi }}_i+b_i^2{\widehat{\sigma }}_v^{2\text{HB}}\right)}\right)}^{-1}{\omega }_i\left({\ddot{\psi }}_i+b_i^2{\widehat{\sigma }}_v^{2\text{HB}}\right)$ for $i\in A,$ and ${\widehat{\theta }}^{*}$ is the benchmark value. The terms ${\omega }_i$ are defined as follows: ${\omega }_i=1$ if the benchmarking is to a total, and ${\omega }_i=N_i/N$ if the benchmarking is for the mean. The ${\ddot{\psi }}_i’s$ are either known or unknown. The ${\widehat{\theta }}^{*}$ can be a value provided by the user that represents the total or mean of the $y$ -values of population $U.$ The benchmarking ensures that ${\displaystyle {\sum }_{i\in A\cup \overline{A}}{\omega }_i{\widehat{\theta }}_i^{\text{HB},b}}={\widehat{\theta }}^{*}.$

• Previous
• Article main page
• Next