Chapter 07 relaxed the assumption that the between-group variances \(\tau^2_k\) are known. This chapter does the same for the within-group variance \(\sigma^2\) — the residual scatter of individual observations around their group’s fitted line.
Three levels of ambition:
| Specification | \(\sigma^2\) | Argument |
|---|---|---|
| Fixed scalar | One number, held at the lmer estimate |
dispersion_ranef = ps$group.dispersion |
| Sampled scalar | One number with a posterior | dispersion_ranef = dGamma(...) |
| Per-group | One \(\sigma^2_j\) per group, each sampled | dispformula = ~group, dGamma_list() |
The per-group case is the interesting one, and it is what this
chapter demonstrates. Because it runs an accept-reject dispersion draw
for every group on every sweep, the fit takes about a minute, so it is
precomputed by data-raw/make_vignette_outputs.R and loaded
here.
data(sleepstudy, package = "lme4")
dat <- sleepstudy
dat$Days_c <- dat$Days - mean(dat$Days)
form <- Reaction ~ Days_c + (1 + Days_c || Subject)
ps <- Prior_Setup_GLMM(form, data = dat, pop.pwt = 0.01)
ps$group.dispersion
#> [1] 709.826
#> attr(,"source")
#> [1] "pooled S_marg (A12 3.3.4)"
#> attr(,"classical")
#> [1] 654.941One number, 709.8, applied to all 180 observations. It assumes every subject’s reaction times scatter by the same amount around their own line — which, for a sleep-deprivation study, is implausible on its face. Some people get erratic when tired and some do not, and that heterogeneity is a finding, not a nuisance.
Ask Prior_Setup_GLMM() for a per-group calibration with
dispformula:
ps_group <- Prior_Setup_GLMM(form, data = dat, pop.pwt = 0.01,
dispformula = ~Subject)
dl <- dGamma_list(ps_group)
length(dl)
dl[[1]]#> 18
#> Call: glmbayesCore::dGamma(shape = g$shape_ING, rate = g$rate, beta = matrix(0,
#> 1, 1, dimnames = list("(Intercept)", NULL)), Inv_Dispersion = TRUE,
#> max_disp_perc = mdp_lev, disp_lower = disp_lower, disp_upper = disp_upper)
#> Prior Family: dGamma
#> Prior List:
#> $shape
#> [1] 1.550505
#> $rate
#> [1] 1501.98
#> $beta
#> [,1]
#> (Intercept) 0
#> $Inv_Dispersion
#> [1] TRUE
#> $max_disp_perc
#> [1] 0.99
#> $disp_lower
#> [1] 1088.082
#> $disp_upper
#> [1] 7277.047
#> attr(,"Prior Type")
#> [1] "dGamma"dGamma_list() returns one dGamma prior per
group. Each carries a shape and rate for a
Gamma prior on the precision \(1/\sigma^2_j\) (hence
Inv_Dispersion = TRUE), plus a truncation window
disp_lower / disp_upper calibrated from that
group’s own data.
The windows differ across groups by more than an order of magnitude,
because they are quantiles of each group’s limiting posterior and the
subjects genuinely differ. Calling dGamma_list() on a
dispformula = ~1 setup is an error rather than a silent
fallback:
set.seed(2024)
fit <- lmerb(
form, data = dat,
pfamily_list = pfamily_list(ps_group),
dispersion_ranef = dl,
dispformula = ~Subject,
n = 200L,
progbar = FALSE
)
fit#> Call:
#> lmerb(formula = form, data = dat, pfamily_list = pfamily_list(ps_group),
#> dispersion_ranef = dl, dispformula = ~Subject, n = 200L,
#> progbar = FALSE)
#> Bayesian linear mixed model [200 draws, two-block Gibbs]
#> Formula: Reaction ~ Days_c + (1 + Days_c || Subject)
#> Random effects (variance components fixed at lmer estimates):
#> Conditional model:
#> Groups Name Std.Dev.
#> Subject (Intercept) 38.265
#> Days_c 6.059
#> Residual NA
#> Number of obs: 180, groups: Subject, 18
#> --- Posterior means (ICM exact, under fixed variance components) ---
#> RE component parameter fixef.mode fixef.means draws SD
#> ------------ ----------- ------------ ------------ ----------
#> (Intercept) (Intercept) 298.4862 297.4782 9.0110
#> Days_c (Intercept) 10.3898 10.0397 1.4591The engine is two-block Gibbs with
m_convergence = 4. Sampling \(\sigma^2\) breaks the joint normality that
made the earlier fits exact, so sweeps are required — the same mechanism
as Chapter 07, applied at the observation level instead of the group
level.
Note the Residual line now reads NA: there
is no single residual standard deviation to report, because there are 18
of them.
summary_sigma2() exists for exactly this case. On a
scalar-dispersion fit it stops with a message; here it produces a full
report.
#> Measurement dispersion (sigma^2, per group)
#> Formula: Reaction ~ Days_c + (1 + Days_c || Subject)
#> Groups: Subject (18 level(s)) [200 draws]
#> Prior and lmer reference:
#> Prior n_prior 1/E[1/sigma2] E[sigma2] sqrt(E[sigma2]) disp_lower
#> 308 dGamma 0.101 968.7038 2728.3675 52.2338 1088.0817
#> 309 dGamma 0.101 28.0287 78.9431 8.8850 31.4828
#> 310 dGamma 0.101 54.8003 154.3459 12.4236 61.5536
#> 330 dGamma 0.101 177.2113 499.1180 22.3409 199.0498
#> 331 dGamma 0.101 197.2977 555.6917 23.5731 221.6116
#> 332 dGamma 0.101 1363.6127 3840.6339 61.9728 1531.6571
#> 333 dGamma 0.101 57.1848 161.0618 12.6910 64.2320
#> 334 dGamma 0.101 160.0552 450.7975 21.2320 179.7795
#> 335 dGamma 0.101 46.2221 130.1852 11.4099 51.9183
#> 337 dGamma 0.101 109.9134 309.5725 17.5947 123.4585
#> 349 dGamma 0.101 75.4226 212.4288 14.5749 84.7173
#> 350 dGamma 0.101 246.1052 693.1587 26.3279 276.4339
#> 351 dGamma 0.101 187.2562 527.4097 22.9654 210.3327
#> 352 dGamma 0.101 249.9618 704.0209 26.5334 280.7658
#> 369 dGamma 0.101 94.1328 265.1264 16.2827 105.7333
#> 370 dGamma 0.101 236.2855 665.5014 25.7973 265.4041
#> 371 dGamma 0.101 231.5002 652.0234 25.5347 260.0290
#> 372 dGamma 0.101 47.8088 134.6543 11.6041 53.7006
#> disp_upper n_data lmer lmer SD
#> 308 7277.0467 10 654.9399 25.5918
#> 309 210.5554 10 654.9399 25.5918
#> 310 411.6681 10 654.9399 25.5918
#> 330 1331.2373 10 654.9399 25.5918
#> 331 1482.1296 10 654.9399 25.5918
#> 332 10243.6613 10 654.9399 25.5918
#> 333 429.5809 10 654.9399 25.5918
#> 334 1202.3581 10 654.9399 25.5918
#> 335 347.2273 10 654.9399 25.5918
#> 337 825.6855 10 654.9399 25.5918
#> 349 566.5858 10 654.9399 25.5918
#> 350 1848.7789 10 654.9399 25.5918
#> 351 1406.6965 10 654.9399 25.5918
#> 352 1877.7504 10 654.9399 25.5918
#> 369 707.1399 10 654.9399 25.5918
#> 370 1775.0120 10 654.9399 25.5918
#> 371 1739.0636 10 654.9399 25.5918
#> 372 359.1472 10 654.9399 25.5918
#> Overview:
#> n_total Post.Mode Post.Mean Post.Sd MC Error Mean SD Pr(Prior_tail)
#> 308 10.101 NA 2416.164 1136.283 80.347 48.040 0.300
#> 309 10.101 NA 76.385 33.520 2.370 8.560 0.350
#> 310 10.101 NA 150.516 67.874 4.799 11.994 0.360
#> 330 10.101 NA 509.060 225.143 15.920 22.070 0.410
#> 331 10.101 NA 551.793 253.699 17.939 22.953 0.360
#> 332 10.101 NA 3484.084 1389.546 98.256 58.000 0.305
#> 333 10.101 NA 162.006 61.802 4.370 12.513 0.420
#> 334 10.101 NA 425.047 197.092 13.937 20.154 0.310
#> 335 10.101 NA 119.281 46.612 3.296 10.739 0.345
#> 337 10.101 NA 266.730 108.863 7.698 16.040 0.250
#> 349 10.101 NA 206.168 91.679 6.483 14.049 0.370
#> 350 10.101 NA 595.623 253.575 17.930 23.936 0.250
#> 351 10.101 NA 481.447 192.553 13.616 21.543 0.300
#> 352 10.101 NA 620.597 251.072 17.753 24.458 0.315
#> 369 10.101 NA 252.266 112.848 7.980 15.537 0.345
#> 370 10.101 NA 600.942 265.415 18.768 24.016 0.295
#> 371 10.101 NA 606.422 291.508 20.613 24.042 0.305
#> 372 10.101 NA 126.839 58.315 4.123 11.009 0.325
#> Cand/draw
#> 308 454.199
#> 309 142.227
#> 310 102.754
#> 330 127.725
#> 331 68.914
#> 332 8.523
#> 333 55.625
#> 334 49.414
#> 335 364.715
#> 337 1443.852
#> 349 79.311
#> 350 626.420
#> 351 63.456
#> 352 102.732
#> 369 51.851
#> 370 326.788
#> 371 32.358
#> 372 67.639
#> Distribution percentiles (sigma^2):
#> 1.0% 2.5% 5.0% Median 95.0% 97.5% 99.0%
#> 308 1132.2951 1163.9960 1252.8164 2043.3335 4611.9647 5586.9382 6360.7105
#> 309 35.3387 39.1729 41.2409 65.6035 146.5303 170.0677 179.0980
#> 310 62.9709 65.7958 68.6617 134.9033 302.8436 344.7183 359.4439
#> 330 224.6992 237.3274 260.3205 455.2978 985.9787 1053.1138 1194.4626
#> 331 231.7889 257.8157 278.4663 472.8253 1118.5123 1246.0137 1331.9463
#> 332 1583.8944 1793.5572 1843.9301 3109.2937 5989.4070 6787.0218 7938.0145
#> 333 70.0109 78.6750 87.0476 148.4158 272.9803 307.1649 338.9115
#> 334 193.7870 199.2195 212.9288 360.8734 828.4528 990.4675 1090.7642
#> 335 60.1902 63.9069 66.4300 103.7073 214.8687 246.7130 264.6717
#> 337 127.2926 133.8304 138.0755 240.6378 448.2460 587.8152 624.8474
#> 349 91.7745 97.5812 101.8211 185.3812 381.0241 419.3905 516.8079
#> 350 296.6683 304.3242 316.1304 529.5270 1184.2007 1294.6091 1398.3420
#> 351 225.1037 237.8322 252.2384 445.7604 872.0049 947.5692 1083.2297
#> 352 292.7784 312.4403 336.9877 579.1063 1107.0397 1176.3035 1363.3322
#> 369 109.2150 110.4281 122.5052 222.1470 477.6992 548.7330 621.7914
#> 370 284.6313 293.9897 303.5037 528.8531 1047.1712 1370.6783 1590.2976
#> 371 275.7320 288.7767 313.9880 529.1661 1252.6591 1497.1955 1572.1433
#> 372 57.2172 59.8910 61.4589 113.9158 246.4217 292.9594 321.2456Three blocks, one per question.
Prior and reference lists, per subject, the Gamma
prior in use, the prior mean, the truncation window, the number of
observations, and the pooled lmer residual variance for
comparison. n_data = 10 everywhere — ten observations is
very little information about a variance, which is why the prior weight
matters here more than anywhere else in the package.
Overview gives the posterior mean and standard
deviation of each \(\sigma^2_j\), Monte
Carlo error, and Cand/draw. That last column is
accept-reject efficiency: values in the tens or hundreds mean the
envelope is loose for that group, usually because its posterior sits
near a truncation boundary.
Percentiles shows the spread directly, and the spread is the finding. The posterior means run from roughly 80 to 3600 — subject 308’s reaction times scatter more than an order of magnitude more than subject 309’s. A single pooled \(\sigma^2\) of 710 splits the difference and describes neither.
That heterogeneity also propagates: groups with large \(\sigma^2_j\) carry less weight in the population estimate and get shrunk harder toward it, which is the statistically correct behaviour and is not something a pooled model can express.
Between the two extremes there is a middle option: one \(\sigma^2\) shared by all groups, but
sampled rather than fixed. Pass a single dGamma() object
instead of a list:
fit_scalar_sampled <- lmerb(
form, data = dat,
pfamily_list = pfamily_list(ps),
dispersion_ranef = dGamma(shape = 2, rate = 2 * ps$group.dispersion,
Inv_Dispersion = TRUE),
n = 200L, progbar = FALSE
)This is the usual Bayesian treatment of a residual variance, and it is what you want when you believe the groups share a scatter but do not want to pretend you know its value. It costs sweeps, like any non-Gaussian target, but far less than the per-group version.
| Situation | Use |
|---|---|
| Many observations per group, homogeneous scatter | Fixed scalar |
| Few groups, homogeneous scatter, honest intervals wanted | Sampled scalar |
| Groups visibly differ in residual spread | dispformula = ~group with
dGamma_list() |
Diagnosing which you are in is easy before fitting: compute the
residual variance per group from an lmer fit and look at
the ratio of largest to smallest. Anything beyond a factor of three or
four is worth modeling.
m <- lme4::lmer(form, data = dat, REML = TRUE)
v <- tapply(residuals(m), dat$Subject, var)
round(c(min = min(v), median = median(v), max = max(v),
ratio = max(v) / min(v)), 1)
#> min median max ratio
#> 91.0 417.9 3296.6 36.2A ratio of that size is not something a single pooled variance can absorb.
vignette("Chapter-B07", package = "lmebayesCore") — the
per-group dispersion algebra and the marginal precision blocks.?dGamma_list, ?summary_sigma2