Classical Small Area Estimation (SAE) models—such as the Fay-Herriot (1979) model—traditionally assume a Gaussian distribution for direct survey estimators with known sampling error variances. However, real-world official statistics frequently encounter non-Gaussian indicators: - Poverty rates, unemployment rates, and prevalence ratios: strictly bounded in the open interval \((0, 1)\). - Disease counts, crime incidents, and rare events: discrete counts subject to Poisson or Negative Binomial overdispersion with population exposures. - Binary survey aggregates: proportions derived from small sample sizes (Binomial trials). - Skewed economic expenditures: strictly positive continuous outcomes (Gamma).
Furthermore, when data are collected across geographical areas and repeated over time (panel/longitudinal surveys), outcomes exhibit both spatial dependency and temporal autocorrelation.
The fastsae package introduces
hb_area(), a unified hierarchical Bayesian framework
powered by Integrated Nested Laplace Approximations
(INLA). It provides fast deterministic approximations of the
posterior marginals, without the chains, warm-up and tuning of Markov
chain Monte Carlo (MCMC); against long Stan runs with the same model and
priors its posterior means and SDs agreed closely (see the Engines
section of ?hb_area).
hb_area() supports 6 probability distributions for the
direct response:
| Family | Typical SAE Use Case | Required Sampling Input | Model Scale Parameter |
|---|---|---|---|
"gaussian" |
Continuous indicators (e.g., mean income, mean years of schooling) | vardir (sampling variance \(D_d\)) |
Inverted sampling variance \(\text{scale} = 1/D_d\) |
"beta" |
Rates and proportions strictly in \((0, 1)\) | vardir or trials |
Sampling dispersion via Janicki (2020): \(\phi_d = \frac{y_d(1 - y_d)}{V_d} - 1\) |
"poisson" |
Counts and rare event rates | exposure (expected count / population
offset) |
Log-link offset \(\log(E_d)\) |
"nbinomial" |
Overdispersed counts and clustered events | exposure |
Log-link offset with negative binomial dispersion |
"binomial" |
Aggregated binary proportions | trials (domain sample size \(n_d\)) |
Number of trials per domain |
"gamma" |
Positive, right-skewed indicators (e.g., expenditures) | vardir |
Shape-rate parameterization with known sampling variance |
The hierarchical Bayesian model underlying hb_area()
extends the classical Fay-Herriot framework to handle non-Gaussian
responses and structured random effects. The model comprises three
levels:
For domain \(d\) with direct estimator \(y_d\) and known sampling variance \(D_d\):
\[\pi(y_d \mid \eta_d) = \text{ExpFamily}(y_d \mid \eta_d, D_d)\]
The link function connects the linear predictor \(\eta_d\) to the mean \(\mu_d = E(y_d)\):
| Family | Link Function | Linear Predictor |
|---|---|---|
| Gaussian | Identity: \(\eta_d = \mu_d\) | \(\eta_d = x_d^\top \beta + u_d\) |
| Beta | Logit: \(\eta_d = \log(\mu_d / (1-\mu_d))\) | \(\eta_d = x_d^\top \beta + u_d\) |
| Poisson | Log: \(\eta_d = \log(\mu_d)\) | \(\eta_d = x_d^\top \beta + u_d + \log(E_d)\) |
| Binomial | Logit: \(\eta_d = \log(\mu_d / (1-\mu_d))\) | \(\eta_d = x_d^\top \beta + u_d\) |
| Negative Binomial | Log: \(\eta_d = \log(\mu_d)\) | \(\eta_d = x_d^\top \beta + u_d + \log(E_d)\) |
| Gamma | Log: \(\eta_d = \log(\mu_d)\) | \(\eta_d = x_d^\top \beta + u_d\) |
\[\beta \sim N(0, \tau_\beta^2 I_p)\]
A weakly informative normal prior with large variance \(\tau_\beta^2\) is placed on regression coefficients.
\[u \sim N(0, \Sigma)\]
The random effects \(u = (u_1, \ldots, u_m)^\top\) follow a zero-mean Gaussian Markov Random Field (GMRF) with precision matrix \(\Sigma^{-1}\) encoding the chosen spatial or temporal structure (BYM2, Besag ICAR, RW1, etc.).
Rather than MCMC sampling, INLA approximates the posterior marginals:
\[\pi(\theta \mid y) \approx \tilde{\pi}(\theta \mid y) = \frac{\pi(\theta) \pi(y \mid \theta)}{\pi(y \mid \theta)}\]
and
\[\pi(u_d \mid y) \approx \tilde{\pi}(u_d \mid y) = \sum_k \tilde{\pi}(u_d \mid \theta_k, y) \tilde{\pi}(\theta_k \mid y) \Delta_k\]
where the integrals are evaluated using Laplace
approximations and numerical integration.
Three computational strategies are available through
hb_area(..., engine = inla_engine(strategy = ...)), which
also takes the priors and any further INLA::inla()
arguments:
strategy = "gaussian": Fastest, best
for well-behaved problemsstrategy = "laplace": Default, good
accuracy/speed balancestrategy = "simplified.laplace":
Slightly faster Laplace, good for large modelsThe same models can be fitted by MCMC on CPU or GPU with
engine = jax_engine() (NumPyro NUTS on JAX; see
?jax_engine).
spatial = "bym2"): The scaled
Besag-York-Mollié model (Simpson et al., 2017) decomposes area random
effects into a spatial structured component and an unstructured IID
component: \[u_d = \frac{1}{\sqrt{\tau_u}}
\left( \sqrt{1 - \phi} \, v_d + \sqrt{\phi} \, u_* d \right)\]
where \(\phi \in [0, 1]\) measures the
proportion of variance explained by spatial structure, regulated by
Penalized Complexity (PC) priors.spatial = "besag"):
Intrinsic Conditional Autoregressive spatial model.spatial = "slm"):
Simultaneous spatial autoregression \((I -
\rho W)^{-1}\).When data are observed longitudinally across \(T \ge 2\) time periods, specify
time: - temporal = "rw1":
First-order Random Walk (\(\Delta t_\tau \sim
N(0, \sigma_t^2)\)). -
temporal = "rw2": Second-order Random Walk
for smooth non-linear temporal trends. -
temporal = "ar1": First-order
Autoregressive process with autocorrelation \(|\rho_t| < 1\). -
temporal = "iid": Unstructured time random
effects.
st_interaction)"none" (Additive): Separate spatial
main effect and temporal trend: \[\eta_{dt} =
x_{dt}^\top \beta + s_d + \gamma_t\]"domain-specific" (matching
tipsae): Domain-specific random walks / AR(1) processes
with shared precision: \[\eta_{dt} =
x_{dt}^\top \beta + s_d + t_{d, t}, \quad t_{d, \cdot} \sim
\text{RW1}(\sigma_t^2)\]"separable" (matching classical
eblup_stfh): Spatial random field evolving dynamically
across time with AR(1) temporal correlation."type1" to "type4":
Knorr-Held (2000) interactions \(\delta_{dt}\) added to the main effects,
\(\eta_{dt} = x_{dt}^\top \beta + s_d +
\gamma_t + \delta_{dt}\):
"type1": \(\delta_{dt}\) independent (unstructured
space \(\times\) unstructured
time);"type2": an independent temporal process
for each domain (unstructured space \(\times\) structured time);"type3": an independent ICAR field for each period
(structured space \(\times\)
unstructured time);"type4": an ICAR field evolving with the
temporal structure (structured space \(\times\) structured time)."type3" and "type4" need a structured
spatial effect ("bym2", "bym",
"besag" or "leroux").We demonstrate hb_area() using the Italian poverty
dataset emilia from the tipsae package (38
health districts in Emilia-Romagna across 5 years, \(N = 190\)). The goal is to estimate the
Head Count Ratio (hcr), which represents the proportion of
households below the poverty line.
library(fastsae)
library(tipsae)
library(spdep)
data("emilia")
data("emilia_shp")
# Construct binary spatial adjacency matrix W from district polygons
nb <- poly2nb(emilia_shp)
W <- nb2mat(nb, style = "B", zero.policy = TRUE)
rownames(W) <- colnames(W) <- as.character(emilia_shp$NAME_DISTRICT)
head(emilia[, c("id", "year", "hcr", "vars", "x")])
#> id year hcr vars x
#> 1 CASALECCHIO DI RENO 2014 0.0404 9.090478e-05 -0.2624
#> 2 CITTA' DI BOLOGNA 2014 0.0825 6.404001e-05 -0.0008
#> 3 IMOLA 2014 0.1033 3.120275e-04 -0.0522
#> 4 PIANURA EST 2014 0.0633 1.025764e-04 -0.4007
#> 5 PIANURA OVEST 2014 0.0625 1.562500e-04 -0.2277
#> 6 PORRETTA TERME 2014 0.1276 6.643609e-04 -0.4434We model the poverty proportion using a Beta likelihood with a Besag
ICAR spatial random effect, a domain-specific RW(1) temporal dynamic,
and known sampling variances vars:
fit_beta_st <- hb_area(
formula = hcr ~ x,
data = emilia,
domain = "id",
time = "year",
vardir = "vars",
family = "beta",
spatial = "besag",
temporal = "rw1",
st_interaction = "domain-specific",
W = W,
print_result = FALSE
)
summary(fit_beta_st)
#>
#> Variance Components:
#> sigma2_u: 0.062051
#> sigma2_t (temporal): 0.006912
#>
#> Coefficients:
#> beta std.error zvalue pvalue ci_lower
#> (Intercept) -2.2517e+00 1.5232e-02 -1.4783e+02 0.0000e+00 -2.2815e+00
#> x 4.0613e-01 6.4372e-02 6.3090e+00 2.8079e-10 2.7947e-01
#> ci_upper
#> (Intercept) -2.2218
#> x 0.5321
#>
#> Hyperparameters:
#> mean sd 0.025quant 0.5quant 0.975quant
#> Precision for ..domain_id.. 18.223 6.626455 8.634412 17.1022 34.35005
#> Precision for ..time_id.. 172.314 78.436863 72.151498 155.4507 373.05592
#> mode
#> Precision for ..domain_id.. 15.0637
#> Precision for ..time_id.. 127.0839
#>
#> Goodness of Fit:
#> DIC
#> -976.34164
#> pD
#> 56.02335
#> WAIC
#> -984.78951
#> pWAIC
#> 38.70805
#> Marginal_LogLik.log marginal-likelihood (integration)
#> 439.32544
#>
#> HB Summary Statistics:
#> hb linear_pred sd mse
#> Min. :0.05270 Min. :-2.894 Min. :0.005185 Min. :2.688e-05
#> 1st Qu.:0.08271 1st Qu.:-2.413 1st Qu.:0.007702 1st Qu.:5.932e-05
#> Median :0.09635 Median :-2.244 Median :0.009371 Median :8.781e-05
#> Mean :0.09802 Mean :-2.251 Mean :0.009566 Mean :9.767e-05
#> 3rd Qu.:0.11011 3rd Qu.:-2.095 3rd Qu.:0.011051 3rd Qu.:1.221e-04
#> Max. :0.15657 Max. :-1.687 Max. :0.016768 Max. :2.812e-04
#> rse
#> Min. : 6.370
#> 1st Qu.: 8.725
#> Median : 9.812
#> Mean : 9.827
#> 3rd Qu.:10.979
#> Max. :13.670The output object contains full domain-period predictions, posterior standard errors, and 95% Credible Intervals:
head(fit_beta_st$df_hb[, c("domain", "time", "y", "hb", "sd", "rse", "ci_lower", "ci_upper")])
#> domain time y hb sd rse ci_lower
#> 1 CASALECCHIO DI RENO 2014 0.0404 0.05554654 0.006127590 11.031452 0.04408741
#> 2 CITTA' DI BOLOGNA 2014 0.0825 0.08225899 0.006131546 7.453952 0.07068307
#> 3 IMOLA 2014 0.1033 0.09343252 0.010214679 10.932681 0.07511427
#> 4 PIANURA EST 2014 0.0633 0.06509831 0.006294766 9.669629 0.05342031
#> 5 PIANURA OVEST 2014 0.0625 0.07112658 0.007870632 11.065669 0.05671280
#> 6 PORRETTA TERME 2014 0.1276 0.09608716 0.012916400 13.442378 0.07314330
#> ci_upper
#> 1 0.06810577
#> 2 0.09473234
#> 3 0.11517419
#> 4 0.07814088
#> 5 0.08759936
#> 6 0.12376122tipsae (Stan
MCMC)In
tipsae::fit_sae(..., spatial_error = TRUE, temporal_error = TRUE),
the model is estimated by Hamiltonian Monte Carlo (HMC) in Stan.
fastsae::hb_area() fits a model with the same structure by
INLA; the priors of the two packages differ.
| Characteristic | tipsae::fit_sae (Stan MCMC) |
fastsae::hb_area (INLA) |
|---|---|---|
| Model | Beta, Besag ICAR + domain-specific RW(1) | spatial = "besag",
temporal = "rw1",
st_interaction = "domain-specific" |
| Estimation | MCMC (NUTS) | Integrated nested Laplace approximation (INLA’s default
"compact" mode for Beta) |
| Execution time | 120 s (4 chains x 2000 iterations, largest \(\hat R\) 1.07); one chain of 200 iterations (about 10 s) has not converged | 1.1 s |
Posterior means vs converged
tipsae |
Reference | \(r = 0.998\), mean absolute difference 0.0016 |
Measured on 2026-10-16 on one machine (Apple silicon); times depend
on the hardware. The accuracy of hb_area() was also checked
against Stan with identical models and priors; see the Engines section
of ?hb_area.
Use the interactive controls below to compare runtime, memory
footprint, and speedup factors between fastsae::hb_area and
tipsae::fit_sae across both Standard Beta
and Spatial Beta (Besag ICAR) 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 count data such as disease incidence or crime events,
hb_area() seamlessly supports Poisson and Negative Binomial
distributions with population offsets:
# Simulate 42-domain count data with known spatial structure
data("sim_area", package = "fastsae")
# The BYM2 model uses a neighbourhood graph: the 4 nearest centroids of each domain
W_knn <- create_weights(cbind(sim_area$x_coord, sim_area$y_coord), method = "knn", k = 4)
# Fit Poisson Spatial Model with population exposure offset
fit_pois <- hb_area(
formula = y_poisson ~ x1 + x2,
data = sim_area,
exposure = "exposure",
family = "poisson",
spatial = "bym2",
W = W_knn,
print_result = FALSE
)
summary(fit_pois)
#>
#> Variance Components:
#> sigma2_u: 0.06618
#> rho (spatial): 0.152
#> phi (spatial fraction): 0.152
#>
#> Coefficients:
#> beta std.error zvalue pvalue ci_lower
#> (Intercept) 1.4181e-01 1.3561e-01 1.0457e+00 2.9568e-01 -1.2539e-01
#> x1 3.8631e-01 4.9074e-02 7.8719e+00 3.4938e-15 2.8964e-01
#> x2 -2.9856e-01 3.2988e-02 -9.0507e+00 1.4208e-19 -3.6392e-01
#> ci_upper
#> (Intercept) 0.4094
#> x1 0.4831
#> x2 -0.2339
#>
#> Hyperparameters:
#> mean sd 0.025quant 0.5quant
#> Precision for ..domain_id.. 16.3138062 4.5697066 9.02117042 15.749009
#> Phi for ..domain_id.. 0.1519622 0.1492371 0.00706639 0.100582
#> 0.975quant mode
#> Precision for ..domain_id.. 26.8605449 14.70759276
#> Phi for ..domain_id.. 0.5678022 0.01687449
#>
#> Goodness of Fit:
#> DIC
#> 317.75652
#> pD
#> 31.80732
#> WAIC
#> 309.86484
#> pWAIC
#> 17.22853
#> Marginal_LogLik.log marginal-likelihood (integration)
#> -185.30458
#>
#> HB Summary Statistics:
#> hb linear_pred sd mse
#> Min. :0.3072 Min. :-1.1919 Min. :0.02957 Min. :0.0008745
#> 1st Qu.:0.7646 1st Qu.:-0.2790 1st Qu.:0.05802 1st Qu.:0.0033665
#> Median :1.3848 Median : 0.3238 Median :0.07676 Median :0.0058921
#> Mean :1.5290 Mean : 0.2384 Mean :0.16054 Mean :0.0659981
#> 3rd Qu.:2.2984 3rd Qu.: 0.8022 3rd Qu.:0.12181 3rd Qu.:0.0148494
#> Max. :3.6013 Max. : 1.2809 Max. :0.74697 Max. :0.5579689
#> rse
#> Min. : 2.790
#> 1st Qu.: 5.159
#> Median : 6.259
#> Mean :10.346
#> 3rd Qu.: 9.880
#> Max. :29.445