--- title: "Bayesian Area-Level Small Area Estimation with INLA" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Bayesian Area-Level Small Area Estimation with INLA} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5, # INLA is not on CRAN: the code is shown but not run without it eval = requireNamespace("INLA", quietly = TRUE) ) ``` ## Introduction 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, and highly accurate analytical posterior approximations without the computational bottlenecks, convergence diagnostics, or chain-tuning associated with Markov Chain Monte Carlo (MCMC). --- ## Supported Distribution Families `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 | --- ## Mathematical Formulation 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: ### Level 1: Sampling Model (Likelihood) 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$ | ### Level 2: Prior on Fixed Effects $$\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. ### Level 3: Prior on Random Effects $$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.). ### INLA Approximation 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 problems - **`strategy = "laplace"`**: Default, good accuracy/speed balance - **`strategy = "simplified.laplace"`**: Slightly faster Laplace, good for large models The same models can be fitted by MCMC on CPU or GPU with `engine = jax_engine()` (NumPyro NUTS on JAX; see `?jax_engine`). --- ## Spatial and Spatio-Temporal Structures ### Spatial Priors 1. **BYM2 (`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. 2. **Besag ICAR (`spatial = "besag"`)**: Intrinsic Conditional Autoregressive spatial model. 3. **Spatial Lag Model (`spatial = "slm"`)**: Simultaneous spatial autoregression $(I - \rho W)^{-1}$. ### Temporal Models 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. ### Space-Time Interaction Structures (`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"`**: Full Knorr-Held (2000) classifications of space-time interactions. --- ## Case Study: Spatio-Temporal Beta SAE on Emilia-Romagna Poverty 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. ### 1. Load Data and Construct Spatial Adjacency Matrix ```{r emilia_setup, message = FALSE, warning = FALSE} 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")]) ``` ### 2. Fit Spatio-Temporal Beta Model We 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`: ```{r fit_beta_st, message = FALSE, warning = FALSE} 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) ``` ### 3. Inspect Posterior Predictions and Uncertainty The output object contains full domain-period predictions, posterior standard errors, and 95% Credible Intervals: ```{r inspect_predictions} head(fit_beta_st$df_hb[, c("domain", "time", "y", "hb", "sd", "rse", "ci_lower", "ci_upper")]) ``` --- ## Comparison and Validation Against `tipsae` (Stan MCMC) In `tipsae::fit_sae(..., spatial_error = TRUE, temporal_error = TRUE)`, the model is estimated using Hamiltonian Monte Carlo (HMC) in Stan. In `fastsae::hb_area()`, the identical structural model is estimated analytically via INLA. | Characteristic | `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 match | | **Estimation Method** | MCMC (NUTS) | Integrated Nested Laplace Approximations | Deterministic, no burn-in needed | | **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 estimates | | **Mean Absolute Error (MAE)** | Baseline | **0.00304** | Negligible numerical discrepancy | ### Interactive Scalability Explorer ($n = 30$ to $n = 1,000$) 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: ```{=html}
Execution time (median, seconds) — Standard Beta · log scale
Speedup Factor & Efficiency Multiplier (fastsae vs tipsae) — Standard Beta
``` #### Inferential Equivalence & Validation Metrics ($n = 1,000$) 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 | ------------------------------------------------------------------------ ## Cross-Sectional Count Data: Poisson & Negative Binomial For count data such as disease incidence or crime events, `hb_area()` seamlessly supports Poisson and Negative Binomial distributions with population offsets: ```{r count_example, message = FALSE, warning = FALSE} # Simulate 42-domain count data with known spatial structure data("sim_area", package = "fastsae") data("mys_proxmat", package = "fastsae") # 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 = mys_proxmat, print_result = FALSE ) summary(fit_pois) ``` --- ## References - De Nicolò, S., & Gardini, A. (2022). tipsae: Mapping proportions and rates using spatio-temporal Beta small area models. *Journal of Statistical Software*. - Janicki, H. (2020). Properties of the Beta-logistic model for small area estimation. *Survey Methodology*, 46(1), 89–112. - Knorr-Held, L. (2000). Bayesian modelling of inhomogeneous spatial and temporal variation in rates. *Statistics in Medicine*, 19(17-18), 2555–2567. - Rue, H., Martino, S., & Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. *Journal of the Royal Statistical Society: Series B*, 71(2), 319–392. - Simpson, D., Rue, H., Riebler, A., Martins, T. G., & Sørbye, S. H. (2017). Penalising model component complexity: A principled, practical approach to constructing priors. *Statistical Science*, 32(1), 1–28.