Model-based small area estimates, such as the EBLUP under the Fay-Herriot model, usually do not add up to the reliable direct estimate of the larger area. Benchmarking adjusts them so that \(\sum_i w_i \hat\theta_i^{B} = \sum_i w_i \hat\theta_i^{DIR}\).
data_eblup contains simulated data for 50 areas with
EBLUP estimates; data_hb contains the same areas with
hierarchical Bayes estimates.
head(data_eblup[, c("area", "weight", "direct", "vardir", "eblup", "mse")])
#> area weight direct vardir eblup mse
#> 1 1 0.11968123 6.742670 0.2758621 6.987922 0.2143232
#> 2 2 0.07528106 8.720773 0.3478261 8.659673 0.2521964
#> 3 3 0.06887719 13.266346 0.3636364 13.439523 0.2700445
#> 4 4 0.05692330 11.149910 0.4000000 10.500418 0.2767443
#> 5 5 0.04358190 9.840026 0.4571429 9.738421 0.3013463
#> 6 6 0.04358190 15.581303 0.4571429 15.024109 0.3180245methods <- c("difference", "ratio", "optimum")
bench <- lapply(methods, function(m) {
sae_benchmarking(method = m, direct = "direct", weight = "weight",
estimate = "eblup", mse = "mse", vardir = "vardir",
phi_source = "mse", data = data_eblup)
})
names(bench) <- methods
sapply(bench, function(x) x$aggregation[, c("Estimate", "Bench", "Target")])
#> difference ratio optimum
#> Estimate 9.331620 9.331620 9.331620
#> Bench 9.307983 9.307983 9.307983
#> Target 9.307983 9.307983 9.307983The difference method shifts every area by the same amount, the ratio method scales every area by the same factor, and the optimum method distributes the discrepancy according to the MSE (or sampling variance) of each area.
Z <- cbind(1, data_eblup$z)
s2v <- attr(data_eblup, "sigma2_v")
mse_diff <- mse_benchmarking(bench$difference, estimator = "eblup",
z = Z, sigma2_v = s2v, fitting_method = "REML")
mse_opt <- mse_benchmarking(bench$optimum, estimator = "eblup",
z = Z, sigma2_v = s2v, B = 100, seed = 2026)
comparison <- data.frame(
EBLUP = data_eblup$mse,
Difference = mse_diff$mse,
Optimum = mse_opt$mse
)
head(round(comparison, 4))
#> EBLUP Difference Optimum
#> Area_1 0.2143 0.2157 0.2193
#> Area_2 0.2522 0.2536 0.2300
#> Area_3 0.2700 0.2715 0.2475
#> Area_4 0.2767 0.2782 0.2746
#> Area_5 0.3013 0.3028 0.3422
#> Area_6 0.3180 0.3194 0.3626A small number of bootstrap replicates is used here to keep the
vignette fast; use a larger B (for example the default of
1000) in practice.
bm_hb <- sae_benchmarking("optimum", direct = "direct", weight = "weight",
estimate = "theta_hb", vardir = "vardir",
phi_source = "vardir", data = data_hb)
mse_hb <- mse_benchmarking(bm_hb, estimator = "hb",
theta_hb = data_hb$theta_hb, V_hb = data_hb$var_hb)
head(data.frame(V_hb = data_hb$var_hb, PMSE = mse_hb$mse))
#> V_hb PMSE
#> Area_1 0.2144676 0.2151610
#> Area_2 0.2376328 0.2380690
#> Area_3 0.2611638 0.2615628
#> Area_4 0.3011694 0.3014992
#> Area_5 0.2891295 0.2893820
#> Area_6 0.3214677 0.3217202Datta, G. S., Ghosh, M., Steorts, R. and Maples, J. (2011). Bayesian benchmarking with applications to small area estimation. TEST, 20(3), 574-588.
Rao, J. N. K. and Molina, I. (2015). Small Area Estimation, 2nd edition. Wiley.
Steorts, R. C. and Ghosh, M. (2013). On estimation of mean squared errors of benchmarked empirical Bayes estimators. Statistica Sinica, 23(2), 749-767.
Sugasawa, S., Tamae, H. and Kubokawa, T. (2017). Bayesian estimators for small area models shrinking both means and variances. Scandinavian Journal of Statistics, 44(1), 150-167.
Wang, J., Fuller, W. A. and Qu, Y. (2008). Small area estimation under a restriction. Survey Methodology, 34(1), 29-36.