The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
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[@]dogado.de.

Getting Started with fastsae

Introduction

Small Area Estimation (SAE) encompasses statistical techniques designed to produce reliable estimates for sub-populations or geographical domains where sample sizes are too small for direct survey estimators to achieve acceptable precision.

The fastsae package provides high-performance C++ implementations (via Rcpp and RcppArmadillo) for standard and advanced SAE models. It offers:

The Fay-Herriot Model

The area-level model introduced by Fay and Herriot (1979) links direct survey estimators \(y_d\) with auxiliary variables \(x_d\):

\[y_d = x_d^\top \beta + u_d + e_d, \quad d = 1, \dots, D\]

where: - \(u_d \sim \text{i.i.d. } N(0, \sigma_u^2)\) represents domain-specific random effects. - \(e_d \sim \text{ind. } N(0, D_d)\) represents sampling errors with known sampling variance \(D_d\) (vardir).

The Empirical Best Linear Unbiased Predictor (EBLUP) is a weighted combination of the direct estimator and the regression-synthetic estimator:

\[\hat{\theta}_d = \gamma_d y_d + (1 - \gamma_d) x_d^\top \hat{\beta}\]

where \(\gamma_d = \frac{\hat{\sigma}_u^2}{\hat{\sigma}_u^2 + D_d}\) is the shrinkage factor (\(0 \le \gamma_d \le 1\)).

Step-by-Step Example

1. Load Package and Dataset

We use the built-in mys dataset (mean years of schooling):

library(fastsae)
library(ggplot2)

data("mys")
head(mys)
#> # A tibble: 6 × 9
#>    area     y vardir   rse    x1    x2    x3      n  weight
#>   <int> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>   <dbl>
#> 1     1  8.36  0.662  9.73   124    24    14  7280. 0.0326 
#> 2     2  7.60  0.837 12.0     89    18     9  2743. 0.0123 
#> 3     3  5.51  0.882 17.0     57    19     5  1706. 0.00764
#> 4     4  3.87  0.658 21.0     88    35    19  3073. 0.0138 
#> 5     5  6.31  1.28  17.9    141    46    29 13400. 0.0600 
#> 6     6  3.93  0.388 15.9     96    29    10  2004. 0.00897

2. Fit Fay-Herriot Model (eblup_fh)

To fit an area-level Fay-Herriot model using Restricted Maximum Likelihood (REML):

# Fit Fay-Herriot model
fit_fh <- eblup_fh(
  formula = y ~ x1 + x2 + x3,
  vardir = ~vardir,
  data = mys,
  method = "REML",
  print_result = FALSE
)

3. Model Summary and Coefficients

The standard S3 summary() method provides comprehensive model diagnostics, variance components, and coefficient tests:

summary(fit_fh)
#> 
#> ── Summary of fastsae Fit ──────────────────────────────────────────────────────
#> Call :
#> eblup_fh(formula = y ~ x1 + x2 + x3, vardir = ~vardir, data = mys, method =
#> "REML", print_result = FALSE)
#> 
#> ✔ Convergence: Yes (in 6 iterations)
#> Model: Fay-Herriot (Area-level)
#> Method: eblup
#> 
#> Variance Components:
#> sigma2_u: 2.608103 
#> 
#> Coefficients:
#>                   beta  std.error     zvalue pvalue
#> (Intercept)  3.1077510  0.7697687  4.0372527 0.0001
#> x1          -0.0019323  0.0098886 -0.1954019 0.8451
#> x2           0.0555184  0.0614129  0.9040187 0.3660
#> x3           0.0335344  0.0580013  0.5781663 0.5632
#> 
#> Goodness of Fit:
#> loglikelihood           AIC           BIC 
#>     -65.14251     140.28502     147.61370 
#> 
#> EBLUP Summary Statistics:
#>      eblup             mse              rse        
#>  Min.   :0.8506   Min.   :0.2037   Min.   : 9.844  
#>  1st Qu.:4.1973   1st Qu.:0.6886   1st Qu.:14.889  
#>  Median :5.0196   Median :1.0156   Median :22.157  
#>  Mean   :5.0746   Mean   :1.5434   Mean   :26.275  
#>  3rd Qu.:6.1360   3rd Qu.:2.4886   3rd Qu.:35.020  
#>  Max.   :9.0220   Max.   :3.7074   Max.   :53.065

You can extract fixed-effects coefficients using coef():

coef(fit_fh)
#>  (Intercept)           x1           x2           x3 
#>  3.107750953 -0.001932259  0.055518370  0.033534410

Fitted EBLUP estimates and residuals can be extracted using standard generics:

# Fitted values (EBLUP)
head(fitted(fit_fh))
#> [1] 7.612738 6.782316 5.187060 4.201545 6.323679 4.048590

# Residuals (direct estimate - EBLUP)
head(residuals(fit_fh))
#> [1]  0.74678890  0.81733383  0.32707611 -0.33221827 -0.01861583 -0.12178363

4. Diagnostic Plots (autoplot)

fastsae extends ggplot2::autoplot() to provide convenient diagnostic and comparison plots.

Direct Estimates vs EBLUP

Comparing the direct estimates against EBLUP demonstrates shrinkage towards the regression synthetic line:

autoplot(fit_fh, type = "estimates")

Mean Squared Error (MSE) Across Domains

Inspect domain-level uncertainty with MSE plots:

autoplot(fit_fh, type = "mse")

Exact Numerical Equivalence with sae

fastsae produces results that are mathematically identical to sae::eblupFH:

if (requireNamespace("sae", quietly = TRUE)) {
  mys_clean <- as.data.frame(na.omit(mys))
  fit_fast <- eblup_fh(y ~ x1 + x2 + x3, vardir = ~vardir, data = mys_clean, print_result = FALSE)
  fit_sae <- sae::eblupFH(y ~ x1 + x2 + x3, vardir = vardir, data = mys_clean)

  # Check EBLUP estimates
  all.equal(fit_fast$df_eblup$eblup, as.vector(fit_sae$eblup))

  # Check regression coefficients
  all.equal(as.vector(coef(fit_fast)), as.vector(fit_sae$fit$estcoef$beta))

  # Check random effect variance (sigma2_u)
  all.equal(fit_fast$random_effect_var, fit_sae$fit$refvar)
}
#> [1] TRUE

References

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.
Health stats visible at Monitor.