Package 'hbsaems'

Province-level Adjacency Matrix

Description

A small example adjacency matrix used for fitting Conditional Autoregressive (CAR) random effects on the **province** level. Pairs with data_fhnorm and data_betalogitnorm, whose province column has matching labels.

Usage

adjacency_matrix_car

Format

A binary symmetric $5 \times 5$ matrix with row- and column-names province_01 .. province_05; entries are 1 for adjacent province pairs and 0 otherwise.

Source

Simulated.

---

Regency-level Adjacency Matrix (Coarse Spatial Cluster)

Description

A small example adjacency matrix used for fitting Conditional Autoregressive (CAR) random effects on the **regency** (coarse spatial-cluster, kabupaten) level. Pairs with data_binlogitnorm and data_lnln, whose regency column has matching labels.

Usage

adjacency_matrix_car_regency

Format

A binary symmetric $5 \times 5$ matrix with row- and column-names regency_01 .. regency_05; entries are 1 for adjacent regency pairs and 0 otherwise.

Note

The naming regency_01..05 (two-digit suffix) is reserved in this package for the COARSE 5-level cluster used by data_binlogitnorm and data_lnln. The 100-level FINE regency column in data_fhnorm and data_betalogitnorm uses three-digit suffixes (regency_001..100) and pairs with the larger spatial_weight_sar matrix. See vignette("hbsaems-spatial") for the naming convention.

Source

Simulated.

---

Loglogistic as a Custom Distribution Family for brms

Description

Returns the custom_family + stanvars pair required to fit a brms model with a Loglogistic response distribution. The Stan code is loaded from inst/stan/loglogistic.stan via build_brms_custom_family; post-processing functions (log_lik, posterior_predict, posterior_epred) are wired in automatically.

Usage

brms_custom_loglogistic()

Details

Two parameters:

mu

Scale (mu > 0, log link).

beta

Shape (beta > 0, log link).

Value

A list with elements custom_family and stanvars_family ready for use with brms::brm(). The returned custom_family object has log_lik, posterior_predict, and posterior_epred registered so that brms::loo(), brms::posterior_predict(), and brms::posterior_epred() work without further setup.

See Also

dloglogistic, build_brms_custom_family, register_hbsae_brms_custom.

Examples

library(hbsaems)
library(brms)
ll <- brms_custom_loglogistic()
# fit <- brm(y ~ x, data = d,
#            family   = ll$custom_family,
#            stanvars = ll$stanvars_family)
# loo(fit)                     # uses log_lik_loglogistic
# posterior_predict(fit)       # uses posterior_predict_loglogistic
# posterior_epred(fit)         # uses posterior_epred_loglogistic

---

Shifted Loglogistic as a Custom Distribution Family for brms

Description

Returns the custom_family + stanvars pair required to fit a brms model with a Shifted Loglogistic response distribution. The Stan code is loaded from inst/stan/shifted_loglogistic.stan via build_brms_custom_family; post-processing functions (log_lik, posterior_predict, posterior_epred) are wired in automatically.

Usage

brms_custom_shifted_loglogistic()

Details

Three parameters:

mu

Location parameter (identity link).

sigma

Scale parameter (sigma > 0, log link).

xi

Shape parameter (identity link).

Value

A list with elements custom_family and stanvars_family ready for use with brms::brm().

See Also

dshifted_loglogistic, build_brms_custom_family, register_hbsae_brms_custom.

Examples

library(hbsaems)
library(brms)
sll <- brms_custom_shifted_loglogistic()
# fit <- brm(y ~ x, data = d,
#            family   = sll$custom_family,
#            stanvars = sll$stanvars_family)

---

Build a brms Custom Family + Stanvars Pair from a Single Spec

Description

Convenience function that combines custom_family and stanvar into the standard list(custom_family, stanvars_family) pair returned by every brms_custom_* helper in hbsaems. The Stan code is read directly from inst/stan/<name>.stan – the user does not need to maintain it as an R string.

Usage

build_brms_custom_family(
  name,
  dpars,
  links,
  lb = NA,
  ub = NA,
  type = c("real", "int"),
  loop = FALSE,
  log_lik = NULL,
  posterior_predict = NULL,
  posterior_epred = NULL,
  use_stan_native = FALSE
)

Arguments

name	Character. Family name; must match a <name>.stan file under inst/stan/.
dpars	Character vector of distributional parameter names. The first MUST be "mu" (brms convention).
links	Character vector of link functions, same length as dpars. Common values: "identity", "log", "logit".
lb, ub	Numeric vectors of lower / upper bounds; NA for none.
type	"real" (continuous) or "int" (discrete).
loop	Logical. FALSE (default) selects the vectorised brms convention (Stan signatures take vectors of y and mu); TRUE selects scalar Stan signatures. The default matches the neodistr convention. Whichever you choose, the corresponding .stan file under inst/stan/ must use the same convention.
log_lik	Optional function for computing observation-level log-likelihoods (used by loo(), waic()). Signature must be function(i, prep). See the brms vignette "Define Custom Response Distributions" for details.
posterior_predict	Optional function for drawing from the posterior predictive distribution (used by pp_check(), posterior_predict()). Signature function(i, prep, ...).
posterior_epred	Optional function returning the conditional expectation $E[Y \mid X]$ (used by fitted(), conditional_effects()). Signature function(prep).
use_stan_native	Logical. If TRUE, no stanvars block is attached – the <name>_lpdf function is assumed to exist as a Stan built-in (e.g.\ loglogistic_lpdf in Stan $\geq$ 2.29). When FALSE (the default), the function definition is loaded from inst/stan/<name>.stan.

Details

Stan code is built once per call; if your code is reused many times, wrap the returned object in a function and call it lazily.

Value

A list with two elements:

custom_family

A brms::customfamily object.

stanvars_family

A brms::stanvars object with the Stan code in the functions block, or NULL when use_stan_native = TRUE.

See Also

read_stan_function, register_hbsae_brms_custom.

Examples

library(hbsaems)
library(brms)

# Build the loglogistic family.  Stan code is loaded from
# inst/stan/hbsae_loglogistic.stan — the `hbsae_` prefix avoids
# collision with Stan's BUILT-IN `loglogistic_lpdf` (Stan >= 2.29).
ll <- build_brms_custom_family(
  name             = "hbsae_loglogistic",
  dpars            = c("mu", "beta"),
  links            = c("log", "log"),
  lb               = c(0,    0),
  ub               = c(NA,   NA),
  type             = "real"
)
class(ll$custom_family)

---

build_spatial_weight: Construct M for CAR / SAR models

Description

Constructs a spatial weight / adjacency matrix $M$ suitable for use as the M argument of hbm. The function accepts either a path to a shapefile (.shp) or an sf / Spatial* object already in memory, and returns a square numeric matrix that is automatically validated against the theoretical requirements of the target model class via check_spatial_weight.

Usage

build_spatial_weight(
  shp,
  for_model = NULL,
  type = NULL,
  style = NULL,
  k = 4L,
  threshold = NULL,
  id_col = NULL,
  longlat = FALSE,
  validate = TRUE
)

Arguments

shp	Either a character path to a .shp file, an sf object, or a Spatial* (sp) object.
for_model	Optional convenience argument: "car" or "sar". When supplied, sets sensible defaults for type and style as documented above. When the user also supplies type or style explicitly, those take precedence.
type	Character. Neighbour definition: "queen" Polygons sharing at least one vertex. Symmetric. "rook" Polygons sharing an edge. Symmetric. "knn" k nearest centroids. Asymmetric in general; auto-symmetrised when style = "B". "distance" Centroids within threshold. Symmetric.
style	Character. "B" (binary, suitable for CAR) or "W" (row-standardised, suitable for SAR).
k	Integer. Number of nearest neighbours when type = "knn" (default 4).
threshold	Numeric. Distance threshold (in CRS units) when type = "distance". When NULL, the maximum first-nearest-neighbour distance is used.
id_col	Optional character. Column name in the spatial object whose values become the matrix dimnames.
longlat	Logical. Whether coordinates are in longitude/latitude.
validate	Logical. Whether to run check_spatial_weight after construction (default TRUE). Set to FALSE to skip diagnostics.

Details

Build a Spatial Weight or Adjacency Matrix from Shapefile

Choosing type and style for your SAE model

The combination of type (neighbour definition) and style (matrix coding) determines whether the resulting matrix is theoretically appropriate for a CAR or SAR model. The recommended combinations are:

Model	type	style	Reference
CAR / ICAR / BYM2	queen or rook	B	Besag (1974, 1991), Riebler et al. (2016)
SAR (lag/error)	knn or distance	W	Anselin (1988), Whittle (1954)
CAR (irregular)	knn or distance	B	kNN auto-symmetrised when style="B"

Use the convenience wrapper for_model = "car" or for_model = "sar" to set both type and style automatically from the recommendation above. When type or style is also supplied, the explicit argument wins (the helper only fills in the missing one).

Mathematical contracts

CAR

$M$ must be binary, symmetric, and have zero diagonal. The full distribution is $u \sim \mathcal{N}\bigl(0, \sigma^2 (D - \rho W)^{-1}\bigr)$ with $D = \mathrm{diag}(W \mathbf{1})$.

SAR

$M$ should be row-standardised so that $\rho$ lies in $(-1, 1)$. Symmetry is not required.

All matrices are post-checked with check_spatial_weight and offending matrices are flagged.

Value

A square numeric matrix with row/column names taken from id_col (if supplied) or sequential integers, plus the following attributes: "hbsaems_type", "hbsaems_style", "hbsaems_for_model", "hbsaems_check" (the result of check_spatial_weight).

Reproducibility note

KNN graphs depend on the order in which ties are broken, so when several centroids are exactly equidistant the resulting matrix may depend on feature ordering.

References

Anselin, L. (1988). Spatial Econometrics: Methods and Models. Kluwer Academic Publishers.

Besag, J. (1974). Spatial interaction and the statistical analysis of lattice systems (with discussion). JRSS B 36(2), 192–236.

Bivand, R. S., Pebesma, E., & Gomez-Rubio, V. (2013). Applied Spatial Data Analysis with R (2nd ed.). Springer.

See Also

check_spatial_weight, hbm, adjacency_matrix_car, spatial_weight_sar

Examples

# Inspect a pre-built CAR adjacency matrix shipped with the package
data("adjacency_matrix_car")
dim(adjacency_matrix_car)
check_spatial_weight(adjacency_matrix_car, spatial_model = "car",
                      verbose = FALSE)$compatible


# Build a CAR matrix from an sf object (requires sf + spdep)
if (requireNamespace("sf", quietly = TRUE) &&
    requireNamespace("spdep", quietly = TRUE)) {
  library(sf)
  # A small 2x2 grid of polygons
  g <- st_sf(
    id = LETTERS[1:4],
    geometry = st_sfc(
      st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0)))),
      st_polygon(list(rbind(c(1,0), c(2,0), c(2,1), c(1,1), c(1,0)))),
      st_polygon(list(rbind(c(0,1), c(1,1), c(1,2), c(0,2), c(0,1)))),
      st_polygon(list(rbind(c(1,1), c(2,1), c(2,2), c(1,2), c(1,1))))
    )
  )
  M_car <- build_spatial_weight(g, for_model = "car")
  M_sar <- build_spatial_weight(g, for_model = "sar", k = 2L)
}

---

Inspect Data Before Fitting an HBSAE Model

Description

Performs three independent checks on the supplied dataset and returns a structured hbsaems_data_check object that summarises the results. This function is intended to be called before hbm or any of the distribution-specific wrappers.

Usage

check_data(
  data,
  response,
  auxiliary = NULL,
  area_var = NULL,
  spatial_var = NULL,
  M = NULL,
  trials = NULL,
  n_var = NULL,
  predictors = NULL,
  group = NULL,
  sre = NULL
)

Arguments

data	A data.frame.
response	Character. Name of the response variable in data.
auxiliary	Character vector. Names of the auxiliary variables (the SAE 'X' covariates).
area_var	Optional character. Name of the area (random-effect grouping) variable.
spatial_var	Optional character. Name of the spatial-grouping variable (column in data that indexes rows of the spatial weight matrix M).
M	Optional spatial weight matrix to dimension-check against data.
trials	Optional character. Name of the trials variable (binomial models).
n_var	Optional character. Name of the sample-size variable (beta / lognormal direct-estimator models).
predictors	Deprecated. Use auxiliary instead.
group	Deprecated. Use area_var instead.
sre	Deprecated. Use spatial_var instead.

Details

1. Variable presence

Verifies that response, every name in auxiliary, and the optional area_var / spatial_var / trials / n_var columns exist in data. Missing variables are reported in $issues.

2. Missing-value pattern

The pattern is one of:

"none"

All listed columns are complete.

"y_only"

Only the response has NAs.

"x_only"

Only the auxiliary variables have NAs.

"both"

Both Y and X have NAs.

Based on the pattern, a strategy is recommended:

"y_only"

First, check whether the NA-Y rows are non-sample (out-of-sample) areas – domains for which a prediction is desired but no direct estimate exists. If yes, do not delete them; fit on the complete-Y subset and pass the NA-Y rows to sae_predict via the newdata argument. If they are merely missing observations within sampled areas, use handle_missing = "deleted".

"x_only"

handle_missing = "multiple" – multiple imputation via mice.

"both" (continuous outcome)

handle_missing = "model" – joint Bayesian imputation via brms::mi().

"both" (discrete outcome, binomial)

handle_missing = "multiple".

3. Dimension check

When M is supplied, verifies that it is square and that nrow(M) matches the number of distinct levels in data[[spatial_var]] (or nrow(data) when spatial_var is NULL).

Value

An object of class hbsaems_data_check with components:

n_obs

Number of rows in the data.

missing_summary

Named integer vector: per-variable count of NA.

missing_pattern

Character: "none", "y_only", "x_only", or "both".

dimension_check

Named list of dimension diagnostics.

non_sample_warning

Character or NULL – a hint to investigate whether NA-Y rows are non-sample (out-of-sample) areas.

recommended_method

Character: suggested handle_missing value ("deleted", "multiple", "model", or NA when no missing values are present).

recommendation_text

Human-readable rationale.

issues

Character vector of fatal errors (length 0 if OK).

See Also

hbm, sae_predict

Examples

data("data_fhnorm")

# 1. Complete data -> no warnings, no recommendation
chk <- check_data(data_fhnorm,
                  response   = "y",
                  auxiliary  = c("x1", "x2", "x3"))
print(chk)

# 2. Missing-Y pattern -> recommends checking for non-sample areas
d <- data_fhnorm
d$y[1:5] <- NA
chk2 <- check_data(d, response = "y",
                      auxiliary  = c("x1", "x2", "x3"))
summary(chk2)

# 3. Missing-X-only -> recommends multiple imputation
d2 <- data_fhnorm
d2$x1[10:15] <- NA
chk3 <- check_data(d2, response = "y",
                      auxiliary  = c("x1", "x2", "x3"))
chk3$recommended_method

# 4. Spatial dimension check
data("adjacency_matrix_car")
chk4 <- check_data(data_fhnorm[1:5, ],
                   response   = "y",
                   auxiliary  = c("x1", "x2", "x3"),
                   M          = adjacency_matrix_car)
chk4$dimension_check

---

Check Shiny App Dependencies

Description

Classifies and inspects the optional packages used by the Shiny dashboard. The dashboard is launched by run_sae_app and the dependencies it touches live in Suggests, not Imports, so users who only use the modelling functions never have to install them.

Usage

check_shiny_deps(verbose = TRUE)

Arguments

verbose

Logical. When TRUE (default), prints a formatted summary of which dependencies are installed and which are missing.

Value

Invisibly, a named list with components:

critical_missing

Character vector of critical packages not installed. When non-empty, the app cannot start.

optional_missing

Character vector of optional packages not installed. When non-empty, the app starts but some panels degrade.

install_cmd

A ready-to-paste install.packages() call covering all missing packages, or NULL when none are missing.

See Also

run_sae_app

Examples

check_shiny_deps()

---

Validate a Spatial Weight Matrix Against CAR/SAR Theory

Description

Runs five theoretical compatibility checks on a spatial weight matrix $M$ and reports any deviations from the standard requirements of the chosen model class. Returns a structured object summarising the results.

Usage

check_spatial_weight(
  M,
  spatial_model = c("car", "sar"),
  verbose = TRUE,
  sre_type = NULL
)

Arguments

M	A square numeric matrix.
spatial_model	Character. "car" or "sar" – the model class the matrix is intended for.
verbose	Logical. When TRUE (default), prints a formatted diagnostic report.
sre_type	Deprecated. Use spatial_model instead. Kept for backward compatibility; will be removed in v2.0.0.

Details

Theoretical requirements

For a CAR model (Besag 1974), the joint distribution

$u \sim \mathcal{N}\bigl(0,\, \sigma^2 (D - \rho W)^{-1}\bigr)$

is well-defined only when $W$ is symmetric with zero diagonal. By convention, $W$ is taken as the binary adjacency matrix (style = "B") so that $D = \mathrm{diag}(\text{row sums})$ has integer entries.

For a SAR model (Whittle 1954, Anselin 1988),

$u = \rho W u + \varepsilon, \quad \varepsilon \sim \mathcal{N}(0, \sigma^2 I)$

and $W$ is conventionally row-standardised (style = "W") so that the spatial autoregressive parameter $\rho$ can be interpreted as a normalised correlation in $(-1, 1)$. Symmetry is not required for SAR.

Style detection heuristic

• "B" (binary): all values in {0, 1}.

• "W" (row-standardised): every non-zero row sums to 1.

• "other": neither of the above.

Value

Invisibly, an object of class hbsaems_spatial_check with components:

is_square

Logical.

has_zero_diag

Logical.

is_symmetric

Logical.

detected_style

Character: "B", "W", or "other".

n_isolated

Integer: number of areas with no neighbours.

n_components

Integer: number of connected components.

issues

Character vector of fatal errors.

warnings

Character vector of soft warnings.

compatible

Logical: TRUE if matrix is theoretically compatible with spatial_model.

See Also

build_spatial_weight, hbm

Examples

# Build a small valid CAR matrix
M <- matrix(c(0, 1, 1, 0,
              1, 0, 0, 1,
              1, 0, 0, 1,
              0, 1, 1, 0), 4, 4)
check_spatial_weight(M, spatial_model = "car")

# An asymmetric matrix flagged for CAR
M2 <- M; M2[1, 2] <- 2
check_spatial_weight(M2, spatial_model = "car", verbose = FALSE)$issues

---

MCMC Convergence Diagnostics for Fitted HBMs

Description

Computes a battery of convergence tests and diagnostic plots for an hbmfit object. This is the primary convergence diagnostic function in hbsaems (supersedes the deprecated hbcc).

Usage

convergence_check(
  model,
  diag_tests = c("rhat", "geweke", "heidel", "raftery"),
  plot_types = c("trace", "dens", "acf", "nuts_energy", "rhat", "neff"),
  ...
)

Arguments

model	An hbmfit or brmsfit object.
diag_tests	Character vector of tests to run. Any subset of c("rhat", "geweke", "heidel", "raftery").
plot_types	Character vector of plot types to generate. Any subset of c("trace", "dens", "acf", "nuts_energy", "rhat", "neff").
...	Additional arguments passed to the underlying brms or coda routines.

Details

For each parameter $\theta_j$, the Gelman-Rubin statistic is computed as

$\widehat{R}_j = \sqrt{\frac{\widehat{\mathrm{var}}^+(\theta_j)}{W_j}}$

where $W_j$ is the within-chain variance and $\widehat{\mathrm{var}}^+$ is the marginal posterior variance estimate. Values close to 1 (typically below 1.1) indicate convergence.

Value

An hbcc_results object containing:

rhat_ess

Matrix with columns Rhat, Bulk_ESS, Tail_ESS.

geweke,raftery,heidel

Outputs from the corresponding coda routines, or NULL when the test fails.

plots

Named list of ggplot/bayesplot objects.

See Also

is_converged, diagnostic_summary, hbm_warnings

Examples

library(hbsaems)
library(brms)
data("data_fhnorm")
# Uses brms's default MCMC settings (chains = 4, iter = 2000,
# warmup = 1000).  For challenging posteriors (e.g. funnel
# geometries in Fay-Herriot with small D_i), consider
# chains = 4, iter = 4000, warmup = 2000 and
# control = list(adapt_delta = 0.99).
model <- hbm(brms::bf(y ~ x1 + x2 + x3),
             data   = data_fhnorm,
             re     = ~ (1 | regency),    # area-level random effect
             chains = 4, iter = 2000, warmup = 1000,
             cores  = 1, seed = 123, refresh = 0)

diag <- convergence_check(model)
summary(diag)
is_converged(model)
is_converged(model, threshold = 1.05)

---

Simulated Beta Logit-Normal Data

Description

A simulated dataset for 100 regencies under a Beta logit-normal model. The response y is a proportion in $(0, 1)$.

Usage

data_betalogitnorm

Format

A data frame with 100 rows and 9 variables:

y

Direct estimator of the area-level proportion.

theta

True area-level proportion.

x1, x2, x3

Auxiliary covariates.

n

Area sample size.

deff

Design effect.

regency

Regency identifier ("regency_001" .. "regency_100").

province

Province identifier ("province_01" .. "province_05") – spatial cluster level for CAR / SAR.

Source

Simulated.

---

Simulated Binomial Logit-Normal Data

Description

A simulated dataset for 100 districts under a Binomial logit-normal model. Each district provides a number of successes (y) out of n trials.

Usage

data_binlogitnorm

Format

A data frame with 100 rows and 14 variables:

n

Number of trials in the district.

y

Number of successes.

p

Direct proportion (y / n).

x1, x2, x3

Auxiliary covariates.

u_true

True area-level random effect (logit scale).

eta_true

True linear predictor (logit scale).

p_true

True success probability.

psi_i

Sampling variance.

y_obs, p_obs

Observed (direct) values.

district

District identifier ("district_001" .. "district_100") – the 100 small areas to estimate.

regency

Regency identifier ("regency_01" .. "regency_05") – coarse spatial-cluster level (5 regencies each containing 20 districts). Pair with adjacency_matrix_car_regency.

Source

Simulated.

---

Simulated Fay-Herriot Normal Data

Description

A simulated dataset for 100 regencies under the Fay-Herriot Normal small-area model. Used as the running example throughout the package documentation and vignettes. The simulation is engineered so that the canonical Fay-Herriot fit (hbm(..., sampling_variance = "D")) converges with default brms / Stan settings – no divergent transitions, no manual tuning required.

Usage

data_fhnorm

Format

A data frame with 100 rows and 9 variables:

y

Direct (survey) estimator of the area mean.

D

Sampling variance of the direct estimator (KNOWN from the survey design; treat as input, not as a parameter).

x1, x2, x3

Auxiliary covariates at the area level, simulated from $\mathcal{N}(0, 1)$.

theta_true

True area-level latent value $\theta_i$.

u

True area-level random effect $u_i$.

regency

Regency identifier ("regency_001" through "regency_100") used as the IID random-effect grouping variable. Use with re = ~ (1 | regency) or area_var = "regency".

province

Province identifier ("province_01" through "province_05") – 20 regencies per province. Used as the spatial random-effect grouping variable for CAR / SAR / BYM2 examples; also serves as the higher level in the hierarchical-area example area_var = c("province", "regency").

Details

Generative model. For each regency $i = 1, \ldots, 100$,

$y_i = \theta_i + \varepsilon_i, \quad \varepsilon_i \sim \mathcal{N}(0, D_i)$

$\theta_i = 10 + 0.8 \, x_{1i} - 0.5 \, x_{2i} + 0.3 \, x_{3i} + u_i, \quad u_i \sim \mathcal{N}(0, \sigma_u^2)$

with auxiliary covariates $x_j \sim \mathcal{N}(0, 1)$ (already standardised), area RE SD $\sigma_u = 1$, and known sampling variances $D_i \sim \mathrm{Gamma}(\mathrm{shape} = 4, \mathrm{rate} = 4)$ – a realistic spread ($\approx [0.2, 3.0]$) that mirrors varying sample sizes across regencies.

Important: pass D as the sampling variance. In any fit on this dataset, supply sampling_variance = "D"; otherwise the residual $\sigma$ and the area-RE $\sigma_u$ compete to explain the same variance, producing weak identifiability and divergent transitions.

Source

Simulated. Reproducible script in data-raw/data_fhnorm.R.

---

Simulated Lognormal-Lognormal Data

Description

A simulated dataset for 100 districts under a Lognormal-Lognormal model. Suitable for strictly positive, right-skewed outcomes.

Usage

data_lnln

Format

A data frame with 100 rows and 13 variables:

district

District identifier ("district_001" .. "district_100") – the 100 small areas to estimate.

x1, x2, x3

Auxiliary covariates.

u_true

True area-level random effect (log scale).

teta_true

True linear predictor (log scale).

mu_orig_true

True mean on the original scale.

n

Area sample size.

y_obs

Observed direct estimator.

lambda_dir

Direct estimator scale parameter.

y_log_obs

Observed value on the log scale.

psi_i

Sampling variance on the log scale.

regency

Regency identifier ("regency_01" .. "regency_05") – coarse spatial-cluster level (5 regencies each containing 20 districts). Pair with adjacency_matrix_car_regency.

Source

Simulated.

---

Deprecated Functions

Description

These functions were the main entry point in hbsaems $\le$ 0.2.x. They are retained for backwards compatibility but now call the new primary functions and emit a deprecation warning via .Deprecated. They will be removed in v2.0.0.

Usage

hbcc(
  model,
  diag_tests = c("rhat", "geweke", "heidel", "raftery"),
  plot_types = c("trace", "dens", "acf", "nuts_energy", "rhat", "neff"),
  ...
)

hbmc(
  model,
  model2 = NULL,
  ndraws_ppc = 100,
  moment_match = FALSE,
  moment_match_args = list(),
  reloo_args = list(),
  plot_types = c("pp_check", "params"),
  comparison_metrics = c("loo", "waic", "bf"),
  run_prior_sensitivity = FALSE,
  sensitivity_vars = NULL,
  ...
)

hbpc(model, data = NULL, response_var = NULL, ndraws_ppc = 50, ...)

hbsae(model, newdata = NULL, ...)

Arguments

model	An hbmfit object.
diag_tests, plot_types, ndraws_ppc, moment_match, moment_match_args, reloo_args, comparison_metrics, run_prior_sensitivity, sensitivity_vars	Forwarded unchanged to the replacement function.
...	Additional arguments forwarded to the replacement.
model2	For hbmc: optional second model.
data	For hbpc: the data data.frame.
response_var	For hbpc: name of the response variable.
newdata	For hbsae: optional new data.

Value

Identical to the corresponding replacement function.

Migration guide

Deprecated	Replacement
hbcc(model, ...)	convergence_check(model, ...)
hbmc(model, ...)	model_compare(model, ...)
hbpc(model, data, response_var)	prior_check(model, data, response_var)
hbsae(model, ...)	sae_predict(model, ...)

---

Extract a Diagnostic Summary

Description

Extract a Diagnostic Summary

Usage

diagnostic_summary(x)

Arguments

x	An hbmfit or hbcc_results object.

Value

A named list of diagnostic statistics.

---

Inspect a Registered HBSAE Model Specification

Description

Inspect a Registered HBSAE Model Specification

Usage

get_hbsae_model(key)

Arguments

key	Character. A model key returned by list_hbsae_models.

Value

The named list spec, or NULL if not found. Useful fields include:

• family – brms family name passed to brmsfamily.

• link – default link function used by the family ("identity", "logit", "log", ...). See brmsfamily for the complete set of supported links per family.

• discrete – whether the response is discrete.

• supports_mi – whether brms-canonical mi() imputation is allowed for this family (FALSE for all discrete responses).

See Also

register_hbsae_model, brmsfamily for the canonical brms family / link reference.

Examples

get_hbsae_model("lognormal")
get_hbsae_model("beta")$link  # "logit"

---

hbm: Hierarchical Bayesian Small Area Models

Description

Fits a hierarchical Bayesian model for Small Area Estimation (SAE) using the brms package (Stan back-end). The function supports fixed effects, random effects, spatial random effects (CAR/SAR), user-defined priors, and three strategies for handling missing data in auxiliary (predictor) variables.

Usage

hbm(
  formula,
  family = NULL,
  hb_sampling = "gaussian",
  hb_link = "identity",
  link_phi = "log",
  re = NULL,
  spatial_var = NULL,
  spatial_model = NULL,
  car_type = NULL,
  sar_type = NULL,
  M = NULL,
  data,
  prior = NULL,
  fixed_params = NULL,
  sampling_variance = NULL,
  prior_type = "default",
  hs_df = 1,
  hs_df_global = 1,
  hs_df_slab = 4,
  hs_scale_global = NULL,
  hs_scale_slab = 2,
  hs_par_ratio = NULL,
  hs_autoscale = TRUE,
  r2d2_mean_R2 = 0.5,
  r2d2_prec_R2 = 2,
  r2d2_cons_D2 = NULL,
  r2d2_autoscale = TRUE,
  nonlinear = NULL,
  nonlinear_type = "spline",
  spline_k = -1L,
  spline_bs = "tp",
  gp_k = NA_integer_,
  gp_cov = "exp_quad",
  gp_c = NULL,
  gp_scale = NULL,
  handle_missing = NULL,
  m = 5L,
  mice_args = list(),
  measurement_error = NULL,
  control = list(),
  chains = 4L,
  iter = 4000L,
  warmup = floor(iter/2),
  cores = 1L,
  sample_prior = "no",
  sre = NULL,
  sre_type = NULL,
  stanvars = NULL,
  ...
)

Arguments

formula	A brmsformula or standard formula object specifying the model structure. For multi-response or imputation sub-models use bf. Examples: formula(y ~ x1 + x2), bf(y ~ x1 + x2), or bf(y | mi() ~ mi(x1) + x2) + bf(x1 | mi() ~ x2).
family	The model family (primary argument, brms-consistent). Accepts a family name string (e.g.\ "gaussian"), a brms family object (e.g.\ gaussian(), Gamma(link = "log")), or a registered custom family object. A link carried on the object is used when hb_link is left at its default. Supply either family or its alias hb_sampling, not both. family is the uniform way to choose the distribution across hbm() and hbm_flex; each function additionally keeps its historical alias (hb_sampling here, family_key in hbm_flex).
hb_sampling	Character string naming the distribution family of the response variable (default: "gaussian"). Any family supported by brmsfamily is accepted.
hb_link	Character string specifying the link function (default: "identity").
link_phi	Character string specifying the link function for the precision/phi parameter (default: "log"). Only used for Beta and related families.
re	An optional one-sided formula specifying group-level (random) effects, e.g. ~(1|area). Must follow standard lme4-style syntax: ~ (1|group1) + (1|group2). If NULL (default) and spatial_model is also NULL, the function emits a warning recommending an area-level random effect, since a Fay-Herriot SAE model with neither IID nor spatial area effects degenerates to fixed-effects regression and does not borrow strength across areas. The warning can be silenced with suppressWarnings() if a fixed-effects-only baseline is intentional.
spatial_var	Character. Name of the column in data that identifies the spatial areas (e.g. "regency" or "province"). Must be supplied together with spatial_model and M; providing only one of them is an error. Distinct from re: re is a formula for IID random effects, whereas spatial_var is a column name (string) for the spatially-structured random effect.
spatial_model	Character. Type of spatial model: "car" (Conditional Autoregressive) or "sar" (Simultaneous Autoregressive). Must be supplied together with spatial_var and M; providing only one of them is an error.
car_type	Character. CAR subtype passed to brms: one of "escar" (exact sparse CAR), "esicar" (exact sparse intrinsic CAR), "icar" (intrinsic CAR), or "bym2". Defaults to "icar" when spatial_model = "car".
sar_type	Character. SAR subtype: "lag" (SAR of the response values) or "error" (SAR of the residuals). Defaults to "lag" when spatial_model = "sar".
M	Spatial matrix supplied as data2 to brms. For CAR this must be a binary adjacency matrix; for SAR a spatial weight matrix. Row names must match the levels of spatial_var.
data	A data frame containing all variables referenced in formula.
prior	Priors specified via prior or a list thereof. If NULL (default), brms default priors are used.
fixed_params	Named list pinning distributional parameters to known values instead of sampling them. Each entry maps a parameter name to one of: A column name (character) The named column in data is used as the fixed values. A scalar (numeric, length 1) Broadcast to all rows. A vector (numeric, length nrow(data)) Used directly. A one-sided formula (e.g. ~ I(n / deff - 1)) Evaluated against data to produce the vector of fixed values. Internally, each pinned parameter <par> is attached to the data as a column .hbsaems_<par>_fixed and added to the brms formula as <par> ~ 0 + offset(.hbsaems_<par>_fixed). Using fixed_params on a parameter for which the user also supplies an explicit prior is an error. Typical use cases include pinning phi for Beta regression from survey design effect (phi = ~ I(n / deff - 1)) and pinning sigma for Fay–Herriot-style models with known sampling variance.
sampling_variance	Optional character. Name of a column in data containing the known sampling variance $D_i$ for each area (the Fay-Herriot sugar). When supplied, $\sigma_i = \sqrt{D_i}$ is pinned via offset. This is the canonical way to fit a Gaussian Fay-Herriot model: without it, the residual $\sigma$ and the area-RE $\sigma_u$ compete to explain the same variance and the model is unidentified, typically producing divergent transitions. Equivalent to fixed_params = list(sigma = sqrt(data[[<col>]])). Family compatibility. sampling_variance requires a continuous family whose response distribution has a residual scale parameter named sigma. Supported: gaussian (the canonical Fay-Herriot case), lognormal (D must be on the log scale; see also hbm_lnln), student, skew_normal, exgaussian, asym_laplace. A helpful error is thrown when an incompatible family is supplied. For non-Gaussian SAE families the analogous pinning mechanism is family-specific: Beta: pin the precision phi via the survey design effect, e.g.\ fixed_params = list(phi = ~ I(n/deff - 1)) (Liu 2009). See hbm_betalogitnorm. Binomial: sampling variability enters through the trials addition term, not through a separate sigma. See hbm_binlogitnorm. Poisson, Gamma, Weibull: variance is tied algebraically to the mean – no separate scale parameter to pin.
prior_type	Character. Global-local shrinkage prior applied to all regression coefficients (class = "b"). One of: "default" No shrinkage prior is added; the prior argument governs everything (default). "horseshoe" Regularised horseshoe prior (Piironen & Vehtari 2017). Encourages exact sparsity while allowing large signals through. Controlled by hs_df, hs_df_global, hs_df_slab, hs_scale_global, hs_scale_slab, hs_par_ratio, and hs_autoscale. "r2d2" R2D2 prior (Zhang et al. 2022). Places a prior directly on the model $R^2$ and distributes explained variance across predictors via a Dirichlet decomposition. Controlled by r2d2_mean_R2, r2d2_prec_R2, r2d2_cons_D2, and r2d2_autoscale. If prior already contains a global class = "b" entry, prior_type is ignored and a warning is issued. Cascading to smooth and GP terms. When the formula contains s() or gp() terms, the shrinkage prior is automatically extended to the corresponding parameter classes "sds" (spline SDs) and/or "sdgp" (GP SDs) using the brms-canonical main = TRUE pattern. The resulting prior regularises ALL components jointly – linear coefficients, nonlinear smooth wiggliness, and GP marginal variance – which is the principled approach to global-local shrinkage in models that combine parametric and nonparametric components.
hs_df	Numeric $> 0$. Local half-$t$ degrees of freedom for the horseshoe prior (default 1 = half-Cauchy).
hs_df_global	Numeric $> 0$. Global half-$t$ degrees of freedom (default 1).
hs_df_slab	Numeric $> 0$. Slab half-$t$ degrees of freedom (default 4).
hs_scale_global	Numeric $> 0$ or NULL. Scale for the global half-$t$ prior. NULL (default) lets brms compute it automatically from the number of predictors via hs_autoscale.
hs_scale_slab	Numeric $> 0$. Scale for the slab component (default 2).
hs_par_ratio	Numeric $> 0$ or NULL. Expected ratio of non-zero to total coefficients. NULL (default) treats all coefficients as potentially non-zero.
hs_autoscale	Logical. Whether brms should auto-scale the horseshoe prior using the residual SD $\sigma$. Default TRUE; set to FALSE for non-continuous responses (binomial, Poisson, ...) where $\sigma$ is not defined.
r2d2_mean_R2	Numeric in $(0, 1)$. Prior mean of the model $R^2$ (default 0.5).
r2d2_prec_R2	Numeric $> 0$. Prior precision of $R^2$ (default 2). Higher values concentrate mass around r2d2_mean_R2.
r2d2_cons_D2	Numeric $> 0$ or NULL. Dirichlet concentration for the D2 component. NULL (default) uses the brms default 0.5.
r2d2_autoscale	Logical. Whether brms should auto-scale the R2D2 prior using $\sigma$. Default TRUE; set to FALSE for non-continuous responses.
nonlinear	Character vector or NULL. Names of predictor variables to model with a smooth nonlinear term. Each listed variable is replaced in the formula RHS with s(var) (spline) or gp(var) (Gaussian process). Variables not listed remain linear. Do not also write s(x) in the formula when using this argument – the modification is applied automatically. Default NULL (all predictors remain linear).
nonlinear_type	Character. Smooth term family to use. One of "spline" (default, penalised regression spline via mgcv::s()) or "gp" (Gaussian process via brms::gp()).
spline_k	Integer. Spline basis dimension (number of knots) passed to mgcv::s(..., k = ...). -1L (default) lets mgcv choose automatically. For SAE typically k = 8 to 15. Ignored when nonlinear_type = "gp".
spline_bs	Character. Spline basis type passed to mgcv::s(..., bs = ...). Default "tp" (thin-plate regression spline, the mgcv default). Common alternatives: "cr" (cubic regression spline; often more stable for SAE with correlated auxiliary variables), "cs" (cubic with shrinkage; allows variable selection), "ps" (P-splines). Ignored when nonlinear_type = "gp".
gp_k	Integer or NA. Number of basis functions for the Hilbert-space approximate GP (Riutort-Mayol et al.\ 2020), passed to brms::gp(..., k = ...). NA (default) = exact GP which scales $O(n^3)$ and is not recommended for $n > 100$ areas; an immediate warning is emitted in that case pointing to this argument. Integer values 10–25 are typical for SAE and dramatically improve convergence and runtime. Ignored when nonlinear_type = "spline".
gp_cov	Character. GP covariance function passed to brms::gp(..., cov = ...): "exp_quad" (squared exponential / RBF, default), "matern15" (Matern 3/2), "matern25" (Matern 5/2; often more numerically stable for SAE than RBF), or "exponential".
gp_c	Numeric $> 0$ or NULL. Hilbert-space GP boundary-scale factor passed to brms::gp(..., c = ...). Default brms value is $5/4$ (= 1.25); increase if the GP appears truncated at the domain boundaries. Only relevant when gp_k is supplied.
gp_scale	Deprecated. Use gp_c instead. The old name suggested a length-scale interpretation but actually mapped to the HSGP boundary-scale factor. Will be removed in v2.0.0.
handle_missing	Character or NULL. Strategy for missing data. One of "deleted", "multiple", or "model" (see Details). If NULL (default) and missing values exist in the data, an informative error is raised.
m	Integer. Number of imputations when handle_missing = "multiple" (default: 5). Ignored for other strategies.
mice_args	A named list of additional arguments forwarded to mice, for example list(method = "pmm", seed = 42). Only used when handle_missing = "multiple".
measurement_error	Optional named list specifying which auxiliary variables are measured with error and where to find their standard errors. The list maps variable names to columns in data containing the SE, e.g. measurement_error = list(x1 = "se_x1", x2 = "se_x2"). Listed variables are wrapped on-the-fly with mi(var, se_col) in the brmsformula so that brms treats them as latent variables with a Gaussian measurement-error structure, following Ybarra and Lohr (2008). Standard errors must be non-negative and have no missing values; measurement_error variables must be a subset of the model's auxiliary (linear) predictors. When the user has already written mi(...) explicitly in the formula, the corresponding entries in measurement_error are detected and not duplicated. Note that ME inflates the parameter space (each x_i becomes latent), which slows sampling and increases the risk of divergent transitions, especially when combined with smooth terms (nonlinear). See the "Measurement error" section of the SAE vignette.
control	A named list of sampler control parameters (default: list()). Passed directly to brm.
chains	Integer. Number of MCMC chains (default: 4).
iter	Integer. Total iterations per chain (default: 4000).
warmup	Integer. Warm-up iterations per chain (default: floor(iter / 2)).
cores	Integer. Number of CPU cores for parallel sampling (default: 1).
sample_prior	Character. Whether to draw from the prior distribution. One of "no" (default), "yes", or "only".
sre	Deprecated. Use spatial_var instead. Kept for backward compatibility; will be removed in v2.0.0.
sre_type	Deprecated. Use spatial_model instead. Kept for backward compatibility; will be removed in v2.0.0.
stanvars	An optional stanvar object (or a list of such objects) supplying additional Stan code, data, or parameters. Passed directly to brm and brm_multiple. Intended for use by wrapper functions such as hbm_betalogitnorm that require custom Stan blocks; end users typically do not need to set this argument. Default: NULL.
...	Additional arguments forwarded to brm or brm_multiple.

Details

Hierarchical Bayesian Model for Small Area Estimation

**Spatial Small Area Estimation Models**

For spatially correlated areas, hbsaems extends the standard area-level SAE model (Fay-Herriot 1979) by adding a spatially structured random effect:

$y_i = x_i^\top \boldsymbol{\beta} + u_i + e_i,$

where $u_i$ is the spatial random effect for area $i$ and $e_i$ the sampling error. Two families of spatial structures are supported.

CAR (Conditional Autoregressive; Besag 1974)

Specified by spatial_model = "car". The joint distribution of the spatial effects is

$u \sim \mathcal{N}\bigl(0,\, \sigma_u^2 (D - \rho W)^{-1}\bigr),$

where $W$ is a binary adjacency matrix (1 if neighbour, 0 otherwise) and $D = \mathrm{diag}(W \mathbf{1})$. Sub-types via car_type: "icar" (intrinsic, $\rho = 1$; Besag 1991); "escar", "esicar" (exact sparse formulations of Morris et al.\ 2019); "bym2" (BYM2 reparameterisation of Riebler et al.\ 2016, recommended for disconnected graphs).

SAR (Simultaneous Autoregressive; Whittle 1954, Anselin 1988)

Specified by spatial_model = "sar". The model is

$u = \rho W u + \varepsilon, \quad \varepsilon \sim \mathcal{N}(0, \sigma_\varepsilon^2 I),$

where $W$ is row-standardised so that $\rho \in (-1, 1)$ carries an interpretable correlation meaning. Sub-types via sar_type: "lag" (spatial lag of the response, $y = \rho W y + X\boldsymbol{\beta} + \varepsilon$); "error" (spatial error model).

Use check_spatial_weight to verify that $M$ satisfies the theoretical requirements (square, zero diagonal, symmetry for CAR, style-appropriate for the model class). Use build_spatial_weight to construct $M$ from a shapefile.

**Missing Data Strategies**

The three strategies differ in scope and statistical assumptions:

"deleted"

Complete-case analysis. Only rows where all response variable(s) are observed are used for model fitting. Auxiliary variables must be complete; otherwise an informative error is raised. Appropriate under MCAR (Missing Completely At Random).

"multiple"

Multiple imputation via mice for auxiliary (predictor) variables only. The response variable $Y$ is never imputed. In a Bayesian model, missing outcomes are naturally marginalised through the posterior predictive distribution:

$p(\theta \mid Y_{\text{obs}}, X) = \int p(\theta \mid Y_{\text{obs}}, Y_{\text{mis}}, X)\, p(Y_{\text{mis}} \mid Y_{\text{obs}}, X, \theta)\, \mathrm{d}Y_{\text{mis}}.$

Imputing $Y$ before fitting would replace this integral with a single point substitute, deflate posterior uncertainty, and potentially bias the estimates if the imputation model is misspecified. If $Y$ has missing values, those rows are excluded from model fitting (a warning is issued) but are retained in the returned object for subsequent prediction via hbsae. Appropriate under MAR (Missing At Random).

"model"

Model-based imputation using brms::mi(). Missing values in auxiliary variables are jointly estimated with the model parameters. The user must specify imputation sub-models explicitly in the formula argument, e.g.: bf(y | mi() ~ mi(x1) + x2) + bf(x1 | mi() ~ x2). Only applicable to continuous distributions. Appropriate under MAR.

If data are Missing Not At Random (MNAR), none of the above strategies applies directly; sensitivity analyses and explicit missingness models are recommended.

Value

An object of class hbmfit, which is a named list containing:

model	The fitted brmsfit object (or brmsfit_multiple when handle_missing = "multiple" with missing predictors).
handle_missing	The missing-data strategy used (NULL if the data were complete).
data	The original data frame passed to hbm() before any row deletion or imputation. This is intentional: hbsae needs all rows – including those with missing $Y$ – to generate predictions for every small area.

How to pin distributional parameters per family

hbsaems exposes three layers for pinning distributional parameters to known values, in increasing order of generality:

1. Family-specific sugar (least typing, most readable). Each wrapper exposes a convenience argument that maps to a well-defined Fay-Herriot-style transformation of survey design quantities:

   Wrapper	Sugar argument	Pinned dpar
   hbm(..., hb_sampling = "gaussian")	sampling_variance = "D"	$\sigma_i = \sqrt{D_i}$
   hbm_lnln()	sampling_variance = "psi"	$\sigma_i = \sqrt{\psi_i}$ (on log scale)
   hbm_betalogitnorm()	n = "n", deff = "deff"	$\phi_i = n_i / \mathrm{deff}_i - 1$ (Liu 2009)
   hbm_binlogitnorm()	trials = "n"	(sampling variance built into the family)

2. Universal fixed_params (works everywhere). A named list pinning any distributional parameter – accepts a column name (character), a scalar (numeric of length 1), a vector of length nrow(data), or a one-sided formula (evaluated in data's environment). Examples:

   • fixed_params = list(sigma = "D") – pin sigma to a column.

   • fixed_params = list(phi = ~ I(n / deff - 1)) – pin phi via formula.

   • fixed_params = list(nu = 4) – pin Student-t df to a scalar.

3. Raw stanvars (for power users authoring custom Stan code). Direct injection of Stan code blocks – see stanvar.

Sugar arguments are simply thin wrappers around layer 2: they validate the survey-design inputs (no NAs, strictly positive) and translate to fixed_params before delegating to the universal machinery. Hence the conflict policy below applies uniformly to both sugar and explicit fixed_params.

Conflict resolution between prior, prior_type, fixed_params, and stanvars

hbm() provides four orthogonal mechanisms to influence the prior / parameter specification of the underlying brms model:

1. prior – explicit brmsprior object(s).

2. prior_type – global-local shrinkage prior on the regression coefficients (cascades to "sds" / "sdgp" when splines or GPs are present).

3. fixed_params – pin distributional parameters to known values via the offset trick.

4. stanvars – inject custom Stan code blocks.

Plus two family-specific sugar arguments that translate to fixed_params internally:

• sampling_variance (continuous families): pins $\sigma_i = \sqrt{D_i}$, equivalent to fixed_params$sigma = sqrt(data$D).

• n + deff (hbm_betalogitnorm): pins $\phi_i = n_i / \mathrm{deff}_i - 1$, equivalent to fixed_params$phi = n / deff - 1.

Combining these without rules in mind can produce unidentified models or compile-time errors. hbm() therefore enforces the following conflict matrix, where each cell describes what happens when the row and column inputs both target the same distributional parameter (e.g.\ both pin sigma):

	fixed_params	prior	prior_type	stanvars
sampling_variance	error	error (transitive)	no overlap	error (transitive)
n + deff	error	error (transitive)	no overlap	error (transitive)
fixed_params	--	error (10b.i)	no overlap	error (10b.ii)
prior	error (10b.i)	--	warning, user wins	no check needed
prior_type	no overlap	warning, user wins	--	no check needed
stanvars	error (10b.ii)	no check needed	no check needed	--

Resolution semantics in detail:

• prior vs prior_type. If the user supplies a global (no coef =) prior on class = "b", "sds", or "sdgp", prior_type is silently dropped for that class and a warning is emitted. Coefficient-specific user priors (coef = "x1") are kept alongside the shrinkage prior without warning.

• fixed_params vs prior. A pinned parameter is removed from the sampler; supplying a prior on that same parameter therefore has no effect and is treated as a user error – an informative stop() is issued.

• fixed_params vs stanvars. Same logic as above: a sampling statement in stanvars that targets a pinned parameter would fail at Stan compile time; hbm() catches this and stops with a clear message.

• Sugar vs fixed_params on the same parameter. The sugar translators (.translate_sampling_variance(), .translate_n_deff_to_phi()) error out if the user has also pre-populated fixed_params for the target parameter – there should never be two pin sources for the same dpar.

• Sugar vs prior / stanvars transitively. After the sugar -> fixed_params translation, the downstream fixed_params-vs-prior / -stanvars checks fire automatically. E.g.\ sampling_variance = "D" plus prior = set_prior("normal(0, 1)", class = "sigma") errors via the fixed_params vs prior rule.

The intent is to fail fast and explicitly rather than silently producing an unidentified or mis-specified model.

Convergence advice for SAE practitioners

Common convergence pathologies in hierarchical Bayesian SAE models and how to address them. Run convergence_check() after fitting to inspect $\hat R$, effective sample size (ESS), and divergent transitions.

1. Default sampler settings (recommended starting point).

• chains = 4, iter = 4000, warmup = 2000

• control = list(adapt_delta = 0.95, max_treedepth = 12)

• cores = parallel::detectCores() - 1

2. Divergent transitions. Most common cause is the funnel geometry of hierarchical variance parameters.

• First-line: increase adapt_delta to 0.99.

• If still diverging, increase warmup and consider a tighter prior on the area-level standard deviation (sd class), e.g.\ set_prior("normal(0, 0.5)", class = "sd").

• For Beta/Binomial logit-normal models, prior on the random intercept SD should also be on the logit scale.

• Gaussian (Fay-Herriot) only. Always supply the known sampling variance via sampling_variance = "<col>". Without this, the residual $\sigma$ and the area-RE $\sigma_u$ compete to explain the same variance component, producing weak identifiability and divergent transitions almost regardless of adapt_delta. This is the single most common cause of divergences in Fay-Herriot fits and should be checked first before any sampler-tuning.

3. Low effective sample size (ESS < 1000).

• Increase iter (e.g.\ to 6000); this is the single most reliable fix.

• Centre and scale the auxiliary variables before fitting.

• Check for prior–data conflict via priorsense; see ?prior_check.

4. Gaussian processes.

• Exact GP scales $O(n^3)$. For $n > 100$ areas, set gp_k to use the Hilbert-space approximate GP (Riutort-Mayol et al.\ 2020). A heuristic is gp_k = ceiling(min(n / 5, 25)).

• Try gp_cov = "matern25" (Matern 5/2) if the default squared-exponential covariance is numerically unstable.

• The boundary-scale factor gp_c (brms default 1.25) may need increasing if the posterior GP is truncated at the domain edges.

5. Splines.

• Start with spline_k = -1 (auto). Increase only if the residual diagnostics suggest under-smoothing.

• For strongly correlated auxiliary variables, try spline_bs = "cr" (cubic regression spline) for better numerical stability than the default thin-plate.

• For variable selection, use spline_bs = "cs" (cubic with shrinkage); coefficients on irrelevant smooths shrink toward zero.

6. Spatial models.

• For CAR, car_type = "bym2" (Riebler et al.\ 2016) is the modern recommendation; it stabilises the spatial/IID decomposition via a single mixing parameter.

• Verify the weight matrix with check_spatial_weight(); isolated areas or multiple disconnected components cause non-identifiability.

7. Prior predictive check first. Always call prior_check() (sample_prior = "only") before the full posterior run. Implausible prior predictives are the single most common cause of slow / divergent sampling.

Author(s)

Achmad Syahrul Choir, Saniyyah Sri Nurhayati, and Sofi Zamzanah

References

Rao, J. N. K., & Molina, I. (2015). Small Area Estimation. John Wiley & Sons.

Burkner, P. C. (2017). brms: An R package for Bayesian multilevel models using Stan. Journal of Statistical Software, 80(1), 1–28.

Riutort-Mayol, G., Burkner, P.-C., Andersen, M. R., Solin, A., & Vehtari, A. (2023). Practical Hilbert space approximate Bayesian Gaussian processes for probabilistic programming. Statistics and Computing, 33, 17. doi:10.1007/s11222-022-10167-2

Riebler, A., Sorbye, 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.

van Buuren, S., & Groothuis-Oudshoorn, K. (2011). mice: Multivariate imputation by chained equations in R. Journal of Statistical Software, 45(3), 1–67.

Examples

library(hbsaems)
library(brms)
data("data_fhnorm")
data <- data_fhnorm

# Standard brms-default MCMC settings used throughout these
# examples (chains = 4, iter = 2000, warmup = 1000).  For tougher
# posteriors (funnel geometry, weakly identified priors), bump to
# iter = 4000, warmup = 2000, control = list(adapt_delta = 0.99).
FAST <- list(chains = 4, iter = 2000, warmup = 1000, cores = 1,
             seed = 123, refresh = 0)

# -- Basic model --------------------------------------------------------------
model <- do.call(hbm, c(
  list(formula     = bf(y ~ x1 + x2 + x3),
       hb_sampling = "gaussian",
       hb_link     = "identity",
       re          = ~(1 | regency),
       data        = data),
  FAST
))
summary(model)

# -- Horseshoe prior (sparse coefficients) ------------------------------------
model_hs <- do.call(hbm, c(
  list(formula     = bf(y ~ x1 + x2 + x3),
       re          = ~(1 | regency),
       data        = data,
       prior_type  = "horseshoe",
       hs_df       = 1),
  FAST
))
summary(model_hs)

# -- R2D2 prior (prior on model R-squared) -------------------------------------
model_r2 <- do.call(hbm, c(
  list(formula      = bf(y ~ x1 + x2 + x3),
       re           = ~(1 | regency),
       data         = data,
       prior_type   = "r2d2",
       r2d2_mean_R2 = 0.5,
       r2d2_prec_R2 = 2),
  FAST
))
summary(model_r2)

# -- Spline smooth for x1 (nonlinear) -----------------------------------------
# x1 is modelled with s(x1); x2 and x3 remain linear.
model_spline <- do.call(hbm, c(
  list(formula        = bf(y ~ x1 + x2 + x3),
       re             = ~(1 | regency),
       data           = data,
       nonlinear      = "x1",
       nonlinear_type = "spline"),
  FAST
))
summary(model_spline)

# -- Gaussian process for x2 (nonlinear) --------------------------------------
model_gp <- do.call(hbm, c(
  list(formula        = bf(y ~ x1 + x2 + x3),
       re             = ~(1 | regency),
       data           = data,
       nonlinear      = "x2",
       nonlinear_type = "gp"),
  FAST
))
summary(model_gp)

# -- Missing data: deletion (Y missing, X complete) ---------------------------
data_miss_y        <- data
data_miss_y$y[3:5] <- NA

model_deleted <- do.call(hbm, c(
  list(formula        = bf(y ~ x1 + x2 + x3),
       re             = ~(1 | regency),
       data           = data_miss_y,
       handle_missing = "deleted"),
  FAST
))
summary(model_deleted)

# -- Missing data: multiple imputation (X only -- Y is never imputed) ----------
data_miss_x          <- data
data_miss_x$x1[6:8]  <- NA

model_multiple <- do.call(hbm, c(
  list(formula        = bf(y ~ x1 + x2 + x3),
       re             = ~(1 | regency),
       data           = data_miss_x,
       handle_missing = "multiple",
       m              = 5),
  FAST
))
summary(model_multiple)

# -- Missing data: model-based imputation with mi() ---------------------------
data_miss_x2         <- data
data_miss_x2$x1[6:7] <- NA

model_model <- do.call(hbm, c(
  list(formula        = bf(y | mi() ~ mi(x1) + x2 + x3) +
                        bf(x1 | mi() ~ x2 + x3),
       re             = ~(1 | regency),
       data           = data_miss_x2,
       handle_missing = "model"),
  FAST
))
summary(model_model)

# -- Spatial: CAR (Conditional Autoregressive) --------------------------------
data("adjacency_matrix_car")
model_car <- do.call(hbm, c(
  list(formula     = bf(y ~ x1 + x2 + x3),
       data        = data,
       spatial_var = "province",
       spatial_model    = "car",
       M           = adjacency_matrix_car),
  FAST
))
summary(model_car)

# -- Spatial: SAR (Simultaneous Autoregressive) -------------------------------
# spatial_weight_sar is a 100x100 row-standardised matrix with row-
# names regency_001..regency_100, so it pairs with the fine-grained
# "regency" column (100 levels) -- NOT with "province" (5 levels).
data("spatial_weight_sar")
model_sar <- do.call(hbm, c(
  list(formula     = bf(y ~ x1 + x2 + x3),
       data        = data,
       spatial_var = "regency",
       spatial_model    = "sar",
       M           = spatial_weight_sar),
  FAST
))
summary(model_sar)

---

Small Area Estimation Under a Beta Likelihood (Logit-Normal Link)

Description

Convenience wrapper around hbm for SAE problems where the response $y_i \in (0, 1)$ is modelled as

$y_i \mid \theta_i, \phi_i \sim \mathrm{Beta}(\theta_i \phi_i, (1 - \theta_i)\phi_i),$

$\mathrm{logit}(\theta_i) = x_i^\top \boldsymbol{\beta} + u_i, \quad u_i \sim \mathcal{N}(0, \sigma_u^2).$

Usage

hbm_betalogitnorm(
  response,
  auxiliary = NULL,
  data,
  n = NULL,
  deff = NULL,
  area_var = NULL,
  area_re_structure = c("nested", "crossed"),
  spatial_var = NULL,
  spatial_model = NULL,
  car_type = NULL,
  sar_type = NULL,
  M = NULL,
  link_phi = NULL,
  prior = NULL,
  stanvars = NULL,
  handle_missing = NULL,
  m = 5L,
  control = list(),
  chains = 4L,
  iter = 4000L,
  warmup = floor(iter/2),
  cores = 1L,
  sample_prior = "no",
  fixed_params = NULL,
  predictors = NULL,
  group = NULL,
  sre = NULL,
  sre_type = NULL,
  ...
)

Arguments

response	Character. Name of the response column (must lie strictly between 0 and 1).
auxiliary	Character vector of auxiliary (fixed-effect) variable names; corresponds to area-level covariates in SAE literature (see Rao & Molina 2015 Ch. 4).
data	A data frame.
n	Character or NULL. Name of the column giving the per-area sample size used to compute the fixed $\phi$. Must be supplied together with deff. When supplied, $\phi_i = n_i / \text{deff}_i - 1$ is pinned via offset. Default: NULL (treats phi as random with brms's default prior).
deff	Character or NULL. Name of the design-effect column. Required when n is supplied (and vice versa).
area_var	Character vector or NULL. Name(s) of the area-grouping column(s). Length 1 adds an IID random intercept (1 | area_var); length $\geq$ 2 supports hierarchical areas (e.g.\ c("province", "regency")) – see ?hbm_flex for the nested vs.\ crossed structures.
area_re_structure	Either "nested" (default) or "crossed"; controls how multiple area columns combine.
spatial_var, spatial_model, car_type, sar_type, M	Spatial random-effect arguments, forwarded to hbm.
link_phi	Character or NULL. Link function for phi. Default NULL resolves automatically: "identity" when phi is pinned via the survey design (the n + deff sugar or fixed_params$phi, both of which produce raw positive offsets) and "log" otherwise (the brms default for phi in the Beta family; keeps $\phi > 0$ while letting NUTS sample on $\mathbb{R}$). Manually setting "identity" together with an estimated phi can let NUTS propose negative values and trigger divergent transitions; an informative warning is emitted in that case.
prior	Optional brmsprior object. If NULL, sensible defaults are filled in: Intercept ~ student_t(4, 0, 10) b ~ student_t(4, 0, 2.5) phi – brms default gamma(0.01, 0.01) (in random mode; ignored in fixed mode). The user may pass a partial prior: missing default classes are filled in automatically. To put a custom prior on phi (random mode), set prior = set_prior("gamma(2, 0.5)", class = "phi") or similar.
stanvars	Optional brmsstanvars object for power users who need to inject custom Stan code blocks (e.g.\ transformed data, generated quantities, model-block statements not expressible via the prior argument). Passed through to brms verbatim. Note: As of v1.0.0, this wrapper no longer declares alpha or beta as Stan parameters, so legacy stanvars blocks containing sampling statements on alpha / beta (left over from the pre-v1.0.0 hierarchical-phi construction) will now raise an informative error.
handle_missing, m, control, chains, iter, warmup, cores, sample_prior, ...	Passed through to hbm.
fixed_params	Optional named list pinning distributional parameters to known values. See hbm for the spec format. Allows power-user access to the same machinery used by the n/deff arguments.
predictors	Deprecated. Use auxiliary instead. Kept for backward compatibility; will be removed in v2.0.0.
group	Deprecated. Use area_var instead.
sre	Deprecated. Use spatial_var instead.
sre_type	Deprecated. Use spatial_model instead.

Details

The precision parameter $\phi_i$ can either be pinned to a survey design effect (n + deff) or sampled with brms's default weakly-informative prior $\phi \sim \mathrm{Gamma}(0.01, 0.01)$.

Value

An object of class hbmfit.

Migration note (v1.0.0)

Earlier versions of this wrapper introduced a hierarchical construction $\phi \sim \mathrm{Gamma}(\alpha, \beta),\, \alpha \sim \mathrm{Gamma}(1,1),\, \beta \sim \mathrm{Gamma}(1,1)$ for the random-phi mode, declaring alpha and beta as Stan parameters with hyperpriors injected via stanvars. Starting v1.0.0, that construction has been removed in favour of brms's own default prior, $\phi \sim \mathrm{Gamma}(0.01, 0.01)$ with lower bound 0 (mean 1, variance 100; weakly informative on the precision scale). Three reasons:

• Brittle priors on alpha/beta. The hyperprior $\mathrm{Gamma}(1,1)$ on $\alpha$ has prior mean 1, which is on the boundary of the parameter space declared as real<lower=1>. This routinely produced divergent transitions when the data were not informative about phi.

• Parameter blow-up. Estimating two extra Stan parameters per area-level model with limited data inflated the effective posterior dimension and slowed convergence.

• Cleaner brms semantics. Letting brms apply its own default $\mathrm{Gamma}(0.01, 0.01)$ means that passing prior = NULL now does exactly what the user expects: “brms defaults”. No surprises.

If you need to reproduce the old behaviour, supply stanvars yourself to declare alpha, beta and their hyperpriors, and pass prior = set_prior("gamma(alpha, beta)", class = "phi"). Pre-v1.0.0 code that supplied stanvars with hyperpriors on alpha/beta will now raise an informative error.

Conflict policy

When the precision parameter $\phi$ is pinned via n + deff (or via fixed_params$phi), the function refuses any additional specification that would also set $\phi$. Specifically, all of the following are rejected with an informative error at construction time:

• n supplied without deff, or vice versa.

• n + deff and fixed_params$phi.

• n + deff and a user prior on class = "phi".

• n + deff and a stanvars hyperprior on alpha or beta (the $\mathrm{Gamma}(\alpha, \beta)$ hyperparameters used in random-$\phi$ mode).

• auxiliary and the deprecated predictors in the same call.

References

Liu, B. (2009). Hierarchical Bayes Estimation and Empirical Best Prediction of Small-Area Proportions. University of Maryland.

Rao, J. N. K., & Molina, I. (2015). Small Area Estimation, 2nd ed. Wiley, p. 390.

See Also

hbm_flex, hbm

Examples

library(hbsaems)
library(brms)
data("data_betalogitnorm")
data <- data_betalogitnorm

# -- 1. Basic model (phi random, with default hyperprior) --------------------
model1 <- hbm_betalogitnorm(
  response   = "y",
  auxiliary  = c("x1", "x2", "x3"),
  area_var   = "regency",
  data       = data,
  chains = 4, iter = 2000, warmup = 1000, refresh = 0
)
summary(model1)

# -- 2. Fixed phi via survey design (n + deff) -------------------------------
model2 <- hbm_betalogitnorm(
  response   = "y",
  auxiliary  = c("x1", "x2", "x3"),
  n          = "n",
  deff       = "deff",
  area_var   = "regency",
  data       = data,
  chains = 4, iter = 2000, warmup = 1000, refresh = 0
)
summary(model2)

# -- 3. Custom prior on phi via the standard `prior` argument ----------------
#
# Starting v1.0.0, hbm_betalogitnorm() no longer declares its own
# alpha / beta hyperparameters; phi uses brms's native default
# Gamma(0.01, 0.01).  To override that prior, pass a `brms::set_prior()`
# entry to the `prior` argument -- the legacy
# stanvar("alpha ~ gamma(...); beta ~ gamma(...)") pattern would now
# error because alpha/beta are no longer Stan parameters.
model3 <- hbm_betalogitnorm(
  response   = "y",
  auxiliary  = c("x1", "x2", "x3"),
  area_var   = "regency",
  data       = data,
  prior      = brms::set_prior("gamma(2, 0.5)", class = "phi"),
  chains = 4, iter = 2000, warmup = 1000, refresh = 0
)

# -- 4. Spatial CAR model ----------------------------------------------------
data("adjacency_matrix_car")
model4 <- hbm_betalogitnorm(
  response   = "y",
  auxiliary  = c("x1", "x2", "x3"),
  n          = "n",
  deff       = "deff",
  spatial_var = "province",
  spatial_model   = "car",
  M          = adjacency_matrix_car,
  data       = data,
  chains = 4, iter = 2000, warmup = 1000, refresh = 0
)

---

Small Area Estimation under a Binomial Logit-Normal Model

Description

Convenience wrapper that fits a Hierarchical Bayesian Small Area Estimation model with a binomial likelihood and logit link. Internally delegates to hbm_flex with family_key = "binomial".

Usage

hbm_binlogitnorm(
  response,
  trials,
  auxiliary = NULL,
  data,
  area_var = NULL,
  area_re_structure = c("nested", "crossed"),
  spatial_var = NULL,
  spatial_model = NULL,
  car_type = NULL,
  sar_type = NULL,
  M = NULL,
  fixed_params = NULL,
  predictors = NULL,
  group = NULL,
  sre = NULL,
  sre_type = NULL,
  ...
)

Arguments

response	Character. Name of the successes variable (non-negative integer, $y \le n$).
trials	Character. Name of the trials variable (positive integer).
auxiliary	Character vector of auxiliary (fixed-effect) variable names; corresponds to area-level covariates in SAE literature (see Rao & Molina 2015 Ch. 4).
data	A data.frame.
area_var	Optional character vector. Name(s) of column(s) in data identifying the small area / domain. Length 1 fits an IID area-level random intercept (1 | area_var); length $\geq$ 2 supports hierarchical areas – see ?hbm_flex for the nested vs.\ crossed structures. Default: NULL.
area_re_structure	Either "nested" (default) or "crossed"; controls how multiple area columns combine.
spatial_var	Optional character. Name of a column identifying the spatial cluster. Must be supplied together with spatial_model and M. Default: NULL.
spatial_model	Optional character. Spatial dependence: "car" (default car_type = "icar") or "sar" (default sar_type = "lag"). Default: NULL.
car_type	Optional character. CAR sub-type passed to brms: "escar", "esicar", "icar", or "bym2".
sar_type	Optional character. SAR sub-type: "lag" or "error".
M	Optional numeric matrix. Spatial weight matrix. Required when spatial_model is supplied.
fixed_params	Optional named list pinning distributional parameters to known values. See hbm for the spec format.
predictors	Deprecated. Use auxiliary instead. Kept for backward compatibility; will be removed in v2.0.0.
group	Deprecated. Use area_var instead.
sre	Deprecated. Use spatial_var instead.
sre_type	Deprecated. Use spatial_model instead.
...	Additional arguments forwarded to hbm_flex (e.g.\ prior_type, handle_missing, sampler controls).

Details

Let $y_i$ denote the number of successes in area $i$ out of $n_i$ trials. The model is

$y_i \mid p_i, n_i \sim \mathrm{Binomial}(n_i, p_i),$

$\mathrm{logit}(p_i) = x_i^\top \boldsymbol{\beta} + u_i, \quad u_i \sim \mathcal{N}(0, \sigma_v^2).$

Value

An object of class hbmfit.

Conflict policy

• auxiliary and the deprecated predictors in the same call are rejected with an informative error.

• handle_missing = "model" is not supported (binomial is a discrete family; see Notes on missing data).

Notes on missing data

The binomial family does not support handle_missing = "model" (joint Bayesian imputation via brms::mi()). When NA values are detected and handle_missing is left NULL, the wrapper auto-selects "multiple" (multiple imputation via mice on the predictors only).

See Also

hbm_flex, hbm

Examples

library(hbsaems)
library(brms)
data("data_binlogitnorm")

# -- 1. Standard binomial logit-normal SAE ---------------------------------
fit <- hbm_binlogitnorm(
  response   = "y",
  trials     = "n",
  auxiliary  = c("x1", "x2", "x3"),
  area_var   = "district",        # area-level random effect (1 | district)
  data       = data_binlogitnorm,
  chains = 4, iter = 2000, warmup = 1000, refresh = 0
)

# -- 2. With spatial CAR random effect -------------------------------------
data("adjacency_matrix_car_regency")
fit_car <- hbm_binlogitnorm(
  response      = "y",
  trials        = "n",
  auxiliary     = c("x1", "x2", "x3"),
  spatial_var   = "regency",
  spatial_model = "car",
  M             = adjacency_matrix_car_regency,
  data          = data_binlogitnorm,
  chains = 4, iter = 2000, warmup = 1000, refresh = 0
)

---

Sampler Configuration for HBSAE Models

Description

Bundles the MCMC sampler arguments of hbm (and the hbm_* family wrappers) into a single named list. Use this when you want a reusable sampler profile or to reduce the size of long hbm() calls.

Usage

hbm_control(
  chains = 4L,
  iter = 4000L,
  warmup = NULL,
  thin = 1L,
  cores = 1L,
  seed = NULL,
  refresh = NULL,
  adapt_delta = NULL,
  max_treedepth = NULL,
  control = NULL
)

Arguments

chains	Integer. Number of Markov chains (default 4L).
iter	Integer. Total iterations per chain (default 4000L).
warmup	Integer. Warm-up iterations per chain. Default floor(iter / 2).
thin	Integer. Thinning interval (default 1L).
cores	Integer. Number of cores for parallel chains (default 1L).
seed	Optional integer seed for reproducibility.
refresh	Integer. Stan progress refresh frequency (default NULL: brms default).
adapt_delta	Numeric in $(0, 1)$. When supplied, included in the control list as control = list(adapt_delta = ...).
max_treedepth	Integer. Max tree depth, included in control when supplied.
control	Optional list of additional NUTS control options. Merged with any adapt_delta / max_treedepth above.

Details

This is entirely opt-in: the flat signatures of hbm(), hbm_lnln(), etc.\ continue to work exactly as before. Pass the result directly to hbm() – it is auto-spliced via ....

Value

A named list whose elements are valid arguments of hbm.

See Also

hbm_priors, hbm_nonlinear, hbm

Examples

# Build a reusable "high-quality" profile
hq <- hbm_control(chains = 4, iter = 8000, cores = 4,
                  adapt_delta = 0.99, seed = 1)
str(hq)

# Quick draft profile
draft <- hbm_control(chains = 2, iter = 1000)

---

Return the Data Used to Fit an hbmfit

Description

Return the Data Used to Fit an hbmfit

Usage

hbm_data(model)

Arguments

model

An hbmfit object.

Value

The original data.frame passed to the fitting function.

Examples

library(hbsaems)
library(brms)
data("data_fhnorm")
model <- hbm(brms::bf(y ~ x1), data = data_fhnorm,
             re = ~ (1 | regency),    # area-level random effect
             chains = 4, iter = 2000, warmup = 1000,
             cores = 1, seed = 1, refresh = 0)
head(hbm_data(model))

---

Fit a Flexible HBSAE Model with Any Registered Family

Description

Flexible factory that fits an HBSAE model using any distribution currently registered in the family registry, together with the full set of cross-cutting features: spatial random effects (CAR/SAR), shrinkage priors (horseshoe, R2D2), smooth terms (penalised splines, Gaussian processes), auxiliary-parameter hyperpriors, and missing-data strategies. Distribution-specific wrappers (hbm_lnln, hbm_binlogitnorm, hbm_betalogitnorm) are thin signature shims around this function with preset family_key values. Advanced users can also call it directly once a custom family has been registered via register_hbsae_model.

Usage

hbm_flex(
  family = NULL,
  response,
  auxiliary = NULL,
  data,
  family_key = NULL,
  addition_var = NULL,
  area_var = NULL,
  area_re_structure = c("nested", "crossed"),
  spatial_var = NULL,
  spatial_model = NULL,
  car_type = NULL,
  sar_type = NULL,
  M = NULL,
  prior = NULL,
  fixed_params = NULL,
  sampling_variance = NULL,
  prior_type = "default",
  hs_df = 1,
  hs_df_global = 1,
  hs_df_slab = 4,
  hs_scale_global = NULL,
  hs_scale_slab = 2,
  hs_par_ratio = NULL,
  hs_autoscale = TRUE,
  r2d2_mean_R2 = 0.5,
  r2d2_prec_R2 = 2,
  r2d2_cons_D2 = NULL,
  r2d2_autoscale = TRUE,
  nonlinear = NULL,
  nonlinear_type = "spline",
  spline_k = -1L,
  spline_bs = "tp",
  gp_k = NA_integer_,
  gp_cov = "exp_quad",
  gp_c = NULL,
  gp_scale = NULL,
  handle_missing = NULL,
  m = 5L,
  mice_args = list(),
  control = list(),
  chains = 4L,
  iter = 4000L,
  warmup = floor(iter/2),
  cores = 1L,
  sample_prior = "no",
  link = NULL,
  aux_args = NULL,
  stanvars = NULL,
  predictors = NULL,
  group = NULL,
  sre = NULL,
  sre_type = NULL,
  ...
)

Arguments

family	The model family (primary argument). Accepts EITHER a character family name (e.g.\ "gaussian", "binomial") OR a brms family object (e.g.\ gaussian(), Gamma()), OR a registered custom family object (e.g.\ brms_custom_loglogistic()$custom_family). The name is matched to a registered spec, so a custom family reuses the same Stan code (stanvars) as the key string would. A custom family object must already be registered via register_hbsae_brms_custom; an unregistered custom family raises an informative error. Supply either family or its alias family_key, not both. family is the uniform way to choose the distribution across hbm and hbm_flex(); each function additionally keeps its historical alias (family_key here, hb_sampling in hbm).
response	Character. Name of the response variable column.
auxiliary	Character vector. Names of auxiliary (fixed-effect) variables; corresponds to area-level covariates in SAE literature.
data	A data.frame.
family_key	Character. Backward-compatible alias for family (the registry key of the desired family). Not deprecated; pass EITHER family (preferred) OR family_key, never both (e.g.\ "lognormal", "binomial", "gamma_log").
addition_var	Character or NULL. Name of the addition term variable (e.g.\ trials for binomial). Required when the family spec has has_addition_term = TRUE.
area_var	Character vector or NULL. Name(s) of the column(s) in data identifying the areas. Three usage modes: Length 1 (default behaviour): a single area-level random intercept (1 | area_var). Length $\geq$ 2 with area_re_structure = "nested" (default): a hierarchy of areas given from the highest to the lowest level, e.g. c("province", "regency") yields (1 | province / regency) which brms expands to (1 | province) + (1 | province:regency). This is the canonical multi-stage SAE setup. Length $\geq$ 2 with area_re_structure = "crossed": non-nested levels, e.g.\ adding separate effects for (1 | province) + (1 | urbanrural). Use only when the levels are truly crossed rather than hierarchically nested.
area_re_structure	Either "nested" (default) or "crossed". Only consulted when area_var has length $\geq$ 2. See above.
spatial_var, spatial_model, car_type, sar_type, M	Spatial random-effect arguments forwarded to hbm. See ?hbm for a full description.
prior, prior_type, hs_df, hs_df_global, hs_df_slab, hs_scale_global, hs_scale_slab, hs_par_ratio, hs_autoscale, r2d2_mean_R2, r2d2_prec_R2, r2d2_cons_D2, r2d2_autoscale	Shrinkage prior arguments forwarded to hbm. When the formula contains s() or gp() terms, the prior is automatically cascaded to the corresponding parameter classes ("sds", "sdgp") via the brms main = TRUE pattern.
fixed_params	Optional named list pinning distributional parameters to known values. See hbm for the spec format. Allows power-user access to the generic fixed-parameter machinery (works with custom families too).
sampling_variance	Optional character. Name of a column in data containing the known sampling variance $D_i$ for each area (the Fay-Herriot sugar). When supplied, $\sigma_i = \sqrt{D_i}$ is pinned via offset. Forwarded to hbm, where a family-compatibility check ensures the family exposes a residual SD parameter named sigma (gaussian, lognormal, student, skew_normal, exgaussian, asym_laplace). Incompatible families (beta, binomial, poisson, etc.) raise an explicit error pointing at the family-specific alternative. See ?hbm for details.
nonlinear	Character vector of variable names to include as smooth/nonlinear terms (rather than linear). Variables listed here that also appear in auxiliary are modelled nonlinearly only.
nonlinear_type	Character. "spline" (penalised regression spline via mgcv, default) or "gp" (Gaussian process via brms).
spline_k	Integer. Spline basis dimension passed to mgcv::s(..., k = ...). -1 (default) lets mgcv choose automatically. For SAE typically k = 8 to 15.
spline_bs	Character. Spline basis type passed to mgcv::s(..., bs = ...). Defaults to "tp" (thin-plate regression spline). Set "cr" (cubic regression) for better numerical stability with correlated auxiliary variables. Other choices: "cs" (cubic with shrinkage), "ps" (P-splines).
gp_k	Integer or NA. Number of basis functions for the Hilbert-space approximate GP (Riutort-Mayol et al.\ 2020). NA (default) = exact GP, scales $O(n^3)$ and is not recommended for $n > 100$ areas. Integer gp_k = 10–25 is typical for SAE applications and dramatically improves convergence and runtime.
gp_cov	Character. Covariance function: "exp_quad" (squared exponential / RBF, default), "matern15" (Matern 3/2), "matern25" (Matern 5/2, often more numerically stable than RBF), "exponential".
gp_c	Numeric. Hilbert-space GP boundary-scale factor passed to brms::gp(c = ...). Default brms value is 5/4 (= 1.25); increase if the GP appears truncated at domain boundaries. Only relevant when gp_k is set.
gp_scale	Deprecated. Use gp_c instead. The old name suggested a length-scale interpretation but actually mapped to the boundary-scale factor.
handle_missing, m, mice_args	Missing-data arguments. When handle_missing = NULL the wrapper auto-selects a strategy based on the family registry's supports_mi flag.
control, chains, iter, warmup, cores, sample_prior, link	Sampler and model-spec arguments forwarded to hbm.
aux_args	Optional named list of family-specific auxiliary arguments (e.g.\ list(n = "n", deff = "deff") for the Beta family's phi hyperprior). Forwarded to the family's aux_param_hyperprior callback if it has one.
stanvars	Optional stanvar object passed through to brms. When the family has an aux_param_hyperprior callback that returns its own stanvars, the two are concatenated.
predictors	Deprecated. Use auxiliary instead. Kept for backward compatibility; will be removed in v2.0.0.
group	Deprecated. Use area_var instead.
sre	Deprecated. Use spatial_var instead.
sre_type	Deprecated. Use spatial_model instead.
...	Additional arguments forwarded to brm.

Details

The factory performs five duties that were previously duplicated across wrappers:

1. Validate that response, auxiliary, and optional variables exist in data.

2. Run the family's response_check on the response and report a human-readable error on failure.

3. Auto-select a missing-data strategy that respects the family's supports_mi flag (e.g.\ binomial cannot use "model").

4. Build the brms formula with optional addition terms and apply spline/GP transformations.

5. Invoke the family's aux_param_hyperprior callback (if defined) so distributions with a hyperprior on auxiliary parameters – e.g.\ phi for the Beta family, shape for Gamma, nu for Student-t – can inject Stan code without writing a thick wrapper file.

Value

An object of class hbmfit.

See Also

register_hbsae_model, list_hbsae_models, hbm

Examples

library(hbsaems)
data("data_lnln")
# Equivalent to hbm_lnln(...)
fit <- hbm_flex(
  family = "lognormal",
  response   = "y_obs",
  auxiliary  = c("x1", "x2", "x3"),
  area_var   = "district",         # area-level random effect: (1 | district)
  data       = data_lnln,
  chains = 4, iter = 2000, refresh = 0
)

---

Get Comprehensive Model Information

Description

Returns a one-page summary of the fitted model's metadata: number of observations, family, link function, formula, MCMC settings, missing-data strategy, and so on.

Usage

hbm_info(model)

Arguments

model

An hbmfit object.

Value

A named list of metadata.

Examples

library(hbsaems)
library(brms)
data("data_fhnorm")
model <- hbm(brms::bf(y ~ x1), data = data_fhnorm,
             re = ~ (1 | regency),    # area-level random effect
             chains = 4, iter = 2000, warmup = 1000,
             cores = 1, seed = 1, refresh = 0)
hbm_info(model)

---

Small Area Estimation under a Lognormal-Lognormal Model

Description

Convenience wrapper that fits a Hierarchical Bayesian Small Area Estimation model with a lognormal likelihood. Internally delegates to hbm_flex with family_key = "lognormal".

Usage

hbm_lnln(
  response,
  auxiliary = NULL,
  data,
  area_var = NULL,
  area_re_structure = c("nested", "crossed"),
  spatial_var = NULL,
  spatial_model = NULL,
  car_type = NULL,
  sar_type = NULL,
  M = NULL,
  sampling_variance = NULL,
  fixed_params = NULL,
  predictors = NULL,
  sampling_var = NULL,
  group = NULL,
  sre = NULL,
  sre_type = NULL,
  ...
)

Arguments

response	Character. Name of the response column (must be strictly positive).
auxiliary	Character vector of auxiliary (fixed-effect) variable names; corresponds to area-level covariates in SAE literature (see Rao & Molina 2015 Ch. 4).
data	A data.frame.
area_var	Optional character vector. Name(s) of a column (or columns) in data identifying the small area / domain. Length 1 adds an IID area-level random intercept (1 | area_var); length $\geq$ 2 supports hierarchical areas – see ?hbm_flex for the nested vs.\ crossed structures. Default: NULL.
area_re_structure	Either "nested" (default) or "crossed". Controls how multiple area columns combine. See ?hbm_flex.
spatial_var	Optional character. Name of a column in data identifying the spatial cluster (e.g. province). Must be supplied together with spatial_model and M. Default: NULL.
spatial_model	Optional character. Type of spatial dependence: "car" (conditional autoregressive) or "sar" (simultaneous autoregressive). Default: NULL.
car_type	Optional character. CAR sub-type passed to brms: "escar", "esicar", "icar" (intrinsic CAR; default when spatial_model = "car"), or "bym2".
sar_type	Optional character. SAR sub-type: "lag" (spatial-lag, default when spatial_model = "sar") or "error" (spatial-error).
M	Optional numeric matrix. Spatial weight matrix. Required when spatial_model is supplied.
sampling_variance	Optional character. Name of a column in data containing the known per-area sampling variance $\psi_i$ on the log scale, i.e.\ the variance of $\log(\hat y_i)$. When supplied, $\sigma_i = \sqrt{\psi_i}$ is pinned via offset, recovering the Fay–Herriot lognormal model. If your survey software produces $\mathrm{Var}(\hat y_i)$ on the original scale, convert with the delta-method approximation $\psi_i \approx \mathrm{Var}(\hat y_i) / \hat y_i^2$ before passing. Default: NULL (sigma is random).
fixed_params	Optional named list pinning distributional parameters to known values. See hbm for the spec format. Allows power-user access to the same machinery used by sampling_variance.
predictors	Deprecated. Use auxiliary instead. Kept for backward compatibility; will be removed in v2.0.0.
sampling_var	Deprecated. Use sampling_variance instead. Kept for backward compatibility; will be removed in v2.0.0.
group	Deprecated. Use area_var instead.
sre	Deprecated. Use spatial_var instead.
sre_type	Deprecated. Use spatial_model instead.
...	Additional arguments forwarded to hbm_flex (e.g.\ prior_type, nonlinear, handle_missing, sampler controls such as chains, iter, cores, seed).

Details

The response $y_i$ in area $i$ is assumed to follow

$y_i \mid \theta_i \sim \mathrm{Lognormal}(\theta_i, \sigma^2),$

with the log-mean linked to auxiliary variables via

$\log(\theta_i) = x_i^\top \boldsymbol{\beta} + u_i, \quad u_i \sim \mathcal{N}(0, \sigma_v^2).$

When the user supplies sampling_variance = "psi_i" (the column name of the known per-area sampling variance $\psi_i$), $\sigma_i = \sqrt{\psi_i}$ is pinned for each area via an offset. This recovers the Fay-Herriot-style lognormal model in which residual variability is fully determined by the survey design.

Value

An object of class hbmfit.

Conflict policy

When the residual standard deviation $\sigma_i = \sqrt{\psi_i}$ is pinned via sampling_variance (or via fixed_params$sigma), the function refuses any additional specification that would also set $\sigma$. Specifically, all of the following are rejected with an informative error at construction time:

• sampling_variance and fixed_params$sigma.

• sampling_variance and a user prior on class = "sigma".

• sampling_variance and a stanvars sampling statement involving sigma.

• auxiliary and the deprecated predictors in the same call.

See Also

hbm_flex, hbm

Examples

library(hbsaems)
data("data_lnln")

# -- 1. Standard lognormal SAE with area random effect ----------------------
fit1 <- hbm_lnln(
  response   = "y_obs",
  auxiliary  = c("x1", "x2", "x3"),
  area_var   = "district",
  data       = data_lnln,
  chains = 4, iter = 2000, warmup = 1000, refresh = 0
)

# -- 2. Fay-Herriot style with known sampling variance ----------------------
#     (assumes psi_i column is available)
fit2 <- hbm_lnln(
  response     = "y_obs",
  auxiliary    = c("x1", "x2", "x3"),
  area_var     = "district",
  sampling_variance = "psi_i",
  data         = data_lnln,
  chains = 4, iter = 2000, warmup = 1000, refresh = 0
)

---

Nonlinear-Term Configuration for HBSAE Models

Description

Bundles the nonlinear-term arguments of hbm into a single named list. Same opt-in pattern as hbm_control.

Usage

hbm_nonlinear(
  terms,
  type = c("spline", "gp"),
  k = -1L,
  spline_bs = "tp",
  gp_cov = "exp_quad",
  gp_c = NULL,
  gp_scale = NULL
)

Arguments

terms	Character vector of predictor names to be wrapped in s() or gp().
type	Character. "spline" (default) or "gp".
k	Integer. Basis dimension. For splines: passed to mgcv::s(..., k = ...); -1L (default) lets mgcv choose. For GP: passed to brms::gp(..., k = ...) for the Hilbert-space approximate GP (Riutort-Mayol et al.\ 2020); NA = exact GP (slow, not recommended for $n > 100$).
spline_bs	Character. Spline basis type ("tp", "cr", "cs", "ps", ...). Default "tp".
gp_cov	Character. GP covariance function: "exp_quad", "matern15", "matern25", "exponential". Default "exp_quad".
gp_c	Numeric or NULL. HSGP boundary-scale factor for brms::gp(c = ...).
gp_scale	Deprecated. Use gp_c instead.

Value

A named list of arguments for hbm, with class c("hbm_config_nonlinear", "hbm_config", "list").

See Also

hbm_control, hbm_priors, hbm

Examples

# Penalised regression spline (auto-chosen basis dimension):
nl_spline <- hbm_nonlinear(c("x1", "x3"), type = "spline")

# Spline with cubic-regression basis and fixed k:
nl_cr <- hbm_nonlinear("x1", type = "spline", k = 8, spline_bs = "cr")

# Hilbert-space approximate GP with Matern 5/2 covariance (recommended):
nl_gp <- hbm_nonlinear("x2", type = "gp", k = 20, gp_cov = "matern25")

---