Two local diagnostics to evaluate the efficiency of the empirical best predictor under the Fay-Herriot model
Section 6. Simulation study

• Previous
• Next

A simulation study was conducted to evaluate the effectiveness of ${\widehat{D}}_{1i}$ and ${\widehat{D}}_{2i}$ in detecting which of the direct and EB estimators is preferable. We considered $m\mathrm{<mo>=</mo>}$ 140 domains representing Canadian cities. In this simulation study, the vector of auxiliary variables is: $z_i^{{\rm T}}\mathrm{<mo>=</mo>}(\mathrm{1,}{z}_{1i}).$ The auxiliary variable $z_{1i}$ is obtained from administrative files and is defined as the ratio of the number of employment insurance beneficiaries in city $i$ to the number of people over 15 years of age in city $i.$ The sample size in city $i,$ $n_i,$ was obtained from the Canadian Labour Force Survey (LFS). Of the 140 cities, 2 have a sample size smaller than 10, 10 have a sample size smaller than 30, 40 have a sample size smaller than 60, and 68 have a sample size smaller than 100, representing almost 50% of the cities. For these 68 cities, the estimated coefficients of variation of the LFS unemployment rates are in most cases too large to publish direct estimates of the unemployment rate; as a result, small area estimation techniques are required for these domains. In contrast, there are also 17 of the 140 cities with a sample size larger than 1,000 for which the direct estimate of the unemployment rate is reliable.

The population parameter ${\theta }_i$ was simulated for the $m$ domains using the actual values of $n_i$ and $z_i.$ It can be interpreted as the proportion of unemployed people in city $i.$ The parameter ${\theta }_i$ was generated using the beta distribution with mean ${\beta }^{{\rm T}}{z}_i$ and variance ${\sigma }_v^2,$ where ${\sigma }_v^2\mathrm{<mo>=</mo>}\text{7}\text{.58}\times {10}^{-5}$ and ${\beta }^{{\rm T}}\mathrm{<mo>=</mo>}$ (0.0484, 0.95). These values of $\beta$ and ${\sigma }_v^2$ were chosen from real data. We set $b_i\mathrm{<mo>=</mo>}\mathrm{1,}i\mathrm{<mo>=</mo>}\mathrm{1,}\dots ,m.$ Then, we manually changed the values of ${\theta }_i$ for four domains (cities) in order to have a local effect $v_i$ equal to $5{\sigma }_v.$ Cities with different sample sizes were chosen: 10, 100, 501 and 3,773. In the rest of this section, the smallest of these four cities is identified by City 1 $\left(n_i\mathrm{<mo>=</mo>}10\right),$ the second smallest by City 2 $\left(n_i\mathrm{<mo>=</mo>}100\right),$ the second largest by City 3 $\left(n_i\mathrm{<mo>=</mo>\; 501}\right)$ and the largest by City 4 $(n_i\mathrm{<mo>=</mo>}$ 3,773).

A stratified simple random sampling with replacement design was considered where strata coincide with domains. The direct estimator ${\widehat{\theta }}_i$ of ${\theta }_i$ is simply the proportion of sampled people in area $i$ who have the characteristic of interest (e.g., being unemployed). Under such a simple design, it is easy to see that the direct estimator can be generated as follows: ${\widehat{\theta }}_i\mathrm{<mo>=</mo>}{n}_i^{-1}$ $\text{Binomial}\left(n_i\mathrm{,}{\theta }_i\right).$ It is therefore not necessary to create the population of people in domain $i$ to generate ${\widehat{\theta }}_i.$ We proceeded in this way in the simulation. The design variance of ${\widehat{\theta }}_i$ is given by ${\psi }_i\mathrm{<mo>=</mo>}{n}_i^{-1}{\theta }_i\left(1-{\theta }_i\right)$ and its estimator by ${\widehat{\psi }}_i\mathrm{<mo>=</mo>}{\left(n_i-1\right)}^{-1}{\widehat{\theta }}_i\left(1-{\widehat{\theta }}_i\right).$ The smoothed variance ${\tilde{\psi }}_i$ is estimated using the smoothing model in Section 5 with $x_i^{{\rm T}}\mathrm{<mo>=</mo>}\left(\mathrm{1,}\log \left(z_{1i}\right)\mathrm{,}\log \left(1-z_{1i}\right)\mathrm{,}\log \left(n_i\right)\right).$

In order to simulate a realistic scenario, the underlying assumptions of the Fay-Herriot model are not entirely satisfied in our simulation. For example, the errors $v_i$ and $e_i$ do not exactly follow normal distributions. We used a beta distribution to generate ${\theta }_i;$ the normality assumption of $v_i$ is therefore not satisfied although the deviation from the normal distribution is not severe in our simulation. The estimates ${\widehat{\theta }}_i$ were generated from a binomial distribution, which can be approximated by a normal distribution for domains with a large value of $n_i.$ The relationship between the simulated estimates ${\widehat{\theta }}_i$ and the auxiliary vectors $z_i$ is similar to the one observed with the real LFS estimates. Moreover, our simulation scenario is such that the assumption ${\tilde{\psi }}_i\mathrm{<mo>=</mo>}{\psi }_i$ is not satisfied since, for this simple design,

                                                        $${\tilde{\psi }}_i\mathrm{<mo>=</mo>}{n}_i^{-1}\left({\beta }^{{\rm T}}{z}_i(1-{\beta }^{{\rm T}}{z}_i)-b_i^2{\sigma }_v^2\right)$$

(see remark in Section 2). However, we note that the correlation coefficient between ${\tilde{\psi }}_i$ and ${\psi }_i$ is 0.98, which indicates that the deviation from the assumption ${\tilde{\psi }}_i\mathrm{<mo>=</mo>}{\psi }_i$ is modest. As mentioned in the previous paragraph, the smoothing model in Section 5 is used to estimate ${\tilde{\psi }}_i.$ This allows us to remain in a realistic framework where the postulated smoothing model is different from the true model used to generate the estimates ${\widehat{\psi }}_i.$

We conducted a design-based simulation study, i.e., the population parameters ${\theta }_i\mathrm{,}i\mathrm{<mo>=</mo>\; 1,}\dots ,m,$ were generated only once. We repeated sample selection $K\mathrm{<mo>=</mo>}$ 10,000 times. For each replicate $k,$ $k\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}K,$ a direct estimate ${\widehat{\theta }}_i\left(k\right)$ was generated and a smoothed variance estimate ${\widehat{\tilde{\psi }}}_i\left(\text{}k\text{}\right)$ was calculated as described above. The EB estimate was then calculated as:

                                                $${\widehat{\theta }}_i^{\text{EB}}(k)\mathrm{<mo>=</mo>}{\widehat{\gamma }}_i(k){\widehat{\theta }}_i(k)+\left(1-{\widehat{\gamma }}_i(k)\right)\widehat{\beta }{(k)}^{{\rm T}}{z}_i\mathrm{,}$$

where ${\widehat{\gamma }}_i\left(\text{}k\text{}\right)\mathrm{<mo>=</mo>}\frac{{\widehat{\sigma }}_v^2\left(\text{}k\text{}\right)}{{\widehat{\sigma }}_v^2\left(\text{}k\text{}\right)+{\widehat{\tilde{\psi }}}_i\left(\text{}k\text{}\right)}$ and $\widehat{\beta }\left(\text{}k\text{}\right)$ and ${\widehat{\sigma }}_v^2\left(\text{}k\text{}\right)$ are calculated as described in Section 5. The generalized least squares method was used to obtain $\widehat{\beta }\left(\text{}k\text{}\right)$ and the restricted maximum likelihood method was used to obtain ${\widehat{\sigma }}_v^2\left(\text{}k\text{}\right).$ Calculations were performed using Statistics Canada’s small area estimation system (Hidiroglou, Beaumont and Yung, 2019).

For each replicate, standardized residuals ${\widehat{\epsilon }}_i\left(\text{}k\text{}\right)$ and diagnostics ${\widehat{D}}_{1i}\left(\text{}k\text{}\right)$ and ${\widehat{D}}_{2i}\left(\text{}k\text{}\right)$ were also calculated for the $m$ domains. We recorded whether the direct estimator was preferred over the EB estimator for each of the two diagnostics. Decision thresholds were used for this purpose. Below the thresholds, the direct estimator is used. For Diagnostic 1, thresholds of 50% and 75% were used and for Diagnostic 2, thresholds of 5% and 25% were used.

From the previous quantities, calculated for each of the 10,000 replicates, the Monte Carlo averages of Diagnostics 1 and 2 were calculated for the $m$ domains: ${\overline{\widehat{D}}}_{1i}$ and ${\overline{\widehat{D}}}_{2i}.$ The selection rate of the direct estimator was also calculated for each of the two diagnostics, i.e., the percentage of times a given diagnostic led to the selection of the direct estimator.

The Monte Carlo approximation of ${\text{MSE}}_p({\widehat{\theta }}_i^{\text{EB}})$ was calculated as:

                                                    $${\text{MSE}}_{\text{MC}}({\widehat{\theta }}_i^{\text{EB}})\mathrm{<mo>=</mo>}\frac{1}{K}{\displaystyle \sum _{k\mathrm{<mo>=</mo>\; 1}}^K{\left({\widehat{\theta }}_i^{\text{EB}}(k)-{\theta }_i\right)}^2}\mathrm{.}$$

From this Monte Carlo MSE, the relative efficiency of the EB estimator was calculated as:

                                                    $$\frac{{\text{MSE}}_{\text{MC}}({\widehat{\theta }}_i^{\text{EB}})-{\psi }_i}{{\psi }_i}\mathrm{.}(6.1)$$

This ratio is positive when the EB estimator is less efficient than the direct estimator under the design. A diagnostic is potentially useful if it is negatively correlated with this ratio.

Figures 6.1 and 6.2 present the Monte Carlo averages of Diagnostic 1 and 2 respectively as a function of the relative efficiency of the EB estimator defined in equation (6.1). The four cities whose values of $v_i$ have been changed, Cities 1 to 4, are shown in purple, orange, green and red. In the legend, the sample size of these cities has been indicated. The values of the parameter ${\gamma }_i$ for Cities 1 to 4 are 0.01, 0.08, 0.35 and 0.81 respectively. All other cities are shown in blue.

First, we can see in Figures 6.1 and 6.2 that the EB estimator is more efficient than the direct estimator for City 1 (in purple) since this city is to the left of the vertical line (negative relative efficiency) despite the strong local effect. The explanation of this phenomenon is shown in Figure 3.1. It shows that the range of values of $v_i$ for which the B estimator is more efficient than the direct estimator increases as ${\gamma }_i$ decreases. Since ${\gamma }_i$ is small for City 1 $\left({\gamma }_i\mathrm{<mo>=</mo>}\text{0}\text{.01}\right),$ it is not surprising to observe a negative relative efficiency despite a pronounced local effect. For City 2 (in orange), the direct estimator is slightly more efficient than the EB estimator. On the other hand, for Cities 3 (in green) and 4 (in red), the direct estimator is much more efficient than the EB estimator. Note also that there are five cities for which the direct estimator is more efficient than the EB estimator: Cities 2 to 4 as well as two other cities whose values of $v_i$ were randomly generated and not manually modified. One of these cities has the smallest value of $v_i$ and the other has the largest value of $v_i$ after excluding the four cities that had their value manually modified. These two cities have large values of ${\gamma }_i$ (0.62 and 0.49).

Figures 6.1 and 6.2 indicate that our two diagnostics seem to be quite effective in detecting cases where the direct estimator is more efficient than the EB estimator except for City 2 $\left(n_i\mathrm{<mo>=</mo>}100\right)$ where the Monte Carlo average of Diagnostic 1 is very high at 0.97. However, this is a domain where choosing the least efficient estimator is not really problematic since there is very little difference between the efficiencies of the two estimators. Apart from this specific case, Diagnostic 1 seems to have better properties than Diagnostic 2. The Monte Carlo average of Diagnostic 1 is very close to 1 when the EB estimator is significantly more efficient than the direct estimator, decreases slowly when the efficiencies of the two estimators approach each other and becomes small when the direct estimator is significantly more efficient than the EB estimator. Not exactly the same behaviour is observed for Diagnostic 2. The Monte Carlo average of Diagnostic 2 is small when the direct estimator is significantly more efficient than the EB estimator but it is not close to 1 when the EB estimator is significantly more efficient than the direct estimator. Furthermore, it seems to increase when the efficiencies of the two estimators come closer, which is counterintuitive.

Figures 6.3 and 6.4 show the selection rate of the direct estimator over the 10,000 replicates for Diagnostics 1 and 2. Similar conclusions can be drawn as those obtained by analyzing Figures 6.1 and 6.2. As expected, the thresholds of 75% for Diagnostic 1 and 25% for Diagnostic 2 allow better detection of cases where the direct estimator is more efficient than the EB estimator, but these thresholds also lead to the direct estimator being chosen a little too often when it was less efficient than the EB estimator. This is particularly notable for Diagnostic 2. This error can be dampened by decreasing the thresholds, but this also reduces the selection rate of the direct estimator when it is more efficient than the EB estimator. As noted earlier, Diagnostic 1 appears to have better properties than Diagnostic 2, regardless of the thresholds chosen, with very small selection rates of the direct estimator when it is significantly less efficient than the EB estimator. This seems to show the limitations of a fully design-based approach, such as the one presented in Section 4.2, to address the challenge of small domain sample sizes.

• Previous
• Article main page
• Next