Small Area Estimation often involves large administrative registries or extensive spatial networks with hundreds or thousands of domains. Traditional implementations in R often rely on interpreted loops, pure R Fisher-scoring routines, or large memory allocations for \(D \times D\) dense covariance matrices.
fastsae fits the frequentist models in compiled C++
(RcppArmadillo; sparse Cholesky factorisations for the SAR
models) and the Bayesian models with INLA. At \(n = 1,000\) domains:
sae
and 42,874x faster than emdi.sae and 51x faster than emdi.sae.| Model (\(n = 1,000\) domains, 5 covariates) | fastsae |
sae |
emdi |
|---|---|---|---|
| Fay-Herriot | 0.00119 s (0.503 MB) | 1.50 s (72.0 MB; 1,263x slower) | 51.0 s (3.69 GB; 42,874x slower) |
| Spatial Fay-Herriot | 1.00 s (0.203 MB) | 67.2 s (1.88 GB; 67x slower) | 51.1 s (3.68 GB; 51x slower) |
| Spatio-Temporal FH (\(T = 5\)) | 75.3 s (3.29 MB) | 1790 s (31.7 GB; 24x slower) | - |
Times are medians of bench::mark(); memory is the total
memory allocated by R during one call
(bench::mark()$mem_alloc), not the peak resident memory.
The fastsae rows were re-measured with fastsae 1.0.0
(frequentist on 2026-10-13,
benchmarks/bench_frequentist_fastsae.R; Bayesian on
2026-10-14, benchmarks/bench_bayes_fastsae.R); the
sae/emdi rows come from an earlier run of the
same design (simulated data, 5 covariates) on different random draws,
the tipsae rows from an earlier run on the same data.
Use the interactive controls below to compare execution time, memory allocated, and iterations per second across sample sizes:
Compiled C++ Linear Solvers: Instead of
interpreting nested loops in R, fastsae implements Fisher
scoring and the analytical MSE in Armadillo C++. The spatial Fay-Herriot
model never forms a \(D \times D\)
inverse: with \(A = (I - \rho W')(I - \rho
W)\) and \(M = A + \sigma^2
E'\Psi^{-1}E\), \(V^{-1}\),
\(g_1, \dots, g_4\) and the kriging
terms follow from sparse Cholesky solves with \(A\) and \(M\) (Eigen, AMD ordering).
Zero-Copy Matrix Operations: Memory allocations for large intermediate structures (\(V\), \(V^{-1}\), \(P\)) are avoided or reused across Fisher iterations rather than reallocated on the heap.
OpenMP Multi-Threaded Bootstrap: Bootstrap
resampling runs natively in parallel across CPU cores using
#pragma omp parallel for, avoiding the serialization
overhead of R worker processes (parallel::makeCluster /
foreach).
Woodbury Identity for Panel Data: In
eblup_stfh, inversion of the block \((DT \times DT)\) covariance matrix \(V\) is reduced to operations on individual
\(D \times D\) and \(T \times T\) blocks via Kronecker and
Woodbury decomposition, transforming an \(O((DT)^3)\) bottleneck into scalable \(O(D^3 + T^3)\) steps.
Integrated Nested Laplace Approximations (INLA):
For complex hierarchical generalized and spatio-temporal models
(hb_area), fastsae utilizes INLA to compute
analytical posterior margins directly from sparse Gaussian Markov Random
Fields, bypassing Monte Carlo sampling entirely.
fastsae (INLA) vs
tipsae (Stan MCMC)To assess Bayesian small area estimation performance,
fastsae::hb_area() was benchmarked against
tipsae::fit_sae(), the state-of-the-art Stan MCMC package
for spatio-temporal Beta small area models.
The test used the official tipsae panel dataset:
emilia (38 health districts in
Emilia-Romagna over 5 years, \(N =
190\) domains \(\times\) years)
with spatial polygon contiguity matrix \(W\) from
emilia_shp.
| Metric | tipsae::fit_sae (Stan MCMC) |
fastsae::hb_area (INLA) |
Advantage |
|---|---|---|---|
| Model Specification | Besag ICAR + Domain RW(1) | spatial = "besag",
temporal = "rw1",
st_interaction = "domain-specific" |
Exact structural equivalence |
| Computational Engine | Hamiltonian Monte Carlo (NUTS Stan) | Integrated Nested Laplace Approximation | Analytical & deterministic |
| Execution Time | 15.8 s (1 chain, 200 iter) / ~120 s (4 chains) | 1.93 s
(simplified.laplace) |
~9x to 50x+ faster |
| Pearson Correlation (\(r\)) | Baseline | 0.9893 | Near-identical point estimates |
| Mean Absolute Error (MAE) | Baseline | 0.00304 | Negligible numerical error |
| Convergence Overhead | Requires \(\hat{R} < 1.05\) checks, warmup, and tuning | None (closed-form Laplace expansions) | Instant convergence |
Use the interactive controls below to compare runtime, memory
footprint, and speedup factors between fastsae::hb_area and
tipsae::fit_sae across Standard Beta,
Spatial Beta (Besag ICAR), and Spatio-Temporal
Beta (Besag + RW1) models:
To establish that the substantial computational accelerations
achieved by fastsae::hb_area preserve full inferential
validity relative to Stan MCMC (tipsae::fit_sae), summary
concordance and error metrics were evaluated across \(n = 1,000\) domains/units for Standard
Beta, Spatial Beta (Besag ICAR), and Spatio-Temporal Beta models:
| Dimension | Metric | Beta | Spatial Beta | Spatio-Temporal Beta | Statistical Implication |
|---|---|---|---|---|---|
| Point Estimates (EBP) | Pearson Correlation (\(r\)) | 0.99934 | 0.99891 | 0.99614 | Near-perfect linear agreement between INLA and MCMC |
| Spearman Rank Correlation (\(\rho\)) | 0.99921 | 0.99875 | 0.99517 | Consistent domain priority and rank ordering | |
| Mean Absolute Error (MAE) | 0.01969 | 0.01456 | 0.01200 | Average estimation discrepancy \(< 0.02\) | |
| Root Mean Squared Difference (RMSD) | 0.02369 | 0.01713 | 0.01458 | Negligible \(L_2\) deviation across domains | |
| Uncertainty / MSE | MSE Pearson Correlation (\(r_{\text{MSE}}\)) | 0.96315 | 0.89254 | 0.67957 | Preserved area-level precision ordering |
| Mean Absolute Difference (MSE) | 0.00018 | 0.00012 | 0.00025 | Virtually identical error variances |
For comprehensive design-based simulation studies evaluating parameter recovery and empirical properties under repeated survey sampling from finite populations (\(D = 60\) and national-scale \(D = 500\)), see the dedicated article: Finite Population Simulation & Validation.