The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

saebenchmarking

R-CMD-check

Overview

saebenchmarking provides difference, ratio, and optimum benchmarking for small area estimates, so that their weighted aggregate agrees with the weighted aggregate of the direct estimates, together with mean squared error (MSE) estimation of the benchmarked estimates.

Installation

You can install the development version of saebenchmarking from GitHub with:

# install.packages("pak")
pak::pak("fionaaudia/saebenchmarking")

Example

The package ships with data_eblup, a simulated Fay-Herriot dataset with EBLUP estimates for 50 areas.

library(saebenchmarking)

bm <- sae_benchmarking(
  method   = "difference",
  direct   = "direct",
  weight   = "weight",
  estimate = "eblup",
  mse      = "mse",
  vardir   = "vardir",
  data     = data_eblup
)

summary(bm)
#> Small area benchmarking summary
#> Method          : difference 
#> Number of areas : 50 
#> 
#>                      Direct Estimate    Bench   Target
#> Weighted aggregate 9.307983  9.33162 9.307983 9.307983
head(coef(bm))
#>    Area_1    Area_2    Area_3    Area_4    Area_5    Area_6 
#>  6.964284  8.636035 13.415885 10.476780  9.714784 15.000472

Ratio and optimum benchmarking

bm_ratio <- sae_benchmarking("ratio", direct = "direct", weight = "weight",
                             estimate = "eblup", vardir = "vardir",
                             data = data_eblup)

bm_opt <- sae_benchmarking("optimum", direct = "direct", weight = "weight",
                           estimate = "eblup", mse = "mse",
                           vardir = "vardir", phi_source = "mse",
                           data = data_eblup)

bm_opt$aggregation
#>                      Direct Estimate    Bench   Target
#> Weighted aggregate 9.307983  9.33162 9.307983 9.307983

MSE of the benchmarked estimates

Z   <- cbind(1, data_eblup$z)
s2v <- attr(data_eblup, "sigma2_v")

mse_diff <- mse_benchmarking(bm, estimator = "eblup", z = Z,
                             sigma2_v = s2v, fitting_method = "REML")
head(mse_diff$mse)
#>    Area_1    Area_2    Area_3    Area_4    Area_5    Area_6 
#> 0.2157345 0.2536077 0.2714558 0.2781556 0.3027576 0.3194359

mse_ratio <- mse_benchmarking(bm_ratio, estimator = "eblup", z = Z,
                              sigma2_v = s2v, B = 100, seed = 1)
head(mse_ratio$mse)
#>    Area_1    Area_2    Area_3    Area_4    Area_5    Area_6 
#> 0.2569607 0.2312479 0.2697205 0.3118135 0.2699223 0.3808344

For hierarchical Bayes estimates, supply the posterior means and variances:

bm_hb <- sae_benchmarking("ratio", direct = "direct", weight = "weight",
                          estimate = "theta_hb", data = data_hb)

mse_hb <- mse_benchmarking(bm_hb, estimator = "hb",
                           theta_hb = data_hb$theta_hb,
                           V_hb = data_hb$var_hb)
head(mse_hb$mse)
#>    Area_1    Area_2    Area_3    Area_4    Area_5    Area_6 
#> 0.2145888 0.2378206 0.2616146 0.3014484 0.2893667 0.3220340

Learn more

See vignette("saebenchmarking", package = "saebenchmarking").

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.