A typical Bayesian SAE workflow with hbsaems v1.0.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 . Raising adapt_delta
from the brms default 0.95 to 0.99 buys a 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.
#> 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.
is_converged(fit_demo)
#> [1] FALSE
summary(diag)
#>
#> ===== Convergence Diagnostics Summary =====
#>
#> R-hat:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.9902 0.9955 1.0038 1.0121 1.0194 1.1280
#>
#> Bulk ESS:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 30.97 100.65 133.85 138.09 188.47 200.00
#>
#> Geweke test (Z-scores):
#> [[1]]
#>
#> Fraction in 1st window = 0.1
#> Fraction in 2nd window = 0.5
#>
#> b_Intercept b_x1
#> -0.37988 -0.88792
#> b_x2 b_x3
#> -2.25456 1.19864
#> sd_regency__Intercept Intercept
#> -0.48291 -0.27155
#> r_regency[regency_001,Intercept] r_regency[regency_002,Intercept]
#> 2.53333 -1.09432
#> r_regency[regency_003,Intercept] r_regency[regency_004,Intercept]
#> -0.34630 -1.03951
#> r_regency[regency_005,Intercept] r_regency[regency_006,Intercept]
#> 0.29714 0.79992
#> r_regency[regency_007,Intercept] r_regency[regency_008,Intercept]
#> 0.96555 0.69001
#> r_regency[regency_009,Intercept] r_regency[regency_010,Intercept]
#> 2.12666 -0.71444
#> r_regency[regency_011,Intercept] r_regency[regency_012,Intercept]
#> 0.78253 0.97046
#> r_regency[regency_013,Intercept] r_regency[regency_014,Intercept]
#> 1.26104 1.69245
#> r_regency[regency_015,Intercept] r_regency[regency_016,Intercept]
#> -0.38183 0.64330
#> r_regency[regency_017,Intercept] r_regency[regency_018,Intercept]
#> 0.69435 -0.23802
#> r_regency[regency_019,Intercept] r_regency[regency_020,Intercept]
#> 1.50055 -0.33108
#> r_regency[regency_021,Intercept] r_regency[regency_022,Intercept]
#> 0.76543 -0.31163
#> r_regency[regency_023,Intercept] r_regency[regency_024,Intercept]
#> 0.95497 0.55438
#> r_regency[regency_025,Intercept] r_regency[regency_026,Intercept]
#> -1.30189 -0.23310
#> r_regency[regency_027,Intercept] r_regency[regency_028,Intercept]
#> -0.27869 -0.37691
#> r_regency[regency_029,Intercept] r_regency[regency_030,Intercept]
#> -1.03580 -0.62288
#> r_regency[regency_031,Intercept] r_regency[regency_032,Intercept]
#> 1.13821 0.68404
#> r_regency[regency_033,Intercept] r_regency[regency_034,Intercept]
#> 1.74722 0.99892
#> r_regency[regency_035,Intercept] r_regency[regency_036,Intercept]
#> -0.18233 -0.56501
#> r_regency[regency_037,Intercept] r_regency[regency_038,Intercept]
#> 0.79287 -0.36884
#> r_regency[regency_039,Intercept] r_regency[regency_040,Intercept]
#> 0.53184 0.22485
#> r_regency[regency_041,Intercept] r_regency[regency_042,Intercept]
#> -0.51635 0.07436
#> r_regency[regency_043,Intercept] r_regency[regency_044,Intercept]
#> -0.24141 -0.03389
#> r_regency[regency_045,Intercept] r_regency[regency_046,Intercept]
#> 0.74541 -0.12837
#> r_regency[regency_047,Intercept] r_regency[regency_048,Intercept]
#> -0.66193 -0.34274
#> r_regency[regency_049,Intercept] r_regency[regency_050,Intercept]
#> 0.96997 -0.72179
#> r_regency[regency_051,Intercept] r_regency[regency_052,Intercept]
#> -0.61008 1.55256
#> r_regency[regency_053,Intercept] r_regency[regency_054,Intercept]
#> 3.14094 1.29244
#> r_regency[regency_055,Intercept] r_regency[regency_056,Intercept]
#> 0.97997 1.76742
#> r_regency[regency_057,Intercept] r_regency[regency_058,Intercept]
#> 0.23508 0.15553
#> r_regency[regency_059,Intercept] r_regency[regency_060,Intercept]
#> -1.38684 -0.03441
#> r_regency[regency_061,Intercept] r_regency[regency_062,Intercept]
#> -0.39761 1.94751
#> r_regency[regency_063,Intercept] r_regency[regency_064,Intercept]
#> 0.88468 -0.74694
#> r_regency[regency_065,Intercept] r_regency[regency_066,Intercept]
#> -1.11521 0.05771
#> r_regency[regency_067,Intercept] r_regency[regency_068,Intercept]
#> 0.54408 0.15825
#> r_regency[regency_069,Intercept] r_regency[regency_070,Intercept]
#> 0.38837 1.89231
#> r_regency[regency_071,Intercept] r_regency[regency_072,Intercept]
#> 0.86007 -1.10758
#> r_regency[regency_073,Intercept] r_regency[regency_074,Intercept]
#> 0.52490 1.19575
#> r_regency[regency_075,Intercept] r_regency[regency_076,Intercept]
#> 0.31622 0.71654
#> r_regency[regency_077,Intercept] r_regency[regency_078,Intercept]
#> -0.02686 1.11578
#> r_regency[regency_079,Intercept] r_regency[regency_080,Intercept]
#> -0.93260 -1.60636
#> r_regency[regency_081,Intercept] r_regency[regency_082,Intercept]
#> -0.73044 -1.85083
#> r_regency[regency_083,Intercept] r_regency[regency_084,Intercept]
#> 0.38941 -0.18867
#> r_regency[regency_085,Intercept] r_regency[regency_086,Intercept]
#> -0.19211 -6.53790
#> r_regency[regency_087,Intercept] r_regency[regency_088,Intercept]
#> 0.06123 -0.77727
#> r_regency[regency_089,Intercept] r_regency[regency_090,Intercept]
#> 1.48321 -0.81653
#> r_regency[regency_091,Intercept] r_regency[regency_092,Intercept]
#> -0.69895 -1.54475
#> r_regency[regency_093,Intercept] r_regency[regency_094,Intercept]
#> -0.73712 0.97211
#> r_regency[regency_095,Intercept] r_regency[regency_096,Intercept]
#> -0.52273 2.09873
#> r_regency[regency_097,Intercept] r_regency[regency_098,Intercept]
#> 0.92111 -3.44677
#> r_regency[regency_099,Intercept] r_regency[regency_100,Intercept]
#> 0.33993 -0.13166
#> lprior lp__
#> 0.67585 -0.36662
#> z_1[1,1] z_1[1,2]
#> 2.60661 -1.15950
#> z_1[1,3] z_1[1,4]
#> -0.44316 -0.83686
#> z_1[1,5] z_1[1,6]
#> 0.32197 0.72187
#> z_1[1,7] z_1[1,8]
#> 0.85500 0.68356
#> z_1[1,9] z_1[1,10]
#> 1.99573 -0.61475
#> z_1[1,11] z_1[1,12]
#> 0.71367 1.11718
#> z_1[1,13] z_1[1,14]
#> 1.45869 1.84480
#> z_1[1,15] z_1[1,16]
#> -0.17716 1.26162
#> z_1[1,17] z_1[1,18]
#> 0.78684 -0.50915
#> z_1[1,19] z_1[1,20]
#> 0.84447 -0.03705
#> z_1[1,21] z_1[1,22]
#> 0.78376 -0.34361
#> z_1[1,23] z_1[1,24]
#> 1.16054 0.83935
#> z_1[1,25] z_1[1,26]
#> -1.82156 0.07895
#> z_1[1,27] z_1[1,28]
#> -0.37286 0.03831
#> z_1[1,29] z_1[1,30]
#> -0.97086 -0.47947
#> z_1[1,31] z_1[1,32]
#> 1.66448 0.47742
#> z_1[1,33] z_1[1,34]
#> 2.42846 1.23626
#> z_1[1,35] z_1[1,36]
#> -0.17968 -0.45039
#> z_1[1,37] z_1[1,38]
#> 0.60664 -0.85224
#> z_1[1,39] z_1[1,40]
#> 0.33977 0.23800
#> z_1[1,41] z_1[1,42]
#> -0.36941 -0.04355
#> z_1[1,43] z_1[1,44]
#> -0.35746 0.13954
#> z_1[1,45] z_1[1,46]
#> 0.73972 0.31833
#> z_1[1,47] z_1[1,48]
#> -1.22041 -0.37691
#> z_1[1,49] z_1[1,50]
#> 0.94695 -0.28747
#> z_1[1,51] z_1[1,52]
#> -0.74239 1.63518
#> z_1[1,53] z_1[1,54]
#> 2.46719 1.52550
#> z_1[1,55] z_1[1,56]
#> 1.01434 1.35748
#> z_1[1,57] z_1[1,58]
#> 0.02757 0.20723
#> z_1[1,59] z_1[1,60]
#> -1.18128 -0.19537
#> z_1[1,61] z_1[1,62]
#> -0.52778 2.07618
#> z_1[1,63] z_1[1,64]
#> 1.17741 -0.82973
#> z_1[1,65] z_1[1,66]
#> -1.40315 -0.24890
#> z_1[1,67] z_1[1,68]
#> 0.68357 0.17158
#> z_1[1,69] z_1[1,70]
#> 0.25941 1.32047
#> z_1[1,71] z_1[1,72]
#> 0.94439 -1.28296
#> z_1[1,73] z_1[1,74]
#> 0.56654 1.18829
#> z_1[1,75] z_1[1,76]
#> 0.20314 0.82644
#> z_1[1,77] z_1[1,78]
#> 0.36700 1.16686
#> z_1[1,79] z_1[1,80]
#> -0.37497 -1.69950
#> z_1[1,81] z_1[1,82]
#> -0.68790 -1.75051
#> z_1[1,83] z_1[1,84]
#> -0.03720 -0.07350
#> z_1[1,85] z_1[1,86]
#> -0.44064 -5.94021
#> z_1[1,87] z_1[1,88]
#> 0.21214 -0.72701
#> z_1[1,89] z_1[1,90]
#> 1.59174 -0.84100
#> z_1[1,91] z_1[1,92]
#> -0.46963 -1.53098
#> z_1[1,93] z_1[1,94]
#> -0.90421 1.50403
#> z_1[1,95] z_1[1,96]
#> -0.42522 1.75905
#> z_1[1,97] z_1[1,98]
#> 1.37087 -3.50478
#> z_1[1,99] z_1[1,100]
#> 0.54213 -0.44517
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.
#> [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 : 12.2 %
#>
#> Predictions:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 6.329 8.936 9.808 9.849 10.691 13.615
#>
#> RSE by area:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 6.729 9.762 11.969 12.204 14.390 21.050
head(est$result_table, 5)
#> Prediction SD RSE_percent
#> 1 8.061463 1.1702579 14.516694
#> 2 9.326271 1.7778278 19.062580
#> 3 8.216403 1.7295135 21.049521
#> 4 10.437027 0.9165638 8.781848
#> 5 10.694208 1.2784067 11.954197To 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.0.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