Bayesian Area-Level Small Area Estimation with INLA

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 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).

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`	Precision $\phi_d = \hat\mu_d(1 - \hat\mu_d)/V_d - 1$ from a first fit (Janicki, 2020), or $n_d - 1$ from `trials`
`"poisson"`	Counts and rare event rates	`exposure` (expected count / population offset)	Log-link offset $\log(E_d)$
`"nbinomial"`	Overdispersed counts and clustered events	`exposure`; optional `vardir`	Log-link offset; size $\hat\mu_d^2/(V_d - \hat\mu_d)$ from a first fit when `vardir` is given
`"binomial"`	Aggregated binary proportions	`trials` (domain sample size $n_d$ , or the effective sample size $n_d/\mathrm{deff}_d$ for survey-weighted proportions; not rounded)	Number of trials per domain
`"gamma"`	Positive, right-skewed indicators (e.g., expenditures)	`vardir`	Precision (shape) $y_d^2/V_d$

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_0 \propto 1, \qquad \beta_j \sim N(0, 1000), \quad j = 1, \ldots, p$

The intercept has INLA’s flat prior and the other coefficients a normal prior with variance 1000, which is nearly flat on the scale of most covariates.

Level 3: Prior on Random Effects

$u \mid \sigma_u \sim N(0, \sigma_u^2 R^{-}), \qquad \sigma_u \sim \text{half-normal}(s)$

The random effects $u = (u_1, \ldots, u_m)^\top$ follow a zero-mean Gaussian Markov random field whose structure matrix $R$ encodes the chosen spatial or temporal dependence ( $R = I$ for iid effects; Besag, BYM2, RW1 and so on otherwise). The scale of the half-normal prior on their standard deviation is $s = \mathrm{sd}(y)$ , the standard deviation of the direct estimates, for Gaussian models and $s = 1$ on the logit or log scale for the other families. The further hyperparameters (the BYM2 mixing parameter, autocorrelations, dispersion parameters) and the reasons for these defaults are listed in vignette("hb-priors"); inla_engine() replaces any of them.

INLA Approximation

Rather than MCMC sampling, INLA approximates the posterior marginals of the hyperparameters $\theta$ and of the latent field $x = (\beta, u)$ (Rue et al., 2009):

$\tilde{\pi}(\theta \mid y) \propto \left. \frac{\pi(x, \theta, y)}{\tilde{\pi}_G(x \mid \theta, y)} \right|_{x = x^*(\theta)}$

where $\tilde{\pi}_G$ is the Gaussian approximation of the full conditional of $x$ and $x^*(\theta)$ its mode, and

$\tilde{\pi}(x_d \mid y) = \sum_k \tilde{\pi}(x_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). An MCMC fit is only usable when its chains have converged. The fit reports the rank-normalized R-hat and the effective sample size (ESS) of every parameter (Vehtari et al., 2021) and warns when R-hat is at least 1.05, the ESS is below 100 per chain, or there are divergent transitions; the trace and rank plots show which parameter is the problem:

fit_jax <- hb_area(y ~ x1 + x2 + x3,
  data = mys, vardir = "vardir",
  engine = jax_engine(chains = 4, samples = 1000)
)
fit_jax$rhat_max # below 1.01 is recommended
fit_jax$ess_min # above 100 per chain
head(fit_jax$mcmc_diagnostics) # R-hat and ESS of each parameter
autoplot(fit_jax, type = "rank") # chains near the dashed line: mixed
autoplot(fit_jax, type = "trace", pars = c("sigma2_u", "hb[5]"))

Spatial and Spatio-Temporal Structures

Spatial Priors

BYM2 (spatial = "bym2"): The scaled Besag-York-Mollié model (Riebler et al., 2016) 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 and has the penalised complexity (PC) prior $P(\phi < 0.5) = 0.5$ (Simpson et al., 2017).
Besag ICAR (spatial = "besag"): Intrinsic Conditional Autoregressive spatial model.
SAR area effect (spatial = "sar"): $u = (I - \rho W)^{-1}\varepsilon$ , the random effect of the spatial Fay-Herriot model of eblup_sfh().
Spatial Lag Model (spatial = "slm"): $\eta = (I - \rho W)^{-1}(X\beta + \varepsilon)$ ; the covariates are filtered too, so this is not the model of eblup_sfh().

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": Knorr-Held (2000) interactions δdt\delta_{dt} added to the main effects, ηdt=xdt⊤β+sd+γt+δdt\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").

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

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.4434

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:

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.074072 
#> sigma2_t (temporal): 0.008529 
#> 
#> Coefficients:
#>                    beta   std.error      zvalue      pvalue    ci_lower
#> (Intercept) -2.2696e+00  1.5370e-02 -1.4766e+02  0.0000e+00 -2.2998e+00
#> x            4.1691e-01  6.6161e-02  6.3015e+00  2.9474e-10  2.8685e-01
#>             ci_upper
#> (Intercept)  -2.2395
#> x             0.5465
#> 
#> Hyperparameters:
#>                                  mean        sd 0.025quant  0.5quant 0.975quant
#> Precision for ..domain_id..  15.25261  5.508501   7.224452  14.33819   28.60609
#> Precision for ..time_id..   136.59559 57.256988  60.484698 125.04499  281.20961
#>                                  mode
#> Precision for ..domain_id..  12.67086
#> Precision for ..time_id..   104.92249
#> 
#> Goodness of Fit:
#>                                                   DIC 
#>                                            -971.94527 
#>                                                    pD 
#>                                              60.69144 
#>                                                  WAIC 
#>                                            -979.78124 
#>                                                 pWAIC 
#>                                              42.46184 
#> Marginal_LogLik.log marginal-likelihood (integration) 
#>                                             432.98608 
#> 
#> HB Summary Statistics:
#>        hb           linear_pred           sd                mse           
#>  Min.   :0.04905   Min.   :-2.970   Min.   :0.004907   Min.   :2.408e-05  
#>  1st Qu.:0.08142   1st Qu.:-2.428   1st Qu.:0.007907   1st Qu.:6.253e-05  
#>  Median :0.09444   Median :-2.266   Median :0.009723   Median :9.454e-05  
#>  Mean   :0.09673   Mean   :-2.268   Mean   :0.009911   Mean   :1.057e-04  
#>  3rd Qu.:0.10925   3rd Qu.:-2.103   3rd Qu.:0.011574   3rd Qu.:1.339e-04  
#>  Max.   :0.15774   Max.   :-1.678   Max.   :0.018014   Max.   :3.245e-04  
#>       rse        
#>  Min.   : 6.545  
#>  1st Qu.: 9.167  
#>  Median :10.224  
#>  Mean   :10.312  
#>  3rd Qu.:11.594  
#>  Max.   :14.895

3. Inspect Posterior Predictions and Uncertainty

The 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.05102819 0.005596831 10.968116 0.04062428
#> 2   CITTA' DI BOLOGNA 2014 0.0825 0.08182495 0.006281323  7.676537 0.07003566
#> 3               IMOLA 2014 0.1033 0.09310372 0.010890484 11.697152 0.07366705
#> 4         PIANURA EST 2014 0.0633 0.06365559 0.006408099 10.066829 0.05185446
#> 5       PIANURA OVEST 2014 0.0625 0.06937249 0.007927315 11.427174 0.05493782
#> 6      PORRETTA TERME 2014 0.1276 0.09446596 0.014070599 14.894889 0.06968708
#>     ci_upper
#> 1 0.06256385
#> 2 0.09467494
#> 3 0.11635694
#> 4 0.07700995
#> 5 0.08602893
#> 6 0.12479142

Comparison and Validation Against `tipsae` (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	22 min (4 chains x 10,000 iterations, largest $\hat R$ 1.007); 4 chains x 2,000 iterations (120 s) have not converged	1.1 s
Posterior means vs converged `tipsae`	Reference	$r = 0.999$ , mean absolute difference 0.0017

Measured on 2026-10-16 (tipsae reference refitted on 2026-10-17) 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.

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:

Execution time (median, seconds) — Standard Beta · log scale

Speedup Factor & Efficiency Multiplier (fastsae vs tipsae) — Standard Beta

Agreement with Stan MCMC (Beta model, $n = 1,000$ )

To check that the speed of fastsae::hb_area does not come at the cost of its estimates, both packages were fitted with their defaults to the $n = 1,000$ data set of the Beta benchmark (independent area effects, sampling variances from the direct estimates; benchmarks/bench_beta_concordance.R). The two packages do not fit identical models (their priors differ, and fastsae fixes the precision at a first-stage estimate while tipsae makes it a function of the model mean), so the comparison is between what a user obtains from each package:

Metric	`fastsae::hb_area` (INLA)	`tipsae::fit_sae` (Stan)
RMSE against the true area means	0.0678	0.0684
Coverage of the 95% intervals	93.3%	91.1%
Mean width of the 95% intervals	0.252	0.259
Time	4.3 s	184 s

The posterior means of the two packages differed by 0.0096 on average (absolute difference), and the posterior SDs of fastsae were 0.99 times those of tipsae in the median area (90% of the areas between 0.88 and 1.17). The largest $\hat{R}$ of the tipsae fit was 1.008.

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:

# 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.069376 
#> rho (spatial): 0.1595 
#> phi (spatial fraction): 0.1595 
#> 
#> Coefficients:
#>                    beta   std.error      zvalue      pvalue    ci_lower
#> (Intercept)  1.4161e-01  1.3846e-01  1.0228e+00  3.0643e-01 -1.3129e-01
#> x1           3.8646e-01  5.0100e-02  7.7137e+00  1.2220e-14  2.8774e-01
#> x2          -2.9875e-01  3.3665e-02 -8.8741e+00  7.0479e-19 -3.6547e-01
#>             ci_upper
#> (Intercept)   0.4149
#> x1            0.4853
#> x2           -0.2327
#> 
#> Hyperparameters:
#>                                   mean        sd  0.025quant   0.5quant
#> Precision for ..domain_id.. 15.6289832 4.4994272 8.488620837 15.0603778
#> Phi for ..domain_id..        0.1595215 0.1562642 0.007221511  0.1055182
#>                             0.975quant        mode
#> Precision for ..domain_id.. 26.0450162 14.01275025
#> Phi for ..domain_id..        0.5933893  0.01700406
#> 
#> Goodness of Fit:
#>                                                   DIC 
#>                                             317.74616 
#>                                                    pD 
#>                                              31.88700 
#>                                                  WAIC 
#>                                             309.67484 
#>                                                 pWAIC 
#>                                              17.15906 
#> Marginal_LogLik.log marginal-likelihood (integration) 
#>                                            -185.87988 
#> 
#> HB Summary Statistics:
#>        hb          linear_pred            sd               mse           
#>  Min.   :0.3061   Min.   :-1.1959   Min.   :0.02962   Min.   :0.0008776  
#>  1st Qu.:0.7648   1st Qu.:-0.2791   1st Qu.:0.05813   1st Qu.:0.0033800  
#>  Median :1.3850   Median : 0.3240   Median :0.07685   Median :0.0059062  
#>  Mean   :1.5294   Mean   : 0.2381   Mean   :0.16304   Mean   :0.0690516  
#>  3rd Qu.:2.3016   3rd Qu.: 0.8023   3rd Qu.:0.12189   3rd Qu.:0.0148681  
#>  Max.   :3.6016   Max.   : 1.2810   Max.   :0.76549   Max.   :0.5859683  
#>       rse        
#>  Min.   : 2.791  
#>  1st Qu.: 5.164  
#>  Median : 6.265  
#>  Mean   :10.478  
#>  3rd Qu.: 9.915  
#>  Max.   :30.137

Covariates Measured with Error

Area-level covariates are often estimates themselves, for instance a mean from another survey with its own sampling variance. Used as if it were exact, such a covariate gets an attenuated coefficient, every area is shrunk alike however noisy its covariate is, and the model-based estimate can be worse than the direct one (Ybarra and Lohr, 2008). var_x gives the variances of these errors: the true values of the covariates then replace the observed ones in the model, with a normal distribution (me_model = "structural", the default; Bell et al., 2019) or a flat prior (me_model = "functional", the Bayesian counterpart of the Ybarra–Lohr model; Arima, Datta and Liseo, 2015). A variance of 0 marks a value without error.

set.seed(2)
D <- 60
x_true <- rnorm(D, 5, 2)                              # the covariate, not observed
var_x1 <- sample(c(0.05, 0.5, 4), D, replace = TRUE)  # variances of its estimates
d_me <- data.frame(
  area = seq_len(D),
  x1 = x_true + rnorm(D, 0, sqrt(var_x1)),            # the estimates
  var_x1 = var_x1,
  vardir = runif(D, 0.5, 1.5)
)
d_me$y <- 1 + 0.8 * x_true + rnorm(D, 0, 0.7) + rnorm(D, 0, sqrt(d_me$vardir))
d_me$y[1:6] <- NA                                     # areas without a direct estimate

fit_naive <- hb_area(y ~ x1, data = d_me, domain = "area", vardir = "vardir", print_result = FALSE)
fit_me <- hb_area(y ~ x1,
  data = d_me, domain = "area", vardir = "vardir",
  var_x = c(x1 = "var_x1"), print_result = FALSE
)

# Coefficient of the covariate: 0.8 in the simulation
rbind(
  naive = fit_naive$estcoef["x1", c("beta", "ci_lower", "ci_upper")],
  structural = fit_me$estcoef["x1", c("beta", "ci_lower", "ci_upper")]
)
#>                 beta  ci_lower  ci_upper
#> naive      0.5855003 0.4338812 0.7374279
#> structural 0.7281791 0.5811399 0.8799180

# Mean posterior SD of the areas by the error variance of their covariate
sds <- data.frame(var_x1 = d_me$var_x1, naive = fit_naive$df_hb$sd, structural = fit_me$df_hb$sd)
aggregate(cbind(naive, structural) ~ var_x1, data = sds, FUN = mean)
#>   var_x1     naive structural
#> 1   0.05 0.8049357  0.6120351
#> 2   0.50 0.7694047  0.6623815
#> 3   4.00 0.7674869  0.7990772

The naive interval for the coefficient misses its true value, and the naive model gives areas with precise and with noisy covariates about the same posterior SD. The measurement error model narrows the intervals of the areas whose covariate is precise and widens those whose covariate is noisy. With me_model = "functional" the area estimates agree with the Ybarra–Lohr EBLUP (emdi::fh(method = "me")), but the coefficient stays close to the naive one: integrating the true values out leaves the regression on the observed values, only with a larger variance. The details, and the validation against Stan, are in the section Covariates measured with error of ?hb_area.

References

Arima, S., Datta, G. S., & Liseo, B. (2015). Bayesian estimators for small area models when auxiliary information is measured with error. Scandinavian Journal of Statistics, 42(2), 518–529.
Bell, W. R., Chung, H. C., Datta, G. S., & Franco, C. (2019). Measurement error in small area estimation: Functional versus structural versus naive models. Survey Methodology, 45(1), 61–80.
De Nicolò, S., & Gardini, A. (2024). The R package tipsae: Tools for mapping proportions and indicators on the unit interval. Journal of Statistical Software, 108(1), 1–36.
Janicki, R. (2020). Properties of the beta regression model for small area estimation of proportions and application to estimation of poverty rates. Communications in Statistics - Theory and Methods, 49(9), 2264–2284.
Knorr-Held, L. (2000). Bayesian modelling of inhomogeneous spatial and temporal variation in rates. Statistics in Medicine, 19(17-18), 2555–2567.
Riebler, A., Sørbye, S. H., Simpson, D., & Rue, H. (2016). An intuitive Bayesian spatial model for disease mapping that accounts for scaling. Statistical Methods in Medical Research, 25(4), 1145–1165.
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.
Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C. (2021). Rank-normalization, folding, and localization: An improved R-hat for assessing convergence of MCMC. Bayesian Analysis, 16(2), 667–718.
Ybarra, L. M. R., & Lohr, S. L. (2008). Small area estimation when auxiliary information is measured with error. Biometrika, 95(4), 919–931.