Development of a small area estimation system at Statistics Canada

Section 3. Area level model

• Previous
• Next

The area level small area estimator first appeared in the seminal paper of Fay and Herriot (1979). Following that paper, let the parameter of interest be ${\theta }_i;$ common examples are totals, $Y_i={\displaystyle {\sum }_{j\in U_i}{y}_j},$ or means, ${\overline{Y}}_i=Y_i/N_i.$ As noted above, the vector of auxiliary variables may differ from the one used in direct estimation and is denoted as $z.$ The area level model can be expressed as two equations.

The first equation, commonly known as the sampling model, is given by

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

and expresses the direct estimate ${\widehat{\theta }}_i$ in terms of the unknown parameter ${\theta }_i$ plus a random error $e_i$ due to sampling. The sampling errors $e_i$ are independently and identically distributed with mean 0 and variance ${\psi }_i:$ that is $E_p\left(e_i|{\theta }_i\right)=0$ and $V_p\left(e_i|{\theta }_i\right)={\psi }_i,$ where $p$ denotes expectation in terms of the sample design. Note that ${\psi }_i$ is also the design variance of ${\widehat{\theta }}_i$ and is typically unknown.

The second equation, known as the linking model, is given by

$${\theta }_i=z_i^T\beta +b_i{v}_i(3.2)$$

and expresses the parameter ${\theta }_i$ as a fixed effect $z_i^T\beta$ plus a random effect $v_i$ multiplied by $b_i.$ In the production system, the $b_i$ term has a default value of one but can be specified by the user to control heteroscedastic errors or the impact of influential observations. The random effects $v_i$ are independently and identically distributed with mean 0 and unknown model variance ${\sigma }_v^2,$ that is $E_m\left(v_i\right)=0$ and $V_m\left(v_i\right)={\sigma }_v^2$ where $E_m$ denotes the model expectation and $V_m$ the model variance. The random errors $e_i$ are independent of the random effects $v_i.$ The combination of the sampling model and linking model results in a single generalized linear mixed model (GLMM) given by

$${\widehat{\theta }}_i=z_i^T\beta +b_i{v}_i+e_i.(3.3)$$

From the Fay-Herriot model (3.3), we observe that $E_{mp}\left({\widehat{\theta }}_i\right)=z_i^T\beta$ and $V_{mp}\left({\widehat{\theta }}_i\right)=b_i^2{\sigma }_v^2+{\tilde{\psi }}_i,$ where ${\tilde{\psi }}_i=E_m\left({\psi }_i\right)$ is the smoothed design variance of ${\widehat{\theta }}_i.$ In general, we cannot treat ${\psi }_i$ as fixed, as it is not strictly a function of auxiliary data. If the ${\sigma }_v^2’s$ and ${\tilde{\psi }}_i’s$ are known, the solution to the GLMM yields the Best Linear Unbiased Predictor (BLUP), ${\tilde{\theta }}_i^{\text{BLUP}}$

$${\tilde{\theta }}_i^{\text{BLUP}}=\{\begin{array}{ll}{\gamma }_i{\widehat{\theta }}_i+\left(1-{\gamma }_i\right)z_i^T\tilde{\beta }\hfill & \text{for}i\in A\hfill \\ z_i^T\tilde{\beta }\hfill & \text{for}i\in \overline{A}\hfill \end{array}(3.4)$$

where ${\gamma }_i=\left(b_i^2{\sigma }_v^2\right)/\left({\tilde{\psi }}_i+b_i^2{\sigma }_v^2\right)$ and $\tilde{\beta }={\left({\displaystyle {\sum }_{i\in A}{z}_i{z}_i^T}/\left({\tilde{\psi }}_i+b_i^2{\sigma }_v^2\right)\right)}^{-1}{\displaystyle {\sum }_{i\in A}{z}_i{\widehat{\theta }}_i}/\left({\tilde{\psi }}_i+b_i^2{\sigma }_v^2\right).$

There are four recursive procedures for estimating ${\sigma }_v^2$ and $\beta$ in the production system. The first three assume that ${\tilde{\psi }}_i$ is known, or that a smoothed version of it is available (see the following section for details). Under this assumption, the variance components can be computed via the Fay-Herriot procedure (FH) as outlined in Fay and Herriot (1979), the restricted maximum likelihood (REML), or the Adjusted Density Maximization (ADM) due to Li and Lahiri (2010). The fourth procedure, WF, due to Wang and Fuller (2003) assumes that ${\psi }_i$ is estimated by ${\widehat{\psi }}_i$ given that $n_i\ge 2.$ The WF procedure does not require any smoothing of the estimated ${\widehat{\psi }}_i$ values before estimating ${\sigma }_v^2.$ Wang and Fuller (2003) carried out simulations with $n_i$ ranging from 9 to 36 and found that their procedure yielded reasonable estimates of ${\theta }_i$ and its estimated mean squared error.

The main difference between these four procedures is how the ${\sigma }_v^2’s$ are computed. They are all based on an iterative scoring algorithm that obtains ${\widehat{\sigma }}_v^2$ as an estimate of the model variance ${\sigma }_v^2.$ The FH, REML, and WF procedures may yield ${\widehat{\sigma }}_v^2’s$ that are smaller than zero. If this occurs, the ${\widehat{\sigma }}_v^2’s$ are set to zero for both the FH and REML procedures. A drawback of truncating the estimated ${\sigma }_v^2$ to zero is that the resulting small area estimator will be synthetic for all areas. Li and Lahiri (2010) suggested the ADM as a way to address the problem of obtaining negative ${\widehat{\sigma }}_v^2$ by maximizing an adjusted likelihood defined as a product of the model variance and a standard likelihood. Although the ADM method always gives a positive solution for ${\sigma }_v^2,$ it should be used cautiously because it overestimates the model variance. The REML, FH and ADM procedures use the smoothed values of the estimated ${\widehat{\psi }}_i$ values obtained from the sample or some estimate provided by the user. For the WF procedure, if ${\widehat{\sigma }}_v^2<0,$ Wang and Fuller (2003) suggested to set ${\widehat{\sigma }}_v^2$ to $0.5\sqrt{\widehat{V}\left({\widehat{\sigma }}_v^2\right)},$ where

$$\widehat{V}\left({\widehat{\sigma }}_v^2\right)={\displaystyle \sum _{i\in A}2{\kappa }_i^2\left[{\left({\widehat{\psi }}_i+b_i^2{\widehat{\sigma }}_v^2\right)}^2+\frac{{\left({\widehat{\psi }}_i\right)}^2}{\left(n_i-1\right)}\right]}$$

and

$${\kappa }_i=\frac{{\left[b_i^2{\widehat{\sigma }}_v^2+\frac{\left(n_i+1\right)}{\left(n_i-1\right)}{\widehat{\psi }}_i\right]}^{-1}}{{\displaystyle {\sum }_{i\in A}{\left[b_i^2{\widehat{\sigma }}_v^2+\frac{\left(n_i+1\right)}{\left(n_i-1\right)}{\widehat{\psi }}_i\right]}^{-1}}}.$$

Plugging ${\widehat{\sigma }}_v^2$ and an estimate of ${\tilde{\psi }}_i’s$ into the ${\tilde{\theta }}_i^{\text{BLUP}},$ defined by equation (3.4), yields the Empirical Best Linear Unbiased Predictor (EBLUP), ${\widehat{\theta }}_i^{\text{EBLUP}}.$ It is given by

$${\widehat{\theta }}_i^{\text{EBLUP}}=\{\begin{array}{ll}{\widehat{\gamma }}_i{\widehat{\theta }}_i+\left(1-{\widehat{\gamma }}_i\right)z_i^T\widehat{\beta }\hfill & \text{for }i\in A\hfill \\ z_i^T\widehat{\beta }\hfill & \text{for }i\in \overline{A}\hfill \end{array}$$

where ${\widehat{\gamma }}_i=\left(b_i^2{\widehat{\sigma }}_v^2\right)/\left({\ddot{\psi }}_i+b_i^2{\widehat{\sigma }}_v^2\right),$ $\widehat{\beta }={\left({\displaystyle {\sum }_{i\in A}{z}_i{z}_i^T}/\left({\ddot{\psi }}_i+b_i^2{\widehat{\sigma }}_v^2\right)\right)}^{-1}{\displaystyle {\sum }_{i\in A}{z}_i{\widehat{\theta }}_i^{\text{DIR}}}/\left({\ddot{\psi }}_i+b_i^2{\widehat{\sigma }}_v^2\right),$ and ${\ddot{\psi }}_i$ is chosen according to the procedure used. For the REML, FH and ADM procedures the ${\ddot{\psi }}_i’s$ are the smoothed values of the estimated ${\widehat{\psi }}_i$ values obtained from the sample or some estimate provided by the user. For the WF procedure, we have that ${\ddot{\psi }}_i={\widehat{\psi }}_i.$ If the estimated model variance $b_i^2{\widehat{\sigma }}_v^2$ is relatively small compared with ${\ddot{\psi }}_i,$ then ${\widehat{\gamma }}_i$ will be small and more weight will be attached to the synthetic estimator $z_i^T\widehat{\beta }.$ Similarly, more weight is attached to the direct estimator, ${\widehat{\theta }}_i,$ if the design variance ${\ddot{\psi }}_i$ is relatively small.

Details of the required computations can be found in the methodology specifications for the production system in Estevao et al. (2015).

3.1  Estimation of the smooth design variance

The design variance, ${\psi }_i,$ could be used as an estimator of the smooth design variance ${\tilde{\psi }}_i=E_m\left({\psi }_i\right)$ if it were known. In most cases, it is unknown. To get around this difficulty, a design-unbiased variance estimator ${\widehat{\psi }}_i$ of ${\psi }_i$ is assumed to be available; i.e., $E_p\left({\widehat{\psi }}_i\right)={\psi }_i.$ Under this assumption, we have that

$$E_{mp}\left({\widehat{\psi }}_i\right)=E_m\left({\psi }_i\right)={\tilde{\psi }}_i.$$

A simple unbiased estimator of the smooth design variance ${\tilde{\psi }}_i$ is ${\widehat{\psi }}_i.$ However, ${\widehat{\psi }}_i$ may be quite unstable when the sample size in domain $i$ is small. A more efficient estimator is obtained by modelling ${\widehat{\psi }}_i$ given $z_i.$ Dick (1995) and Rivest and Belmonte (2000) considered smoothing models given by

$$\log \left({\widehat{\psi }}_i\right)=x_i^T\alpha +{\epsilon }_i,$$

where $x_i$ is a vector of explanatory variables that are functions of $z_i,$ $\alpha$ is a vector of unknown model parameters to be estimated, and ${\epsilon }_i$ is a random error with $E_{mp}\left({\epsilon }_i\right)=0$ and constant variance ${\sigma }_{\epsilon }^2=V_{mp}\left({\epsilon }_i\right).$ We also assume that the errors ${\epsilon }_i$ are identically distributed conditionally on $z_i,$ $i=1,\dots ,m.$ From the above model, we observe that

$${\tilde{\psi }}_i=E_{mp}\left({\widehat{\psi }}_i\right)=\exp \left(x_i^T\alpha \right)\Delta ,$$

where $\Delta =E_{mp}\left(\exp \left({\epsilon }_i\right)\right).$ Dick (1995) estimated ${\tilde{\psi }}_i$ by omitting the factor $\Delta .$ Rivest and Belmonte (2000) estimated $\Delta$ by assuming that the errors ${\epsilon }_i$ are normally distributed. However, we observed empirically that the resulting estimator of $\Delta$ is sensitive to deviations from the normality assumption. This assumption is avoided by using a method of moments (see Beaumont and Bocci, 2016). This leads to the unbiased estimator of $\Delta$ given by

$$\widehat{\Delta }\left(\alpha \right)=\frac{{\displaystyle {\sum }_{i=1}^m{\widehat{\psi }}_i}}{{\displaystyle {\sum }_{i=1}^m\exp \left(x_i^T\alpha \right)}}.$$

An estimator $\widehat{\alpha }$ of the vector of unknown model parameters $\alpha$ is necessary to estimate ${\tilde{\psi }}_i.$ It is obtained using the ordinary least squares method as

$$\widehat{\alpha }={\left({\displaystyle \sum _{i=1}^m{x}_i{x}_i^T}\right)}^{-1}{\displaystyle \sum _{i=1}^m{x}_i\log \left({\widehat{\psi }}_i\right)}.$$

The estimator ${\widehat{\tilde{\psi }}}_i$ of ${\tilde{\psi }}_i$ is then given by

$${\widehat{\tilde{\psi }}}_i=\exp \left(x_i^T\widehat{\alpha }\right)\widehat{\Delta }\left(\widehat{\alpha }\right).$$

A nice property of ${\widehat{\tilde{\psi }}}_i$ is that the average of the smooth design variance estimator, ${\widehat{\tilde{\psi }}}_i,$ is equal to the average of the direct variance estimator, ${\widehat{\psi }}_i;$ i.e.,

$$\frac{{\displaystyle {\sum }_{i=1}^m{\widehat{\tilde{\psi }}}_i}}{m}=\frac{{\displaystyle {\sum }_{i=1}^m{\widehat{\psi }}_i}}{m}.$$

This ensures that ${\widehat{\tilde{\psi }}}_i$ does not systematically overestimate or underestimate ${\tilde{\psi }}_i=E_{mp}\left({\widehat{\psi }}_i\right).$

3.2  Benchmarking

If the parameter of interest ${\theta }_i$ is a total $\left({\theta }_i=Y_i\right),$ the user may wish to have the sum of the small area estimates, $\widehat{\theta }={\displaystyle {\sum }_{i\in A\cup \overline{A}}{\widehat{\theta }}_i^{\text{EBLUP}}},$ agree with the estimated totals $\widehat{Y}={\displaystyle {\sum }_{i\in A}{\widehat{Y}}_i}$ at the overall sample level $s;$ i.e., $\widehat{\theta }=\widehat{Y}.$ In the case of a mean, ${\theta }_i={\overline{Y}}_i,$ this benchmarking condition becomes ${\displaystyle {\sum }_{i\in A\cup \overline{A}}{N}_i{\widehat{\theta }}_i^{\text{EBLUP}}}={\displaystyle {\sum }_{i\in A}{N}_i{\widehat{\theta }}_i},$ where ${\widehat{\theta }}_i={\widehat{\overline{Y}}}_i.$

Two methods are available in the production system to ensure benchmarking for area based small area estimates. The first one is based on a difference adjustment and the second one is based on an augmented vector. They are valid for any method used to compute ${\widehat{\theta }}_i^{\text{EBLUP}}$ or whether the variance estimate ${\ddot{\psi }}_i$ has been smoothed or not. The benchmarking based on a difference adjustment is an adaptation of the benchmarking given in Battese et al. (1988). The benchmarking based on an augmented vector is due to Wang, Fuller and Qu (2008).

Difference adjustment: For this method, the ${\widehat{\theta }}_i^{\text{EBLUP}}$ estimator is adjusted only for those areas where the realized sample size $n_i\ge 1,i\in A$ and the synthetic estimates $z_i^T\widehat{\beta }$ for $i\in \overline{A}$ are left as is. The resulting benchmarked estimator is given by ${\widehat{\theta }}_i^{\text{EBLUP},b}$ and is defined as follows

$${\widehat{\theta }}_i^{\text{EBLUP},b}=\{\begin{array}{ll}{\widehat{\theta }}_i^{\text{EBLUP}}+{\alpha }_i\left({\widehat{\theta }}^{*}-{\displaystyle \sum _{d\in A}{\omega }_d{\widehat{\theta }}_d^{\text{EBLUP}}}\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 U_A}{\omega }_i^2\left({\ddot{\psi }}_i+b_i^2{\widehat{\sigma }}_v^2\right)}\right\}}^{-1}{\omega }_i\left({\ddot{\psi }}_i+b_i^2{\widehat{\sigma }}_v^2\right)$ for $i\in A,$ ${\omega }_i=1,$ if the benchmarking is to a total, and ${\omega }_i=N_i/N,$ if the benchmarking is for the mean. The estimator ${\widehat{\theta }}^{*}$ is 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{EBLUP},b}}={\widehat{\theta }}^{*}.$

Augmented vector: The vector $z_i^T$ is augmented with ${\omega }_i{\ddot{\psi }}_i,$ to form $z_i^{*T}=\left(z_i^T\text{},{\omega }_i{\ddot{\psi }}_i\right)$ with ${\omega }_i$ and ${\ddot{\psi }}_i$ as previously defined. The resulting augmented generalized linear mixed model (GLMM) equation is given by

$${\widehat{\theta }}_i=z_i^{*T}{\beta }^{\text{*}}+b_i{v}_i^{*}+e_i(3.5)$$

where $E_m\left(v_i^{*}\right)=0$ and $V_m\left(v_i^{*}\right)={\sigma }_v^{*2}.$ The estimates for ${\beta }^{*}$ and ${\sigma }_v^{*2}$ are once more solved recursively for the four EBLUP procedures that we denote as ${\widehat{\theta }}_i^{\text{EBLUP*}}.$

The resulting benchmarked estimator ${\widehat{\theta }}_i^{\text{EBLUP}*,b}$ is given by

$${\widehat{\theta }}_i^{{\text{EBLUP}}^{*}\text{},b}=\{\begin{array}{ll}{\widehat{\gamma }}_i^{*}{\widehat{\theta }}_i^{{\text{EBLUP}}^{*}}+\left(1-{\widehat{\gamma }}_i^{*}\right)z_i^{*T}{\widehat{\beta }}^{\text{*}}\hfill & \text{for}i\in A\hfill \\ z_i^{*T}{\widehat{\beta }}^{\text{*}}\hfill & \text{for}i\in \overline{A}\hfill \end{array}$$

where ${\widehat{\gamma }}_i^{*}=\left(b_i^2{\widehat{\sigma }}_v^{*2}\right)/\left({\ddot{\psi }}_i+b_i^2{\widehat{\sigma }}_v^{*2}\right),$ and ${\widehat{\beta }}^{*}={\left({\displaystyle {\sum }_{i\in A}{z}_i^{*}{z}_i^{*T}}/\left({\ddot{\psi }}_i+b_i^2{\widehat{\sigma }}_v^{*2}\right)\right)}^{-1}{\displaystyle {\sum }_{i\in A}{z}_i^{*}{\widehat{\theta }}_i}/\left({\ddot{\psi }}_i+b_i^2{\widehat{\sigma }}_v^{*2}\right).$

All the components of ${\widehat{\theta }}_i^{{\text{EBLUP}}^{*}\text{},b}$ are computed using the augmented model given by (3.5). It can be shown that ${\displaystyle {\sum }_{i\in A\cup \overline{A}}{\omega }_i{\widehat{\theta }}_i^{{\text{EBLUP}}^{*}\text{},b}}={\displaystyle {\sum }_{i\in A}{\omega }_i{\widehat{\theta }}_i},$ and hence the benchmarking holds.

The difference adjustment and augmented vector methods are two ways that benchmarking can be satisfied. Wang et al. (2008) suggested other procedures that can be used. Specifically, they adapted the self-calibrated estimator You and Rao (2002) developed in the context of the unit level model to the area level model. You, Rao and Hidiroglou (2013) obtained an estimator of the mean squared prediction error and its bias under a misspecified model.

3.3  Mean squared error estimation

The reliability of the EBLUP estimators is obtained as $\text{MSE}\left({\widehat{\theta }}_i^{\text{EBLUP}}\right)=E{\left({\widehat{\theta }}_i^{\text{EBLUP}}-{\theta }_i\right)}^2.$ The expectation is with respect to models (3.3) for the non-benchmarked estimator, and (3.5) for the benchmarked estimator.

The estimated Mean Squared Errors (MSEs) of the area level estimators are given in Table 3.1. The specific form of the $g$ terms and the estimated variances can be found in Rao and Molina (2015) or in Estevao et al. (2015). For the benchmarked estimators, the estimated MSE for the difference adjustment approach uses the non-benchmarked MSE formulas. For the case of the augmented vector approach, the MSE is based on augmenting the vector $z_i^T$ with ${\omega }_i{\ddot{\psi }}_i.$

The various $g$ terms in Table 3.1 can be interpreted as follows. The $g_{0i}$ is a bias correction term for FH and ADM. The $g_{1i}$ term given by $g_{1i}={\widehat{\gamma }}_i{\ddot{\psi }}_i,$ accounts for most of the MSE if the number of areas is large. The $g_{2i}$ term accounts for the estimation of $\beta ,$ and $2g_{3i}$ accounts for the estimation of ${\sigma }_v^2.$ The $g_{4i}$ term in the WF procedure reflects that the estimated value of ${\psi }_i,$ ${\widehat{\psi }}_i,$ has been used. The estimated variance of $\widehat{\beta },$ given by $\text{var}\left(\widehat{\beta }\right)={\left({\displaystyle {\sum }_{i\in A}\frac{z_i{z}_i^T}{{\ddot{\psi }}_i+b_i^2{\widehat{\sigma }}_v^2}}\right)}^{-1}$ is dependent on the particular procedure used to estimate ${\sigma }_v^2.$

• Previous
• Article main page
• Next