A typical Bayesian SAE workflow with hbsaems v1.1.0
follows seven steps, each backed by a single primary function. This
vignette walks through the canonical pipeline.
For deeper coverage of topics beyond this introductory workflow, see the articles on the package website:
https://madsyair.github.io/hbsaems/articles/
The articles cover:
data("data_fhnorm")
str(data_fhnorm, max.level = 1)
#> 'data.frame': 100 obs. of 9 variables:
#> $ y : num 6.98 9.01 6.22 10.47 10.62 ...
#> $ D : num 0.921 2.643 1.704 0.659 1.017 ...
#> $ x1 : num -0.9026 -0.7527 -0.2606 0.0759 0.2118 ...
#> $ x2 : num 0.712 -0.878 -0.607 -1.077 -0.502 ...
#> $ x3 : num 1.431 -0.318 -1.029 -1.071 -0.496 ...
#> $ theta_true: num 8.5 9.73 8.54 9.68 11.86 ...
#> $ u : num -0.8531 -0.0145 -1.2505 -0.5989 1.5845 ...
#> $ regency : chr "regency_001" "regency_002" "regency_003" "regency_004" ...
#> $ province : chr "province_01" "province_01" "province_01" "province_01" ...
#> - attr(*, "datalabel")= chr "data_fhnorm"
#> - attr(*, "var.labels")= chr [1:9] "" "" "" "" ...
head(data_fhnorm[, c("regency", "province", "y", "D", "x1", "x2", "x3")], 4)
#> regency province y D x1 x2 x3
#> 1 regency_001 province_01 6.982249 0.9213655 -0.9026345 0.7116302 1.4305563
#> 2 regency_002 province_01 9.005801 2.6432183 -0.7526745 -0.8779797 -0.3182451
#> 3 regency_003 province_01 6.216509 1.7035092 -0.2605613 -0.6066989 -1.0293040
#> 4 regency_004 province_01 10.469687 0.6586857 0.0759456 -1.0766939 -1.0711613data_fhnorm is a simulated Fay-Herriot dataset: 100
regencies nested within 5 provinces, with a known sampling variance
D per area and three covariates.
# This is the production call. Replace `chains`, `iter` with your
# usual values; the lighter settings below are a quick demonstration.
model_prior <- hbm(
formula = brms::bf(y ~ x1 + x2 + x3),
data = data_fhnorm,
re = ~ (1 | regency), # area-level random effect (Fay-Herriot)
sampling_variance = "D", # KNOWN sampling variance D_i
sample_prior = "only",
prior = c(
brms::prior(normal(0, 1), class = "b"),
brms::prior(normal(10, 5), class = "Intercept"),
brms::prior(normal(0, 2), class = "sd")
),
chains = 2, iter = 1000
)
pc <- prior_check(model_prior,
data = data_fhnorm,
response_var = "y")
plot(pc)The sampling_variance = "D" argument is the
defining feature of a Fay-Herriot model: the sampling
variance is treated as known from the design (it is not
estimated), so brms pins the observation-level to . Omitting this
argument makes the model unidentified because the residual and the area
random-effect would compete to explain the same variance, typically
producing divergent transitions and very wide credible intervals.
model <- hbm(
formula = brms::bf(y ~ x1 + x2 + x3),
data = data_fhnorm,
re = ~ (1 | regency),
sampling_variance = "D",
control = list(adapt_delta = 0.99),
chains = 4, iter = 4000, warmup = 2000, cores = 4
)
summary(model)Two settings deserve attention:
sampling_variance = "D" – the column
D in data_fhnorm holds the known sampling
variance for each area; hbsaems translates this into an
offset on the observation-level standard deviation so brms / Stan never
tries to estimate it.adapt_delta = 0.99 – Fay-Herriot
likelihoods can produce mild funnel geometry between and the area random
effects when many areas have small . As of v1.1.0 hbsaems
already raises the default from brms’s 0.8 to
0.95; pushing it further to 0.99 buys an even
more conservative leapfrog step and eliminates the occasional divergent
transition without re-parameterising the model.For the remainder of the vignette we use a tiny model (very few iterations) so that the diagnostic and prediction chunks below have a real object to operate on. In your own analysis, use the full production settings shown above, not the toy settings here.
# Mini fit -- iter = 200, chains = 1 -- for vignette demonstration only.
# Do NOT use these settings for inference: the chains have not
# converged at this length and the posterior will be biased.
fit_demo <- suppressWarnings(
hbm(
formula = brms::bf(y ~ x1 + x2 + x3),
data = data_fhnorm,
re = ~ (1 | regency),
sampling_variance = "D",
chains = 1,
iter = 200,
warmup = 100,
refresh = 0,
seed = 1
)
)# Operate on the mini fit_demo above (NOT a substitute for production
# diagnostics on full chains).
diag <- convergence_check(fit_demo)
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> No divergences to plot.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
is_converged(fit_demo)
#> [1] FALSE
summary(diag)
#>
#> ===== Convergence Diagnostics Summary =====
#>
#> Divergent transitions: 0 (0.00%)
#> E-BFMI (min over chains): 0.930
#> Max-treedepth hit rate: 0.00%
#>
#> R-hat:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.9904 0.9961 1.0030 1.0126 1.0175 1.1435
#>
#> Bulk ESS:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 43.85 131.49 171.06 160.25 200.00 200.00
#>
#> Tail ESS:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 21.15 64.12 76.15 75.61 90.20 116.78
#>
#> Geweke test (Z-scores):
#> [[1]]
#>
#> Fraction in 1st window = 0.1
#> Fraction in 2nd window = 0.5
#>
#> b_Intercept b_x1
#> 0.439675 0.675512
#> b_x2 b_x3
#> -0.110359 -1.329367
#> sd_regency__Intercept Intercept
#> -0.351228 0.518596
#> r_regency[regency_001,Intercept] r_regency[regency_002,Intercept]
#> 1.482982 0.578184
#> r_regency[regency_003,Intercept] r_regency[regency_004,Intercept]
#> 0.064829 -1.694880
#> r_regency[regency_005,Intercept] r_regency[regency_006,Intercept]
#> -0.086819 -1.550906
#> r_regency[regency_007,Intercept] r_regency[regency_008,Intercept]
#> -0.225457 -1.634221
#> r_regency[regency_009,Intercept] r_regency[regency_010,Intercept]
#> -0.573915 0.352415
#> r_regency[regency_011,Intercept] r_regency[regency_012,Intercept]
#> 0.329977 -1.152600
#> r_regency[regency_013,Intercept] r_regency[regency_014,Intercept]
#> -0.184377 0.994029
#> r_regency[regency_015,Intercept] r_regency[regency_016,Intercept]
#> -0.734035 -0.459930
#> r_regency[regency_017,Intercept] r_regency[regency_018,Intercept]
#> 1.866025 -0.857966
#> r_regency[regency_019,Intercept] r_regency[regency_020,Intercept]
#> -0.910251 -0.124925
#> r_regency[regency_021,Intercept] r_regency[regency_022,Intercept]
#> -0.638350 -1.228459
#> r_regency[regency_023,Intercept] r_regency[regency_024,Intercept]
#> -2.227878 0.009446
#> r_regency[regency_025,Intercept] r_regency[regency_026,Intercept]
#> -0.570614 0.706367
#> r_regency[regency_027,Intercept] r_regency[regency_028,Intercept]
#> -0.153144 0.666069
#> r_regency[regency_029,Intercept] r_regency[regency_030,Intercept]
#> -0.675936 0.019550
#> r_regency[regency_031,Intercept] r_regency[regency_032,Intercept]
#> 0.157159 0.068315
#> r_regency[regency_033,Intercept] r_regency[regency_034,Intercept]
#> -1.454501 -1.038916
#> r_regency[regency_035,Intercept] r_regency[regency_036,Intercept]
#> -1.488091 -1.626934
#> r_regency[regency_037,Intercept] r_regency[regency_038,Intercept]
#> 0.353194 0.786427
#> r_regency[regency_039,Intercept] r_regency[regency_040,Intercept]
#> -1.751016 0.171068
#> r_regency[regency_041,Intercept] r_regency[regency_042,Intercept]
#> 0.223230 -1.462420
#> r_regency[regency_043,Intercept] r_regency[regency_044,Intercept]
#> -1.376935 -0.370675
#> r_regency[regency_045,Intercept] r_regency[regency_046,Intercept]
#> -1.063761 -0.341718
#> r_regency[regency_047,Intercept] r_regency[regency_048,Intercept]
#> 0.665556 -1.282668
#> r_regency[regency_049,Intercept] r_regency[regency_050,Intercept]
#> -1.338728 -0.110090
#> r_regency[regency_051,Intercept] r_regency[regency_052,Intercept]
#> -0.292501 -0.653364
#> r_regency[regency_053,Intercept] r_regency[regency_054,Intercept]
#> 0.761260 -0.265161
#> r_regency[regency_055,Intercept] r_regency[regency_056,Intercept]
#> -0.781514 -1.031518
#> r_regency[regency_057,Intercept] r_regency[regency_058,Intercept]
#> -1.326952 -0.026405
#> r_regency[regency_059,Intercept] r_regency[regency_060,Intercept]
#> 0.556657 1.998405
#> r_regency[regency_061,Intercept] r_regency[regency_062,Intercept]
#> -0.076548 -0.893495
#> r_regency[regency_063,Intercept] r_regency[regency_064,Intercept]
#> 1.830265 1.063698
#> r_regency[regency_065,Intercept] r_regency[regency_066,Intercept]
#> -1.320349 -2.125862
#> r_regency[regency_067,Intercept] r_regency[regency_068,Intercept]
#> 0.086755 0.826772
#> r_regency[regency_069,Intercept] r_regency[regency_070,Intercept]
#> -2.425150 -1.021848
#> r_regency[regency_071,Intercept] r_regency[regency_072,Intercept]
#> -0.921946 -0.061025
#> r_regency[regency_073,Intercept] r_regency[regency_074,Intercept]
#> 1.882800 -2.702286
#> r_regency[regency_075,Intercept] r_regency[regency_076,Intercept]
#> -1.796570 0.426268
#> r_regency[regency_077,Intercept] r_regency[regency_078,Intercept]
#> 1.079833 0.048222
#> r_regency[regency_079,Intercept] r_regency[regency_080,Intercept]
#> -0.276625 -0.434901
#> r_regency[regency_081,Intercept] r_regency[regency_082,Intercept]
#> 0.373386 -1.262929
#> r_regency[regency_083,Intercept] r_regency[regency_084,Intercept]
#> -1.747763 0.778250
#> r_regency[regency_085,Intercept] r_regency[regency_086,Intercept]
#> 0.027099 0.458924
#> r_regency[regency_087,Intercept] r_regency[regency_088,Intercept]
#> 1.064617 -1.144627
#> r_regency[regency_089,Intercept] r_regency[regency_090,Intercept]
#> -0.268693 0.005188
#> r_regency[regency_091,Intercept] r_regency[regency_092,Intercept]
#> 0.625051 0.007719
#> r_regency[regency_093,Intercept] r_regency[regency_094,Intercept]
#> -1.305631 -0.166948
#> r_regency[regency_095,Intercept] r_regency[regency_096,Intercept]
#> -0.715510 0.749165
#> r_regency[regency_097,Intercept] r_regency[regency_098,Intercept]
#> -1.672694 -2.087388
#> r_regency[regency_099,Intercept] r_regency[regency_100,Intercept]
#> -0.694299 -1.305847
#> lprior lp__
#> -0.052890 -1.287801
#> z_1[1,1] z_1[1,2]
#> 1.025456 0.657983
#> z_1[1,3] z_1[1,4]
#> -0.007537 -1.689116
#> z_1[1,5] z_1[1,6]
#> -0.185234 -1.505061
#> z_1[1,7] z_1[1,8]
#> -0.189963 -1.297634
#> z_1[1,9] z_1[1,10]
#> -0.611099 0.524402
#> z_1[1,11] z_1[1,12]
#> 0.250522 -1.047651
#> z_1[1,13] z_1[1,14]
#> -0.179149 1.020036
#> z_1[1,15] z_1[1,16]
#> -0.693134 -0.585218
#> z_1[1,17] z_1[1,18]
#> 2.100592 -0.860425
#> z_1[1,19] z_1[1,20]
#> -1.204257 -0.143581
#> z_1[1,21] z_1[1,22]
#> -0.639256 -1.453809
#> z_1[1,23] z_1[1,24]
#> -2.248067 -0.017670
#> z_1[1,25] z_1[1,26]
#> -0.369461 0.547181
#> z_1[1,27] z_1[1,28]
#> -0.231902 0.524318
#> z_1[1,29] z_1[1,30]
#> -0.648019 -0.013206
#> z_1[1,31] z_1[1,32]
#> 0.128067 0.077570
#> z_1[1,33] z_1[1,34]
#> -1.271948 -0.932256
#> z_1[1,35] z_1[1,36]
#> -1.468568 -1.599269
#> z_1[1,37] z_1[1,38]
#> 0.301756 0.989599
#> z_1[1,39] z_1[1,40]
#> -1.740444 0.216227
#> z_1[1,41] z_1[1,42]
#> 0.416618 -1.274363
#> z_1[1,43] z_1[1,44]
#> -1.133581 -0.041793
#> z_1[1,45] z_1[1,46]
#> -1.152071 -0.337986
#> z_1[1,47] z_1[1,48]
#> 0.762589 -1.349635
#> z_1[1,49] z_1[1,50]
#> -1.280454 -0.225432
#> z_1[1,51] z_1[1,52]
#> -0.303364 -0.562894
#> z_1[1,53] z_1[1,54]
#> 0.639453 -0.043853
#> z_1[1,55] z_1[1,56]
#> -0.564568 -1.157184
#> z_1[1,57] z_1[1,58]
#> -1.152060 -0.099215
#> z_1[1,59] z_1[1,60]
#> 0.935305 1.811214
#> z_1[1,61] z_1[1,62]
#> -0.092633 -0.761527
#> z_1[1,63] z_1[1,64]
#> 1.720061 1.020080
#> z_1[1,65] z_1[1,66]
#> -1.753982 -2.212438
#> z_1[1,67] z_1[1,68]
#> 0.202992 0.872749
#> z_1[1,69] z_1[1,70]
#> -2.673858 -1.243029
#> z_1[1,71] z_1[1,72]
#> -0.910726 -0.115829
#> z_1[1,73] z_1[1,74]
#> 1.713765 -2.525547
#> z_1[1,75] z_1[1,76]
#> -1.877583 0.737048
#> z_1[1,77] z_1[1,78]
#> 1.104834 -0.037961
#> z_1[1,79] z_1[1,80]
#> -0.308535 -0.377443
#> z_1[1,81] z_1[1,82]
#> 0.281008 -1.239269
#> z_1[1,83] z_1[1,84]
#> -1.444329 0.529094
#> z_1[1,85] z_1[1,86]
#> 0.062075 0.368325
#> z_1[1,87] z_1[1,88]
#> 1.414949 -1.230895
#> z_1[1,89] z_1[1,90]
#> -0.334978 -0.094349
#> z_1[1,91] z_1[1,92]
#> 0.474061 -0.059493
#> z_1[1,93] z_1[1,94]
#> -1.482723 0.052240
#> z_1[1,95] z_1[1,96]
#> -0.933766 0.583424
#> z_1[1,97] z_1[1,98]
#> -1.645358 -2.388447
#> z_1[1,99] z_1[1,100]
#> -0.720579 -1.331146
hbm_warnings(fit_demo)
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> [1] "R-hat > 1.1: model may not have converged."The full set of diagnostics on your production fit would also include trace plots, R-hat / ESS tables, and pp-checks. These all follow the same calling convention:
model2 <- hbm(brms::bf(y ~ x1 + x2), data = data_fhnorm,
re = ~ (1 | regency),
sampling_variance = "D",
control = list(adapt_delta = 0.99),
chains = 4, iter = 4000)
model3 <- hbm(brms::bf(y ~ x1), data = data_fhnorm,
re = ~ (1 | regency),
sampling_variance = "D",
control = list(adapt_delta = 0.99),
chains = 4, iter = 4000)
model_compare(model, model2)
model_compare_all(full = model, medium = model2, simple = model3)The simplest call uses the data the model was fit on – no
newdata argument needed:
# In-sample prediction (default): operates on the data stored
# inside the fitted object.
est <- sae_predict(fit_demo)
summary(est)
#>
#> ===== Small Area Estimation Summary =====
#>
#> Areas : 100
#> Overall RSE : 7.2 %
#>
#> Predictions:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 6.336 8.938 9.769 9.855 10.637 13.537
#>
#> RSE by area:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 4.416 5.981 7.062 7.196 8.171 11.258
head(est$result_table, 5)
#> Prediction SD RSE_percent
#> 1 8.131821 0.6843505 8.415711
#> 2 9.340465 0.7375581 7.896375
#> 3 8.286397 0.7710245 9.304701
#> 4 10.378301 0.7482991 7.210228
#> 5 10.442262 0.7160301 6.857041To predict at fresh data points, pass newdata. When you
do, any internal offset columns the original fit needed
(e.g. .hbsaems_sigma_fixed = sqrt(D) for the Fay-Herriot
sugar) are repopulated automatically when the new data frame has the
same number of rows as the training data – the standard “predict at the
same areas with updated covariates” use case.
For model averaging, fit several candidate models, then either let
model_average() weight them by stacking weights (LOO-based)
or supply your own weights:
# Stacking weights (default, derived from LOO)
avg_stacked <- model_average(model, model2, model3)
# User-specified weights -- must sum to 1
avg_manual <- model_average(model, model2, weights = c(0.7, 0.3))
# The returned object is an `hbsae_results`, identical in shape to
# what sae_predict() returns, so all post-processing helpers
# (sae_transform, sae_scale, sae_filter, sae_aggregate) work on it.
plot(avg_stacked, type = "predictions")All sae_* post-processing helpers operate on the
hbsae_results object returned by either
sae_predict() or model_average():
sessionInfo()
#> R version 4.6.0 (2026-04-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.3 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Asia/Jakarta
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] hbsaems_1.1.0
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.1 dplyr_1.2.1 farver_2.1.2
#> [4] loo_2.9.0 S7_0.2.2 fastmap_1.2.0
#> [7] TH.data_1.1-3 tensorA_0.36.2.1 rpart_4.1.27
#> [10] digest_0.6.39 estimability_1.5.1 lifecycle_1.0.5
#> [13] StanHeaders_2.32.10 survival_3.8-6 processx_3.9.0
#> [16] magrittr_2.0.5 posterior_1.7.0 compiler_4.6.0
#> [19] rlang_1.2.0 sass_0.4.10 tools_4.6.0
#> [22] yaml_2.3.12 knitr_1.51 bridgesampling_1.2-1
#> [25] pkgbuild_1.4.8 curl_7.1.0 plyr_1.8.9
#> [28] RColorBrewer_1.1-3 multcomp_1.4-28 abind_1.4-8
#> [31] withr_3.0.2 purrr_1.2.2 nnet_7.3-20
#> [34] grid_4.6.0 stats4_4.6.0 jomo_2.7-6
#> [37] xtable_1.8-8 mice_3.19.0 inline_0.3.21
#> [40] ggplot2_4.0.3 emmeans_1.11.2 scales_1.4.0
#> [43] iterators_1.0.14 MASS_7.3-65 dichromat_2.0-0.1
#> [46] cli_3.6.6 mvtnorm_1.3-7 rmarkdown_2.31
#> [49] reformulas_0.4.4 generics_0.1.4 otel_0.2.0
#> [52] RcppParallel_5.1.11-2 rstudioapi_0.17.1 reshape2_1.4.5
#> [55] minqa_1.2.8 cachem_1.1.0 rstan_2.32.7
#> [58] stringr_1.6.0 splines_4.6.0 bayesplot_1.15.0
#> [61] parallel_4.6.0 matrixStats_1.5.0 brms_2.23.0
#> [64] vctrs_0.7.3 boot_1.3-32 V8_8.2.0
#> [67] glmnet_5.0 Matrix_1.7-5 sandwich_3.1-1
#> [70] jsonlite_2.0.0 callr_3.7.6 mitml_0.4-5
#> [73] foreach_1.5.2 jquerylib_0.1.4 tidyr_1.3.2
#> [76] glue_1.8.1 pan_1.9 nloptr_2.2.1
#> [79] codetools_0.2-20 ps_1.9.3 distributional_0.7.0
#> [82] stringi_1.8.7 gtable_0.3.6 shape_1.4.6.1
#> [85] QuickJSR_1.9.2 lme4_2.0-1 tibble_3.3.1
#> [88] pillar_1.11.1 htmltools_0.5.9 Brobdingnag_1.2-9
#> [91] R6_2.6.1 Rdpack_2.6.6 evaluate_1.0.5
#> [94] lattice_0.22-9 rbibutils_2.4.1 backports_1.5.1
#> [97] broom_1.0.13 bslib_0.10.0 rstantools_2.6.0
#> [100] Rcpp_1.1.1-1.1 coda_0.19-4.1 gridExtra_2.3
#> [103] nlme_3.1-169 checkmate_2.3.4 xfun_0.57
#> [106] zoo_1.8-14 pkgconfig_2.0.3