Hierarchical Bayesian Small Area Estimation
in ‘hbsaems’
Introduction
This vignette demonstrates how to perform Bayesian modeling
for Small Area Estimation (SAE) using the tools provided in the
hbsaems package. It follows the
Bayesian workflow, a structured approach that includes
prior specification, prior predictive checks, model fitting using MCMC,
convergence diagnostics, model evaluation, and prediction.
The modeling process is centered around the hbsaems
package, which provides the following core functions to support each
stage of the workflow:
hbpc() — for prior predictive
checking, allowing users to evaluate whether the chosen priors
generate plausible data before observing the actual data.hbm() — for fitting
hierarchical Bayesian models, including fixed and random
effects, spatial structures, and handling of missing data. The
hbsaems package includes two key functions
for modeling: hbm() and
hbm_model().
hbm() offers greater flexibility in
writing formulas but requires more familiarity with model specification
and structure.
hbm_model(), on the other hand, is
more user-friendly, designed for those who may not be familiar with
writing complex formulas, offering a simpler interface for model
specification.
hbcc() — for convergence
diagnostics, assessing whether the MCMC algorithm has properly
converged using R-hat statistics, trace plots, and effective sample
sizes.hbmc() — for model evaluation
and comparison, including information criteria such as WAIC and
LOO.hbsae() — for generating
predictions at the small area level, incorporating posterior
uncertainty.
By integrating these functions, this vignette provides a complete
guide to modeling within a Bayesian framework, specifically tailored for
SAE problems. The goal is not only to fit a model, but to ensure it is
well-calibrated, interpretable, and useful for prediction.
Small Area Estimation (SAE) is a statistical
approach used to produce reliable estimates for small geographic areas
or subpopulations where direct survey estimation is limited by small
sample sizes. SAE models combine survey data with auxiliary information
to enhance precision and reduce the variability of estimates. This
vignette provides an overview of the SAE methodology, including model
formulations and implementation using the hbsaems package. The
area-level model is one of the common SAE approaches, incorporating
auxiliary data at the area level to improve estimation accuracy.
Bayesian SAE methods provide a flexible framework for incorporating
prior information and quantifying uncertainty. Bayesian models can be
specified using Markov Chain Monte Carlo (MCMC)
methods. Bayesian methods are often used to estimate the area-level
model parameters.
Model Overview
The hbm function fits a hierarchical Bayesian model
using the brms package. It supports various modeling
options, including:
- Fixed effects for covariates.
- Random effects for spatial or hierarchical structures.
- Spatial random effects using CAR/SAR structures.
- Handling missing data through deletion, model-based, or multiple
imputation.
- Allows user-defined priors for model flexibility.
- Provides model diagnostics, including convergence checks and model
comparison.
- Outputs include fitted model objects (
hbmfit),
diagnostic summaries, and predictions for small areas.
The function is designed for Small Area Estimation (SAE)
applications, providing flexible modeling capabilities for analyzing
complex hierarchical data.
Simulated Data Example
library(hbsaems)
#> Loading required package: brms
#> Loading required package: Rcpp
#> Loading 'brms' package (version 2.22.0). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
#>
#> Attaching package: 'brms'
#> The following object is masked from 'package:stats':
#>
#> ar
data <- data_fhnorm
head(data)
#> y D x1 x2 x3 theta_true u group
#> 1 4.877812 0.2951779 3.879049 7.868780 23.79524 4.352293 -0.5057526 1
#> 2 3.429696 0.1980448 4.539645 10.770651 20.24965 3.477883 -0.5322315 2
#> 3 4.595225 0.1778341 8.117417 9.259924 13.93942 4.889751 -0.6636471 3
#> 4 3.494359 0.1835110 5.141017 8.957372 17.17278 3.612569 -0.7442393 4
#> 5 4.359012 0.1185852 5.258575 7.145144 13.34264 3.975170 -0.3091185 5
#> 6 4.828993 0.1323618 8.430130 9.864917 13.09501 4.628878 0.2341790 6
#> sre
#> 1 1
#> 2 1
#> 3 1
#> 4 1
#> 5 1
#> 6 1
The data_fhn dataset provides simulated data suitable
for demonstrating the classical Fay-Herriot model under the Gaussian
likelihood. It contains information for 70 small areas, including:
y: the direct estimate for each area, assumed to
follow a Normal distribution with known sampling variance
D: the known sampling variance for each area’s
direct estimate
x1, x2, x3: auxiliary
variables that may explain variation in the outcome
theta_true: the true underlying value of the
area-specific parameter (latent mean)
u: the area-level random effect drawn from a normal
distribution
This dataset is useful for testing and illustrating small area
estimation techniques, particularly in settings where the direct
estimates are continuous and modeled using a Gaussian distribution.
Initial Model
We begin by specifying the model with informative priors on the
coefficients and the intercept. The sample_prior = "only"
argument ensures that the model ignores the data and simulates
predictions based on the priors alone. This is particularly useful for
performing a prior predictive check, which involves
generating data purely from the prior distributions to evaluate whether
the priors lead to plausible values of the outcome variable.
model_prior_pred <- hbm(
formula = bf(y ~ x1 + x2 + x3),
data = data,
hb_sampling = "gaussian",
hb_link = "log",
chains = 4,
iter = 500,
warmup = 250,
sample_prior = "only",
prior = c(
prior("normal(1, 0.2)", class = "Intercept"),
prior("normal(0, 0.1)", class = "b"),
prior("exponential(5)", class = "sd")
)
)
#> Compiling Stan program...
#> Start sampling
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 1.4e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 500 [ 0%] (Warmup)
#> Chain 1: Iteration: 50 / 500 [ 10%] (Warmup)
#> Chain 1: Iteration: 100 / 500 [ 20%] (Warmup)
#> Chain 1: Iteration: 150 / 500 [ 30%] (Warmup)
#> Chain 1: Iteration: 200 / 500 [ 40%] (Warmup)
#> Chain 1: Iteration: 250 / 500 [ 50%] (Warmup)
#> Chain 1: Iteration: 251 / 500 [ 50%] (Sampling)
#> Chain 1: Iteration: 300 / 500 [ 60%] (Sampling)
#> Chain 1: Iteration: 350 / 500 [ 70%] (Sampling)
#> Chain 1: Iteration: 400 / 500 [ 80%] (Sampling)
#> Chain 1: Iteration: 450 / 500 [ 90%] (Sampling)
#> Chain 1: Iteration: 500 / 500 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 0.029 seconds (Warm-up)
#> Chain 1: 0.01 seconds (Sampling)
#> Chain 1: 0.039 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 9e-06 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 500 [ 0%] (Warmup)
#> Chain 2: Iteration: 50 / 500 [ 10%] (Warmup)
#> Chain 2: Iteration: 100 / 500 [ 20%] (Warmup)
#> Chain 2: Iteration: 150 / 500 [ 30%] (Warmup)
#> Chain 2: Iteration: 200 / 500 [ 40%] (Warmup)
#> Chain 2: Iteration: 250 / 500 [ 50%] (Warmup)
#> Chain 2: Iteration: 251 / 500 [ 50%] (Sampling)
#> Chain 2: Iteration: 300 / 500 [ 60%] (Sampling)
#> Chain 2: Iteration: 350 / 500 [ 70%] (Sampling)
#> Chain 2: Iteration: 400 / 500 [ 80%] (Sampling)
#> Chain 2: Iteration: 450 / 500 [ 90%] (Sampling)
#> Chain 2: Iteration: 500 / 500 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 0.03 seconds (Warm-up)
#> Chain 2: 0.017 seconds (Sampling)
#> Chain 2: 0.047 seconds (Total)
#> Chain 2:
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
#> Chain 3:
#> Chain 3: Gradient evaluation took 6e-06 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3:
#> Chain 3:
#> Chain 3: Iteration: 1 / 500 [ 0%] (Warmup)
#> Chain 3: Iteration: 50 / 500 [ 10%] (Warmup)
#> Chain 3: Iteration: 100 / 500 [ 20%] (Warmup)
#> Chain 3: Iteration: 150 / 500 [ 30%] (Warmup)
#> Chain 3: Iteration: 200 / 500 [ 40%] (Warmup)
#> Chain 3: Iteration: 250 / 500 [ 50%] (Warmup)
#> Chain 3: Iteration: 251 / 500 [ 50%] (Sampling)
#> Chain 3: Iteration: 300 / 500 [ 60%] (Sampling)
#> Chain 3: Iteration: 350 / 500 [ 70%] (Sampling)
#> Chain 3: Iteration: 400 / 500 [ 80%] (Sampling)
#> Chain 3: Iteration: 450 / 500 [ 90%] (Sampling)
#> Chain 3: Iteration: 500 / 500 [100%] (Sampling)
#> Chain 3:
#> Chain 3: Elapsed Time: 0.029 seconds (Warm-up)
#> Chain 3: 0.017 seconds (Sampling)
#> Chain 3: 0.046 seconds (Total)
#> Chain 3:
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
#> Chain 4:
#> Chain 4: Gradient evaluation took 6e-06 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4:
#> Chain 4:
#> Chain 4: Iteration: 1 / 500 [ 0%] (Warmup)
#> Chain 4: Iteration: 50 / 500 [ 10%] (Warmup)
#> Chain 4: Iteration: 100 / 500 [ 20%] (Warmup)
#> Chain 4: Iteration: 150 / 500 [ 30%] (Warmup)
#> Chain 4: Iteration: 200 / 500 [ 40%] (Warmup)
#> Chain 4: Iteration: 250 / 500 [ 50%] (Warmup)
#> Chain 4: Iteration: 251 / 500 [ 50%] (Sampling)
#> Chain 4: Iteration: 300 / 500 [ 60%] (Sampling)
#> Chain 4: Iteration: 350 / 500 [ 70%] (Sampling)
#> Chain 4: Iteration: 400 / 500 [ 80%] (Sampling)
#> Chain 4: Iteration: 450 / 500 [ 90%] (Sampling)
#> Chain 4: Iteration: 500 / 500 [100%] (Sampling)
#> Chain 4:
#> Chain 4: Elapsed Time: 0.029 seconds (Warm-up)
#> Chain 4: 0.017 seconds (Sampling)
#> Chain 4: 0.046 seconds (Total)
#> Chain 4:
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess
summary(model_prior_pred)
#>
#> ======= Hierarchical Bayesian Model Summary (from hbmfit) =======
#>
#> --- Full brms Model Summary ---
#> Family: gaussian
#> Links: mu = log; sigma = identity
#> Formula: y ~ x1 + x2 + x3 + (1 | group)
#> Data: data (Number of observations: 100)
#> Draws: 4 chains, each with iter = 500; warmup = 250; thin = 1;
#> total post-warmup draws = 1000
#>
#> Multilevel Hyperparameters:
#> ~group (Number of levels: 100)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.20 0.21 0.00 0.80 1.01 546 206
#>
#> Regression Coefficients:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept 1.04 1.83 -2.58 4.72 1.00 2211 793
#> x1 0.00 0.10 -0.20 0.19 1.01 1937 716
#> x2 -0.00 0.11 -0.20 0.21 1.01 2546 758
#> x3 -0.00 0.10 -0.19 0.19 1.01 2091 623
#>
#> Further Distributional Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 2.65 2.55 0.13 9.61 1.00 1125 692
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).
#>
#>
#> ===== Detailed Group-Level Effects (Random Effects) Summary =====
#>
#> --- Grouping Factor: group ---
#> The table below (in the next console output) contains group-level estimates
#> for the factor 'group', with dimensions 1 rows x 7 columns.
#>
#> Columns: Estimate, Est.Error, l-95% CI, u-95% CI, Rhat, Bulk_ESS, Tail_ESS. Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> sd(Intercept) 0.2046175 0.2058147 0.004169071 0.8039563 1.00677 545.6605
#> Tail_ESS
#> sd(Intercept) 206.1716
#>
#>
#> =========== Missing Data Handling (specified in hbm) ============
#> No specific missing data handling method was applied through the hbm function.
#> (brms internal defaults for NA handling in predictors/response may have applied if NAs were present and not pre-processed).
Prior Predictive Check
Before fitting a Bayesian model to the data, it is important to
evaluate whether the chosen priors lead to sensible predictions. This
process is called a prior predictive check.
Prior predictive checks help to:
- Identify whether the priors are overly vague or unrealistically
informative.
- Ensure that the prior assumptions do not produce unreasonable
predictions before any actual data are used.
- Avoid unintentionally strong or misleading prior beliefs that could
affect posterior inference.
This step is especially important in Bayesian modeling to build
confidence that the model structure and prior choices are coherent
before observing any data.
result_hbpc <- hbpc(model_prior_pred)
#> Using 'y' as the response variable for pp_check.
summary(result_hbpc)
#>
#> ================= Prior Predictive Check Summary =================
#>
#> -------------- Priors Used in the Model --------------
#> prior class coef group resp dpar nlpar lb ub
#> normal(1, 0.2) Intercept
#> normal(0, 0.1) b
#> normal(0, 0.1) b x1
#> normal(0, 0.1) b x2
#> normal(0, 0.1) b x3
#> exponential(5) sd 0
#> exponential(5) sd group 0
#> exponential(5) sd Intercept group 0
#> student_t(3, 0, 2.5) sigma 0
#> source
#> user
#> user
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> user
#> (vectorized)
#> (vectorized)
#> default
#>
#> ----- Summary of Parameter Draws from Prior Model -----
#> Family: gaussian
#> Links: mu = log; sigma = identity
#> Formula: y ~ x1 + x2 + x3 + (1 | group)
#> Data: data (Number of observations: 100)
#> Draws: 4 chains, each with iter = 500; warmup = 250; thin = 1;
#> total post-warmup draws = 1000
#>
#> Multilevel Hyperparameters:
#> ~group (Number of levels: 100)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.20 0.21 0.00 0.80 1.01 546 206
#>
#> Regression Coefficients:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept 1.04 1.83 -2.58 4.72 1.00 2211 793
#> x1 0.00 0.10 -0.20 0.19 1.01 1937 716
#> x2 -0.00 0.11 -0.20 0.21 1.01 2546 758
#> x3 -0.00 0.10 -0.19 0.19 1.01 2091 623
#>
#> Further Distributional Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 2.65 2.55 0.13 9.61 1.00 1125 692
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).
#>
#> --------------- Priors Predictive Plot ----------------
#>
#> Note: The prior predictive plot is stored in '$prior_predictive_plot'.
#> This plot compares observed data (if applicable) with predictions generated solely from the priors.
result_hbpc$prior_predictive_plot
#> Warning in transformation$transform(x): NaNs produced
#> Warning in ggplot2::scale_x_log10(): log-10 transformation introduced infinite
#> values.
#> Warning in transformation$transform(x): NaNs produced
#> Warning in ggplot2::scale_x_log10(): log-10 transformation introduced infinite
#> values.
#> Warning: Removed 571 rows containing non-finite outside the scale range
#> (`stat_density()`).
#> Warning: Removed 1 row containing non-finite outside the scale range
#> (`stat_density()`).
The prior predictive plot indicates that the priors used in the model
are reasonable and appropriate. Here’s why:
The prior predictive distributions (shown in
blue) cover a wide but plausible range of values for the response
variable (medv), which aligns well with the distribution of
the observed data (shown in black).
The shape of the simulated predictions is consistent with
the general pattern of the observed outcome, without being
overly concentrated or too diffuse.
There are no extreme or unrealistic values
produced by the priors, suggesting that the priors are informative
enough to guide the model without being too restrictive.
Fitting the Model
After prior predictive checks confirm the priors are suitable, we
proceed to fit the model using sample_prior = "no". Why
sample_prior = "no"?
We already conducted a prior predictive check (previous
explanation), sample_prior = "no" tells brms
to:
This option is standard for final model fitting after priors have
been validated.
1. Basic Model
model <- hbm(
formula = bf(y ~ x1 + x2 + x3),
data = data,
hb_sampling = "gaussian",
hb_link = "log",
chains = 4,
iter = 500,
warmup = 250,
sample_prior = "no",
prior = c(
prior("normal(1, 0.2)", class = "Intercept"),
prior("normal(0, 0.1)", class = "b"),
prior("exponential(5)", class = "sd")
)
)
#> Compiling Stan program...
#> Start sampling
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 2.8e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 500 [ 0%] (Warmup)
#> Chain 1: Iteration: 50 / 500 [ 10%] (Warmup)
#> Chain 1: Iteration: 100 / 500 [ 20%] (Warmup)
#> Chain 1: Iteration: 150 / 500 [ 30%] (Warmup)
#> Chain 1: Iteration: 200 / 500 [ 40%] (Warmup)
#> Chain 1: Iteration: 250 / 500 [ 50%] (Warmup)
#> Chain 1: Iteration: 251 / 500 [ 50%] (Sampling)
#> Chain 1: Iteration: 300 / 500 [ 60%] (Sampling)
#> Chain 1: Iteration: 350 / 500 [ 70%] (Sampling)
#> Chain 1: Iteration: 400 / 500 [ 80%] (Sampling)
#> Chain 1: Iteration: 450 / 500 [ 90%] (Sampling)
#> Chain 1: Iteration: 500 / 500 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 0.266 seconds (Warm-up)
#> Chain 1: 0.038 seconds (Sampling)
#> Chain 1: 0.304 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 1.7e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 500 [ 0%] (Warmup)
#> Chain 2: Iteration: 50 / 500 [ 10%] (Warmup)
#> Chain 2: Iteration: 100 / 500 [ 20%] (Warmup)
#> Chain 2: Iteration: 150 / 500 [ 30%] (Warmup)
#> Chain 2: Iteration: 200 / 500 [ 40%] (Warmup)
#> Chain 2: Iteration: 250 / 500 [ 50%] (Warmup)
#> Chain 2: Iteration: 251 / 500 [ 50%] (Sampling)
#> Chain 2: Iteration: 300 / 500 [ 60%] (Sampling)
#> Chain 2: Iteration: 350 / 500 [ 70%] (Sampling)
#> Chain 2: Iteration: 400 / 500 [ 80%] (Sampling)
#> Chain 2: Iteration: 450 / 500 [ 90%] (Sampling)
#> Chain 2: Iteration: 500 / 500 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 0.425 seconds (Warm-up)
#> Chain 2: 0.078 seconds (Sampling)
#> Chain 2: 0.503 seconds (Total)
#> Chain 2:
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
#> Chain 3:
#> Chain 3: Gradient evaluation took 1.5e-05 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3:
#> Chain 3:
#> Chain 3: Iteration: 1 / 500 [ 0%] (Warmup)
#> Chain 3: Iteration: 50 / 500 [ 10%] (Warmup)
#> Chain 3: Iteration: 100 / 500 [ 20%] (Warmup)
#> Chain 3: Iteration: 150 / 500 [ 30%] (Warmup)
#> Chain 3: Iteration: 200 / 500 [ 40%] (Warmup)
#> Chain 3: Iteration: 250 / 500 [ 50%] (Warmup)
#> Chain 3: Iteration: 251 / 500 [ 50%] (Sampling)
#> Chain 3: Iteration: 300 / 500 [ 60%] (Sampling)
#> Chain 3: Iteration: 350 / 500 [ 70%] (Sampling)
#> Chain 3: Iteration: 400 / 500 [ 80%] (Sampling)
#> Chain 3: Iteration: 450 / 500 [ 90%] (Sampling)
#> Chain 3: Iteration: 500 / 500 [100%] (Sampling)
#> Chain 3:
#> Chain 3: Elapsed Time: 0.456 seconds (Warm-up)
#> Chain 3: 0.043 seconds (Sampling)
#> Chain 3: 0.499 seconds (Total)
#> Chain 3:
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
#> Chain 4:
#> Chain 4: Gradient evaluation took 1.6e-05 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.16 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4:
#> Chain 4:
#> Chain 4: Iteration: 1 / 500 [ 0%] (Warmup)
#> Chain 4: Iteration: 50 / 500 [ 10%] (Warmup)
#> Chain 4: Iteration: 100 / 500 [ 20%] (Warmup)
#> Chain 4: Iteration: 150 / 500 [ 30%] (Warmup)
#> Chain 4: Iteration: 200 / 500 [ 40%] (Warmup)
#> Chain 4: Iteration: 250 / 500 [ 50%] (Warmup)
#> Chain 4: Iteration: 251 / 500 [ 50%] (Sampling)
#> Chain 4: Iteration: 300 / 500 [ 60%] (Sampling)
#> Chain 4: Iteration: 350 / 500 [ 70%] (Sampling)
#> Chain 4: Iteration: 400 / 500 [ 80%] (Sampling)
#> Chain 4: Iteration: 450 / 500 [ 90%] (Sampling)
#> Chain 4: Iteration: 500 / 500 [100%] (Sampling)
#> Chain 4:
#> Chain 4: Elapsed Time: 0.439 seconds (Warm-up)
#> Chain 4: 0.075 seconds (Sampling)
#> Chain 4: 0.514 seconds (Total)
#> Chain 4:
#> Warning: There were 28 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: The largest R-hat is 1.09, indicating chains have not mixed.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#r-hat
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess
summary(model)
#>
#> ======= Hierarchical Bayesian Model Summary (from hbmfit) =======
#>
#> --- Full brms Model Summary ---
#> Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
#> careful when analysing the results! We recommend running more iterations and/or
#> setting stronger priors.
#> Warning: There were 28 divergent transitions after warmup. Increasing
#> adapt_delta above 0.8 may help. See
#> http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> Family: gaussian
#> Links: mu = log; sigma = identity
#> Formula: y ~ x1 + x2 + x3 + (1 | group)
#> Data: data (Number of observations: 100)
#> Draws: 4 chains, each with iter = 500; warmup = 250; thin = 1;
#> total post-warmup draws = 1000
#>
#> Multilevel Hyperparameters:
#> ~group (Number of levels: 100)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.06 0.04 0.00 0.15 1.06 75 27
#>
#> Regression Coefficients:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept 0.86 0.18 0.47 1.20 1.06 49 23
#> x1 0.11 0.02 0.08 0.15 1.07 42 26
#> x2 -0.09 0.01 -0.11 -0.07 1.01 1018 528
#> x3 0.04 0.01 0.03 0.06 1.05 82 157
#>
#> Further Distributional Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 0.95 0.08 0.79 1.11 1.06 61 45
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).
#> Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
#> careful when analysing the results! We recommend running more iterations and/or
#> setting stronger priors.
#> Warning: There were 28 divergent transitions after warmup. Increasing
#> adapt_delta above 0.8 may help. See
#> http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#>
#>
#> ===== Detailed Group-Level Effects (Random Effects) Summary =====
#>
#> --- Grouping Factor: group ---
#> The table below (in the next console output) contains group-level estimates
#> for the factor 'group', with dimensions 1 rows x 7 columns.
#>
#> Columns: Estimate, Est.Error, l-95% CI, u-95% CI, Rhat, Bulk_ESS, Tail_ESS. Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> sd(Intercept) 0.06188201 0.04405645 0.001418994 0.1549819 1.05605 74.96048
#> Tail_ESS
#> sd(Intercept) 27.12551
#>
#>
#> =========== Missing Data Handling (specified in hbm) ============
#> No specific missing data handling method was applied through the hbm function.
#> (brms internal defaults for NA handling in predictors/response may have applied if NAs were present and not pre-processed).
By default, if we do not explicitly define a random effect, the model
will still generate one based on the natural random variations between
individual records (rows) in the dataset. However, we can also
explicitly define a random effect to account for variations at a
specific hierarchical level, such as neighborhoods or residential
blocks. We can use parameter re.
Note: Some chunk is not executed in this vignette. It is provided as
an example and will produce output similar to the model results.
model_with_defined_re <- hbm(
formula = bf(y ~ x1 + x2 + x3),
data = data,
hb_sampling = "gaussian",
hb_link = "log",
re = ~(1|group),
chains = 4,
iter = 4000,
warmup = 2000,
sample_prior = "no",
prior = c(
prior("normal(1, 0.2)", class = "Intercept"),
prior("normal(0, 0.1)", class = "b"),
prior("exponential(5)", class = "sd")
)
)
summary(model_with_defined_re)
2. Model With Missing Data
The hbm function supports three strategies for handling missing
data:
“deleted”: Removes rows with missing values.
“multiple”: Performs multiple imputation using mice.
“model”: Uses the mi() function to model missingness directly
(not available for discrete outcomes).
data_missing <- data
data_missing$y[3:5] <- NA
Handling missing data by deleted (Only if missing in
response)
model_deleted <- hbm(
formula = bf(y ~ x1 + x2 + x3),
hb_sampling = "gaussian",
hb_link = "log",
re = ~(1|group),
data = data_missing,
handle_missing = "deleted",
chains = 4,
iter = 4000,
warmup = 2000,
sample_prior = "no",
prior = c(
prior("normal(1, 0.2)", class = "Intercept"),
prior("normal(0, 0.1)", class = "b"),
prior("exponential(5)", class = "sd")
)
)
Handling missing data before model
fitting using multiple imputation
model_multiple <- hbm(
formula = bf(y ~ x1 + x2 + x3),
hb_sampling = "gaussian",
hb_link = "log",
re = ~(1|group),
data = data_missing,
handle_missing = "multiple",
chains = 4,
iter = 4000,
warmup = 2000,
sample_prior = "no",
prior = c(
prior("normal(1, 0.2)", class = "Intercept"),
prior("normal(0, 0.1)", class = "b"),
prior("exponential(5)", class = "sd")
)
)
Handle missing data during model fitting using
mi()
If we want to handle missing data in the model, we must define the
formula correctly. In this case, we need to specify missing data
indicators using the mi() function.
data_missing$y[3:5] <- NA
data_missing$x1[6:7] <- NA
model_during_model <- hbm(
formula = bf(y | mi() ~ mi(x1) + x2 + x3) + bf(x1 | mi() ~ x2 + x3),
hb_sampling = "gaussian",
hb_link = "log",
re = ~(1|group),
data = data_missing,
handle_missing = "model",
chains = 4,
iter = 4000,
warmup = 2000,
sample_prior = "no",
prior = c(
prior("normal(1, 0.2)", class = "Intercept", resp = "y"),
prior("normal(0, 0.1)", class = "b", resp = "y"),
prior("exponential(5)", class = "sd", resp = "y"),
prior("normal(1, 0.2)", class = "Intercept", resp = "x1"),
prior("normal(0, 0.1)", class = "b", resp = "x1"),
prior("exponential(5)", class = "sd", resp = "x1")
)
)
summary(model_during_model)
3. Model With Spatial Effect
M <- adjacency_matrix_car
model_spatial <- hbm(
formula = bf(y ~ x1 + x2 + x3),
data = data,
hb_sampling = "gaussian",
hb_link = "log",
sre = "sre",
sre_type = "car",
car_type = "icar",
M = M,
chains = 4,
iter = 4000,
warmup = 2000,
sample_prior = "no",
prior = c(
prior("normal(1, 0.2)", class = "Intercept"),
prior("normal(0, 0.1)", class = "b"),
prior("exponential(5)", class = "sd")
)
)
Model Diagnostics
The hbcc function is designed to evaluate the
convergence and quality of a Bayesian hierarchical model. It performs
several diagnostic tests and generates various plots to assess Markov
Chain Monte Carlo (MCMC) performance.
result_hbcc <- hbcc(model)
#> Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
#> careful when analysing the results! We recommend running more iterations and/or
#> setting stronger priors.
#> Warning: There were 28 divergent transitions after warmup. Increasing
#> adapt_delta above 0.8 may help. See
#> http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
#> careful when analysing the results! We recommend running more iterations and/or
#> setting stronger priors.
#> Warning: There were 28 divergent transitions after warmup. Increasing
#> adapt_delta above 0.8 may help. See
#> http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
#> careful when analysing the results! We recommend running more iterations and/or
#> setting stronger priors.
#> Warning: There were 28 divergent transitions after warmup. Increasing
#> adapt_delta above 0.8 may help. See
#> http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
summary(result_hbcc)
#>
#> =========== Model Convergence Diagnostics Summary ===========
#>
#> ----- R-hat and Effective Sample Size (ESS) -----
#> The following table (shown in the next console output) summarizes convergence diagnosticsincluding R-hat, Bulk ESS, and Tail ESS for each parameter.
#> Table dimensions: 5 rows x 3 columns.
#> Columns: Rhat, Bulk_ESS, Tail_ESS.
#>
#> Rhat Bulk_ESS Tail_ESS
#> Intercept 1.063166 49.29274 22.99907
#> x1 1.070539 42.27035 26.28425
#> x2 1.007382 1018.47012 528.38172
#> x3 1.046599 81.96881 157.41408
#> sd(Intercept) 1.056050 74.96048 27.12551
#>
#> ----------------- Geweke Diagnostic -----------------
#> b_Intercept b_x1 b_x2 b_x3
#> -2.1094770 1.9336186 -0.4812567 1.9520558
#> sd_group__Intercept sigma
#> 1.3397922 -1.6356583
#> Note: Displaying Z-scores. Values outside approx. +/- 2 may indicate non-convergence.
#>
#> -------------- Raftery-Lewis Diagnostic --------------
#> [1] "Error" "3746"
#>
#> ----------- Heidelberger-Welch Diagnostic -----------
#> stest start pvalue htest mean halfwidth
#> b_Intercept 1 101 0.4114651 1 0.88742777 0.0090560432
#> b_x1 1 101 0.9017306 1 0.10846705 0.0007712603
#> b_x2 1 1 0.8844321 1 -0.08933621 0.0005987697
#> b_x3 1 101 0.2006029 1 0.04356016 0.0003787719
#> sd_group__Intercept 1 101 0.2207911 1 0.05709131 0.0047540289
#> sigma 1 201 0.6337887 1 0.95793269 0.0068149924
#> Note: First column indicates stationarity test result (1=pass, 0=fail).
#>
#> ------------------- Diagnostic Plot --------------------
#>
#> Diagnostic plots are stored in the '$plots' element and can be printed separately (e.g., object$plots$trace).
The update_hbm() function allows you to continue
sampling from an existing fitted model by increasing the number of
iterations. This is particularly useful when initial model fitting did
not achieve convergence, for example due to insufficient iterations or
complex posterior geometry.
When update_hbm() is called with additional iterations,
the sampler resumes from the previous fit, effectively increasing the
total number of posterior samples. This helps improve convergence
diagnostics such as Rhat, effective sample size, and chain mixing.
model <- update_hbm(model, iter = 8000)
result_hbcc <- hbcc(model)
#> Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
#> careful when analysing the results! We recommend running more iterations and/or
#> setting stronger priors.
#> Warning: There were 28 divergent transitions after warmup. Increasing
#> adapt_delta above 0.8 may help. See
#> http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
#> careful when analysing the results! We recommend running more iterations and/or
#> setting stronger priors.
#> Warning: There were 28 divergent transitions after warmup. Increasing
#> adapt_delta above 0.8 may help. See
#> http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
#> careful when analysing the results! We recommend running more iterations and/or
#> setting stronger priors.
#> Warning: There were 28 divergent transitions after warmup. Increasing
#> adapt_delta above 0.8 may help. See
#> http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
summary(result_hbcc)
#>
#> =========== Model Convergence Diagnostics Summary ===========
#>
#> ----- R-hat and Effective Sample Size (ESS) -----
#> The following table (shown in the next console output) summarizes convergence diagnosticsincluding R-hat, Bulk ESS, and Tail ESS for each parameter.
#> Table dimensions: 5 rows x 3 columns.
#> Columns: Rhat, Bulk_ESS, Tail_ESS.
#>
#> Rhat Bulk_ESS Tail_ESS
#> Intercept 1.063166 49.29274 22.99907
#> x1 1.070539 42.27035 26.28425
#> x2 1.007382 1018.47012 528.38172
#> x3 1.046599 81.96881 157.41408
#> sd(Intercept) 1.056050 74.96048 27.12551
#>
#> ----------------- Geweke Diagnostic -----------------
#> b_Intercept b_x1 b_x2 b_x3
#> -2.1094770 1.9336186 -0.4812567 1.9520558
#> sd_group__Intercept sigma
#> 1.3397922 -1.6356583
#> Note: Displaying Z-scores. Values outside approx. +/- 2 may indicate non-convergence.
#>
#> -------------- Raftery-Lewis Diagnostic --------------
#> [1] "Error" "3746"
#>
#> ----------- Heidelberger-Welch Diagnostic -----------
#> stest start pvalue htest mean halfwidth
#> b_Intercept 1 101 0.4114651 1 0.88742777 0.0090560432
#> b_x1 1 101 0.9017306 1 0.10846705 0.0007712603
#> b_x2 1 1 0.8844321 1 -0.08933621 0.0005987697
#> b_x3 1 101 0.2006029 1 0.04356016 0.0003787719
#> sd_group__Intercept 1 101 0.2207911 1 0.05709131 0.0047540289
#> sigma 1 201 0.6337887 1 0.95793269 0.0068149924
#> Note: First column indicates stationarity test result (1=pass, 0=fail).
#>
#> ------------------- Diagnostic Plot --------------------
#>
#> Diagnostic plots are stored in the '$plots' element and can be printed separately (e.g., object$plots$trace).
The trace plot
(result_hbcc$plots$trace) shows the sampled values of each
parameter across MCMC iterations for all chains. A well-mixed and
stationary trace (with overlapping chains) indicates good convergence
and suggests that the sampler has thoroughly explored the posterior
distribution.
The density plot
(result_hbcc$plots$dens) displays the estimated posterior
distributions of the model parameters. When the distributions from
multiple chains align closely, it supports the conclusion that the
chains have converged to the same target distribution.
The autocorrelation plot (ACF)
(result_hbcc$plots$acf) visualizes the correlation between
samples at different lags. High autocorrelation can indicate inefficient
sampling, while rapid decay in autocorrelation across lags suggests that
the chains are generating nearly independent samples.
The NUTS energy plot
(result_hbcc$plots$nuts_energy) is specific to models
fitted using the No-U-Turn Sampler (NUTS). This plot examines the
distribution of the Hamiltonian energy and its changes between
iterations. A well-behaved energy plot indicates stable dynamics and
efficient exploration of the parameter space.
result_hbcc$plots$nuts_energy
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
The Rhat plot (result_hbcc$plots$rhat)
presents the potential scale reduction factor (R̂) for each parameter. R̂
values close to 1.00 (typically < 1.01) are a strong indicator of
convergence, showing that the between-chain and within-chain variances
are consistent.
The effective sample size (n_eff) plot
(result_hbcc$plots$neff) reports how many effectively
independent samples have been drawn for each parameter. Higher effective
sample sizes are desirable, as they indicate that the posterior
estimates are based on a substantial amount of independent information,
even if the chains themselves are autocorrelated.
Model Comparison
The hbmc function is used to compare Bayesian
hierarchical models and assess their posterior predictive performance.
It allows for model comparison, posterior predictive checks, and
additional diagnostic visualizations.
result_hbmc <- hbmc(
model = list(model),
comparison_metrics = c("loo", "waic", "bf"),
run_prior_sensitivity= TRUE,
sensitivity_vars = c("b_Intercept", "b_x1")
)
#> WarningModel 'model1': 4 high Pareto k values (> 0.7). Consider setting moment_match = TRUE.
#> Warning:
#> 7 (7.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.
summary(result_hbmc)
#>
#> =============== Model Goodness of Fit & Comparison Summary9 ===============
#>
#> =============== Primary Model Diagnostics ===============
#> Diagnostics for Primary Model model1
#>
#> ------------- Leave-One-Out (LOO) -------------
#> Estimate SE
#> elpd_loo -141.84834 7.468959
#> p_loo 11.42222 1.744477
#> looic 283.69668 14.937918
#> >> Pareto k diagnostic summary:
#> NULL
#> >> Number of high Pareto k values (>0.7): 4
#>
#>
#> -- Widely Applicable Information Criterion (WAIC) --
#> Estimate SE
#> elpd_waic -141.48803 7.410546
#> p_waic 11.06191 1.668086
#> waic 282.97607 14.821091
#>
#> ----------- Posterior Predictive Check Plot ----------
#> Posterior predictive check plot for primary model stored in
#> '$primary_model_diagnostics$pp_check_plot'.
#>
#> ------------------- Parameter Plot -------------------
#> Parameter plots for primary model stored in
#> '$primary_model_diagnostics$params_plot'.
#>
#> =============== Model Comparison Results ================
#> Model comparison metrics was not run or no results are available.
#>
#> ========== Prior Sensitivity Analysis Results ===========
#>
#> --- Prior Sensitivity Analysis Summary ---
#>
#>
#> - For model 'model1':
#>
#> > Sensitivity result:
#> Sensitivity based on cjs_dist
#> Prior selection: all priors
#> Likelihood selection: all data
#>
#> variable prior likelihood diagnosis
#> b_Intercept 0.034 0.551 -
#> b_x1 0.045 0.519 -
#>
#> > Plot:
#> Plot may appear in a separate graphics window.
The posterior predictive check plot
(result_hbmc$primary_model_diagnostics$pp_check_plot)
compares the observed data to replicated datasets generated from the
posterior predictive distribution. This visualization helps evaluate how
well the model reproduces the observed data. A good model fit is
indicated when the simulated data closely resemble the actual data in
distribution and structure.
result_hbmc$primary_model_diagnostics$pp_check_plot
The marginal posterior distributions plot
(result_hbmc$primary_model_diagnostics$params_plot)
displays the estimated distributions of the model parameters based on
the posterior samples. This plot is useful for interpreting the
uncertainty and central tendency of each parameter. Peaks in the
distribution indicate likely values, while wider distributions suggest
greater uncertainty.
result_hbmc$primary_model_diagnostics$params_plot
#> [[1]]
#>
#> [[2]]
Small Area Estimation Predictions
The hbsae() function implements Hierarchical Bayesian
Small Area Estimation (HBSAE) using the
brms package. It generates small area estimates while
accounting for uncertainty and hierarchical structure in the data. The
function calculates Mean Predictions, Relative Standard Error (RSE),
Mean Squared Error (MSE), and Root Mean Squared Error (RMSE) based on
the posterior predictive sample from the fitted Bayesian model.
result_hbsae <- hbsae(model)
summary(result_hbsae)
#>
#> =============== Summary of hbsae Prediction ===============
#>
#> >> Model Performance Metrics (based on provided data for prediction)
#> Mean RSE (%) of Predictions: 8.65
#> Mean MSE of Predictions: 0.11
#>
#> =================== Prediction Table =====================
#> Access the full data frame with $result_table
#> The table below (in the next console output) shows the model results summary.
#> Dimensions: 100 rows x 4 columns.
#> Columns: Prediction, Actual, RSE_percent, MSE.
#> Prediction Actual RSE_percent MSE
#> 1 5.102709 4.8778119 8.046214 0.16857173
#> 2 3.617815 3.4296960 8.177675 0.08752908
#> 3 4.677166 4.5952246 7.495683 0.12291023
#> 4 3.961038 3.4943592 7.140165 0.07998977
#> 5 4.073889 4.3590125 7.428945 0.09159519
#> 6 4.518972 4.8289930 8.447144 0.14571330
#> 7 3.838032 3.0599966 7.655485 0.08633017
#> 8 3.561153 3.9521435 9.298994 0.10966126
#> 9 4.304443 5.5650290 8.524342 0.13463431
#> 10 2.366728 3.5308691 8.677110 0.04217420
#> 11 5.246489 6.0562380 7.106629 0.13901602
#> 12 3.203247 3.7808892 8.030323 0.06616784
#> 13 6.492811 4.8399155 8.177926 0.28193682
#> 14 3.113419 2.9036814 7.527171 0.05492103
#> 15 2.148369 2.3764529 8.474434 0.03314663
#> 16 6.184351 6.8071839 9.108348 0.31729811
#> 17 3.348325 3.9543712 7.457229 0.06234623
#> 18 2.225650 1.2465812 10.542800 0.05505867
#> 19 3.952788 4.9474920 7.848781 0.09625241
#> 20 3.168714 3.6085892 8.756751 0.07699314
#> 21 2.249341 1.4324134 9.150211 0.04236164
#> 22 4.570818 5.2077723 8.149595 0.13875860
#> 23 3.643015 3.6752604 7.398340 0.07264247
#> 24 3.405333 3.6455150 7.895956 0.07229838
#> 25 1.630391 1.1445572 12.555379 0.04190279
#> 26 2.713533 2.8511037 9.678446 0.06897338
#> 27 3.269619 2.8919835 7.652978 0.06261168
#> 28 2.972887 4.1663995 8.531174 0.06432421
#> 29 4.056641 6.1395101 9.345597 0.14372990
#> 30 3.701698 3.8942044 8.101696 0.08994020
#> 31 3.604830 4.8580533 11.457359 0.17058414
#> 32 2.649566 1.5490382 8.366809 0.04914386
#> 33 4.216215 5.0748393 7.769368 0.10730425
#> 34 4.095208 5.3787647 9.750867 0.15945511
#> 35 6.279294 6.7605075 7.463862 0.21965894
#> 36 2.159042 -0.5118425 11.517351 0.06183402
#> 37 5.094396 4.1515834 9.807737 0.24964508
#> 38 2.896800 3.2823407 7.915392 0.05257530
#> 39 1.988919 2.3050308 9.835305 0.03826572
#> 40 3.840756 3.5818581 8.913177 0.11719212
#> 41 2.007398 1.1392960 9.832038 0.03895417
#> 42 3.136701 4.0108705 7.199847 0.05100267
#> 43 5.072484 5.9332207 8.980068 0.20749166
#> 44 6.537032 6.9054779 7.970681 0.27148886
#> 45 6.244763 5.4711920 6.859044 0.18346747
#> 46 4.026407 2.7076652 8.215362 0.10941801
#> 47 4.404316 4.5140762 7.748649 0.11646859
#> 48 1.950932 1.6188507 9.459195 0.03405592
#> 49 1.954209 1.0960787 11.375510 0.04941784
#> 50 4.964691 5.0004754 7.440176 0.13644287
#> 51 2.694261 3.3110774 9.293260 0.06269248
#> 52 2.335392 0.2881463 9.610154 0.05037097
#> 53 2.822280 2.7036679 7.978623 0.05070564
#> 54 5.923933 6.1059180 7.116244 0.17771415
#> 55 4.270969 3.9536459 7.875915 0.11315008
#> 56 4.989620 5.6606685 7.901139 0.15542267
#> 57 2.467381 2.8434503 9.089107 0.05029386
#> 58 4.727701 5.4163079 8.278457 0.15317885
#> 59 2.475503 0.2073796 9.038195 0.05005993
#> 60 3.081217 4.7841618 14.318255 0.19463673
#> 61 2.438325 1.7657814 9.258370 0.05096267
#> 62 3.607512 3.8841889 7.843497 0.08006361
#> 63 4.334188 4.5876794 6.989129 0.09176172
#> 64 1.396952 1.6275968 13.846745 0.03741611
#> 65 4.380603 4.8066269 8.492796 0.13841050
#> 66 4.225297 4.0449047 7.428866 0.09852796
#> 67 3.010980 3.6502501 8.264953 0.06192934
#> 68 2.775453 2.0869592 9.064980 0.06329966
#> 69 3.385389 4.4870445 8.699349 0.08673424
#> 70 4.730218 4.5224078 8.313433 0.15464043
#> 71 3.667125 3.8753188 7.999793 0.08606149
#> 72 2.278604 0.8468650 10.575118 0.05806414
#> 73 4.770172 5.5516865 7.166181 0.11685401
#> 74 1.242402 0.2532807 12.396102 0.02371892
#> 75 4.090563 4.8120917 7.790494 0.10155380
#> 76 4.931991 3.6745948 7.988115 0.15521479
#> 77 3.111525 2.4792518 8.053793 0.06279825
#> 78 2.306213 1.7187768 10.565164 0.05936783
#> 79 3.342215 4.1796409 7.843280 0.06871701
#> 80 3.590973 3.3917136 7.575098 0.07399475
#> 81 3.258022 3.1442226 9.142001 0.08871367
#> 82 2.907908 3.3797754 8.862334 0.06641366
#> 83 3.624350 4.3456497 7.160286 0.06734744
#> 84 4.455195 3.3624261 7.157910 0.10169649
#> 85 3.508964 4.4581308 7.774824 0.07442847
#> 86 3.911967 4.5261084 8.934598 0.12216319
#> 87 3.248348 3.7134602 8.653721 0.07901891
#> 88 4.624796 4.2361325 7.556700 0.12213766
#> 89 2.376503 1.9335444 9.433513 0.05026014
#> 90 5.007472 5.2173241 6.877004 0.11858661
#> 91 4.675099 4.1459538 7.600523 0.12626080
#> 92 4.655806 3.0236563 7.974070 0.13783193
#> 93 4.370818 6.4145220 10.218396 0.19947610
#> 94 3.270847 2.6626921 7.984881 0.06821145
#> 95 8.635957 7.8759563 7.028863 0.36846066
#> 96 1.697907 0.5480250 10.743735 0.03327656
#> 97 6.211099 5.8510226 8.707761 0.29251619
#> 98 5.212860 5.8498897 9.303803 0.23521942
#> 99 3.586353 2.4521186 8.582096 0.09473112
#> 100 4.536686 5.0931048 8.044694 0.13319762
Conclusion
The hbm function in the hbsaems package offers a comprehensive
framework for Small Area Estimation hierarchical Bayesian modeling. It
supports flexible model specifications, advanced handling of missing
data, and integrates well with diagnostic and comparison tools such as
hbcc, hbmc, and hbsae.