Every model so far has assumed a single measurement variance \(\sigma^2\) shared by all groups. This appendix relaxes that to one \(\sigma^2_j\) per group and works through what has to change: the reference fit, the prior calibration, the truncation window, and the Block 1 draw.
The user-facing view is
vignette("Chapter-08", package = "lmebayes"). This chapter
is about how the numbers are produced.
dispformula selectsPooled (dispformula = ~1, the default):
\[ y_j \mid \beta_j \sim N(D_j\beta_j,\ \sigma^2 I_{n_j}), \qquad \beta_j \mid \gamma, \Psi \sim N(\mathcal{W}_j\gamma, \Psi). \]
Per group (dispformula = ~<group_name>):
\[ y_j \mid \beta_j, \Omega_j \sim N(D_j\beta_j,\ \Omega_j^{-1} I_{n_j}), \qquad \Omega_j = 1/\sigma^2_j \sim \mathrm{Gamma}(a_j^0, r_j^0) \ \text{(truncated)}, \]
with the second level untouched. The likelihood already factorized by group; now the prior does too, and each group carries its own ING prior on \((\beta_j, \sigma^2_j)\).
dispformula accepts only ~1 or
~<group_name> naming the random-effects grouping
factor, and only for Gaussian responses — a Poisson or binomial GLMM has
no free observation-level dispersion to calibrate, so
.lmebayes_dispformula_kind() rejects it.
Switching it on changes five things at once:
dispformula = ~1 |
dispformula = ~group |
|
|---|---|---|
| Reference fit | lme4::lmer() |
glmmTMB::glmmTMB() with the same
dispformula |
calibration_source |
"lme4" |
"glmmTMB" |
group.ing_prior |
one flat list | named list, one entry per level |
group.dispersion.pwt |
scalar | scalar or length-\(J\) |
group.max_disp_perc |
scalar | scalar or length-\(J\) |
data(sleepstudy, package = "lme4")
dat <- sleepstudy
dat$Days_c <- dat$Days - mean(dat$Days)
form <- Reaction ~ Days_c + (1 + Days_c || Subject)
ps_pooled <- Prior_Setup_GLMM(form, data = dat, pop.pwt = 0.01)
ps_group <- Prior_Setup_GLMM(form, data = dat, pop.pwt = 0.01,
dispformula = ~Subject)
c(pooled = ps_pooled$calibration_source,
group = ps_group$calibration_source)
#> pooled group
#> "lme4" "glmmTMB"The pooled setup carries one number; the grouped setup carries one list per subject:
The obvious estimator for group \(j\) is its own residual variance,
\[ \hat\sigma^2_{\mathrm{ols},j} = \frac{\mathrm{RSS}_{\mathrm{ols},j}}{n_j - p_{re}}, \qquad \mathrm{RSS}_{\mathrm{ols},j} = \sum_i (y_{ij} - x_{ij}'\hat\beta_j^{\mathrm{ols}})^2, \]
and it is the wrong target here, for a reason specific to the hierarchical setting. \(\mathrm{RSS}_{\mathrm{ols},j}\) is computed at group \(j\)’s own least-squares fit. The sampler never uses that \(\beta_j\): it uses a shrunk one, pulled toward \(\mathcal{W}_j\gamma\). Residuals around a shrunk coefficient are larger than residuals around the OLS coefficient, so calibrating from the OLS RSS systematically understates the dispersion the sampler will actually encounter.
The fix is to integrate \(\beta_j\) out before measuring the residual. With prior mean \(\mu_j\), dispersion-free prior covariance \(\Sigma_{0j}\), and
\[ M_j = \bigl(\Sigma_{0j} + (X_j'X_j)^{-1}\bigr)^{-1}, \]
the marginal sum of squares is
\[ S_{\mathrm{marg},j} = \mathrm{RSS}_{\mathrm{ols},j} + (\hat\beta_j - \mu_j)'\,M_j\,(\hat\beta_j - \mu_j), \qquad \hat\sigma^2_j = \frac{S_{\mathrm{marg},j}}{n_j - p_{re}}. \]
The added term is a quadratic form in a positive definite matrix, hence always non-negative, so \(\hat\sigma^2_j \ge \hat\sigma^2_{\mathrm{ols},j}\) with equality only when \(\hat\beta_j = \mu_j\). The correction is exactly the shrinkage penalty: groups whose own fit disagrees most with the population get the largest upward adjustment, which is right, because those are the groups whose shrunk residuals will be largest.
With \(n_{\mathrm{prior},j} =
\dfrac{w_j}{1-w_j}\,n_j\) and \(w_j
=\) group.dispersion.pwt (default
0.01):
\[ \mathrm{shape}_j = \frac{n_{\mathrm{prior},j} + 1}{2} + \frac{p_{re}}{2}, \qquad \mathrm{rate}_j = \frac{S_{\mathrm{marg},j}}{2} \cdot \frac{n_{\mathrm{prior},j} + p_{re} - 1}{n_j - p_{re}}. \]
These satisfy the identity that makes the calibration interpretable:
\[ \frac{\mathrm{rate}_j}{\mathrm{shape}_j - 1} = \hat\sigma^2_j \quad\text{exactly.} \]
The Gamma prior mean is the marginal dispersion estimate. Check it:
rate versus rate_gammaEach entry stores two rates, and using the wrong one is a real bug rather than a stylistic choice.
rate_gamma is the fixed-\(\beta\) alternative: the RSS evaluated at
the Zellner blend \(\beta^\star_j =
(1-w_j)\hat\beta_j + w_j\mu_j\) rather than the marginal \(S_{\mathrm{marg},j}\). It is a coherent
quantity, but it answers a different question.
Which one is correct depends on what the sampler does, and here the
sampler is decisive. Block 1 draws \((\beta_j,
\sigma^2_j)\) jointly through
rindepNormalGamma_reg(), which integrates \(\beta_j\) out before drawing \(\sigma^2_j\). A prior calibrated at a fixed
\(\beta\) is then mismatched to the
draw it parameterizes. So rate is what reaches the
sampler; rate_gamma is retained for diagnostics
only.
The gap is not academic — on groups whose \(\hat\beta_j\) departs sharply from \(\mu_j\), differences of 40–55% in the rate have been observed.
rates <- t(vapply(ps_group$group.ing_prior, function(g) {
c(rate = g$rate, rate_gamma = g$rate_gamma,
pct_diff = 100 * (g$rate - g$rate_gamma) / g$rate_gamma)
}, numeric(3)))
round(head(rates[order(-abs(rates[, "pct_diff"])), ], 5), 2)
#> rate rate_gamma pct_diff
#> 308 1501.98 1257.02 19.49
#> 350 381.59 326.91 16.73
#> 337 170.42 146.85 16.06
#> 370 366.36 320.41 14.34
#> 349 116.94 108.20 8.08There is a second correction, applied unconditionally in current releases. \(\mu_j\) is a stand-in for \(\mathcal{W}_j\gamma\), and \(\gamma\) is itself uncertain. Propagating that uncertainty inflates the prior covariance:
\[ \Omega_j[k,k] = \mathcal{W}_j[k]\,\Sigma_{\mathrm{fixef},k}\,\mathcal{W}_j[k]', \qquad M_j' = \bigl(\Sigma_{0j} + \Omega_{0,j} + (X_j'X_j)^{-1}\bigr)^{-1}, \]
with \(\Omega_{0,j} = \Omega_j / \mathrm{dispersion}_j^{\mathrm{classical}}\), and \(S_{\mathrm{marg},j}'\) recomputed with \(M_j'\). The shape is unchanged; only the rate moves. A larger \(\Omega_j\) means a larger \(M_j^{-1}\), hence a smaller \(M_j'\), hence a smaller penalty — being less certain about \(\gamma\) makes a given deviation \(\hat\beta_j - \mu_j\) less surprising, so less of it is charged to \(\sigma^2_j\).
The \(\Omega_j\) matrix is stored per group for inspection:
Chapter B05 introduced
disp_lower/disp_upper as the interval that
keeps the ING chain geometrically ergodic. For measurement dispersion
there is a wrinkle worth knowing, because it is a genuine asymmetry in
the package.
The window is not a quantile interval of the prior \(\mathrm{Gamma}(\mathrm{shape}_j, \mathrm{rate}_j)\). It uses a posterior-shape Gamma, mean-matched at the same center:
\[ \mathrm{shape}^{\mathrm{post}}_j = \mathrm{shape}_j + \frac{n_j}{2}, \qquad \mathrm{rate}^{\mathrm{post}}_j = \hat\sigma^2_j\,(\mathrm{shape}^{\mathrm{post}}_j - 1), \]
and then, on the precision scale with \(Q\) the qgamma quantile
function,
\[ \mathrm{disp\_lower}_j = 1 / Q(\mathtt{max\_disp\_perc}), \qquad \mathrm{disp\_upper}_j = 1 / Q(1 - \mathtt{max\_disp\_perc}). \]
The + n_j/2 matches the envelope builder’s
shape2 = Shape + n_w/2 in
EnvelopeDispersionBuild.cpp, so the window brackets where
the posterior will live rather than where the diffuse prior
would. With \(w_j = 0.01\), the prior
alone carries almost no information and its quantiles would be uselessly
wide.
To be explicit, since the two are easy to conflate:
shape and rate handed to dGamma()
are the prior parameters; only disp_lower
and disp_upper are computed from
shape_post/rate_post.
w <- t(vapply(ps_group$group.ing_prior, function(g) {
c(n_j = g$n_j, sigma2_hat = g$sigma2_hat,
lower = g$disp_lower, upper = g$disp_upper,
shape_post = g$shape + g$n_j / 2)
}, numeric(5)))
round(head(w, 5), 2)
#> n_j sigma2_hat lower upper shape_post
#> 308 10 2728.37 1088.08 7277.05 5.55
#> 309 10 78.94 31.48 210.56 5.55
#> 310 10 154.35 61.55 411.67 5.55
#> 330 10 499.12 199.05 1331.24 5.55
#> 331 10 555.69 221.61 1482.13 5.55Two consequences of the + n_j/2 inflation:
dispformula = ~1 the window is still built from prior
quantiles with no + n/2 term. This is a documented gap, not
a design principle.dGamma_list()dGamma_list() turns a calibrated per-group setup into
one dGamma() pfamily per group — the object Block 1
consumes.
dl <- dGamma_list(ps_group)
class(dl)
#> [1] "dGamma_list" "list"
length(dl)
#> [1] 18
dl[[1L]]
#>
#> Call: 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"Its one interesting argument is
max_disp_perc_measurement. Left NULL, stored
bounds are reused verbatim. Given a scalar or a length-\(J\) vector in \((0.5, 1)\), the window is recomputed for
every group whose stored max_disp_perc differs:
dl_tight <- dGamma_list(ps_group, max_disp_perc_measurement = 0.80)
wd <- attr(dl, "window_diagnostics")
wdt <- attr(dl_tight, "window_diagnostics")
head(data.frame(
group = wd$group,
lower_99 = round(wd$disp_lower), upper_99 = round(wd$disp_upper),
lower_80 = round(wdt$disp_lower), upper_80 = round(wdt$disp_upper)
), 5)
#> group lower_99 upper_99 lower_80 upper_80
#> 1 308 1088 7277 1771 3474
#> 2 309 31 211 51 101
#> 3 310 62 412 100 197
#> 4 330 199 1331 324 636
#> 5 331 222 1482 361 708Tightening from 0.99 to 0.80 narrows every
window. That is a real trade-off rather than a free improvement: a
narrower window gives a better ergodicity bound and a smaller
m_convergence, but truncates more genuine posterior mass.
0.99 is the default because it is close to no truncation at
all.
The remaining attributes are diagnostics:
names(attributes(dl))
#> [1] "names" "window_diagnostics"
#> [3] "group.dispersion.fit" "calibration_source"
#> [5] "measurement_prior_group" "class"
str(attr(dl, "measurement_prior_group"), max.level = 1)
#> List of 4
#> $ shape_group : Named num [1:18] 1.55 1.55 1.55 1.55 1.55 ...
#> ..- attr(*, "names")= chr [1:18] "308" "309" "310" "330" ...
#> $ rate_group : Named num [1:18] 1502 43.5 85 274.8 305.9 ...
#> ..- attr(*, "names")= chr [1:18] "308" "309" "310" "330" ...
#> $ disp_lower_group: Named num [1:18] 1088.1 31.5 61.6 199 221.6 ...
#> ..- attr(*, "names")= chr [1:18] "308" "309" "310" "330" ...
#> $ disp_upper_group: Named num [1:18] 7277 211 412 1331 1482 ...
#> ..- attr(*, "names")= chr [1:18] "308" "309" "310" "330" ...Note what measurement_prior_group deliberately omits: no
mu and no Sigma. The Block 1 prior
mean is \(\mathcal{W}_j\gamma\), recomputed every
sweep from the current Block 2 state, and the RE covariance comes from
pfamily_list. Only the dispersion half of the ING prior is
fixed at calibration time.
The per-group route keeps the v1 joint structure: for each group,
each sweep, one call to rindepNormalGamma_reg() draws \((\beta_j, \sigma^2_j)\) together from the
ING conditional, dispatched through
.two_block_block1_ing_group_draw_one_chain().
A v2 partition is derived in
inst/DERIVATION_sigma2_with_block2_v2.md — draw \(\beta \mid (\gamma, \Psi, \Omega)\) as a
plain Gaussian, then
\[ \Omega_j \mid \beta_j \sim \mathrm{Gamma}\Bigl(a_j^0 + \frac{n_j}{2},\ r_j^0 + \frac{1}{2}\|y_j - D_j\beta_j\|^2\Bigr), \]
which would sidestep envelope rejection entirely. It is theory and stubs only; v1 remains the production route. Its cost is that it needs the extended \(P_{11}\) blocks for \(\Omega\) from Chapter B06 before the rate could be certified.
\(\Omega_j\) enters the ergodicity story through \(P_{22}\) alone:
\[ P_{22} = \mathrm{blockdiag}_j\bigl(D_j'\Omega_j D_j + \Psi^{-1}\bigr). \]
This is a structural fact worth isolating. The measurement residual \(e_j = y_j - D_j\beta_j\) depends on \(\beta_j\) only, never on \(\gamma\), so \(\Omega\) is Hessian-orthogonal to the Block 2 coordinates and contributes nothing to \(P_{11}\) or \(P_{12}\). Contrast \(\lambda_p\) from Chapter B06, which enters through the RE link \(u_p = \beta_{\cdot p} - \mathcal{W}_p\gamma_p\) and therefore touches all four blocks.
Because \(\lambda^*\) is
Loewner-non-increasing in \(\Omega\),
the corner is at \(\Omega_{\min}\),
i.e. at the largest \(\sigma^2_j\) the window allows —
disp_upper. That is what
.rLMM_measurement_disp_upper_for_rate() supplies to
two_block_rate_from_pfamily_list(), and it completes the
pair with Chapter B06’s disp_lower corner on \(\tau^2\).
Under the marginal route the plug-in becomes a Student-\(t\)-type Hessian with a rank-one correction,
\[ H_j^{(\Omega)}(\beta_j) = E_t[\Omega_j]\,D_j'D_j - \mathrm{Var}_t[\Omega_j]\,(D_j'e_j)(D_j'e_j)', \]
positive definite only inside the ellipsoid \(q_j(\beta_j) \le E_t[\Omega_j]/\mathrm{Var}_t[\Omega_j]\) where \(q_j = \|\beta_j - \hat\beta_j^{\mathrm{ols}}\|^2_{D_j'D_j}\). Tightening the window shrinks \(\mathrm{Var}_t[\Omega_j]\) and so raises the threshold, which is the third job the window does. See Chapter B06.
Worth stating plainly, since per-group dispersion is the most approximation-laden path in the package.
Empirical-Bayes double counting. Each group’s prior is calibrated from that group’s own data, so its data are used twice. The default \(w_j = 0.01\) keeps the effect at \(O(1/J)\) but does not remove it.
Calibration is conditional on Block 2. The marginal integrates \(\beta_j\) out given the current \(\mathrm{sd}_\tau\) and \(\Sigma_{\mathrm{fixef}}\); it is not joint over \(\tau^2\) uncertainty. This matches the two-block convention, but it is a conditional statement.
The glmmTMB reference fit drives
everything. When calibration_source is
"glmmTMB", that fit’s estimates are the
calibration. glmmTMB can converge to a boundary \(\tau^2 \approx 0\) where lmer
finds a positive estimate, and there is no automatic fallback. Comparing
against the pooled setup’s Sigma_ranef is the cheap
check:
round(rbind(pooled = diag(ps_pooled$group.Sigma),
group = diag(ps_group$group.Sigma)), 3)
#> (Intercept) Days_c
#> pooled 1408.730 35.072
#> group 1464.195 36.708Small groups. With \(n_j\) barely above \(p_{re}\), the RE covariance can be
near-singular and the ING envelope can converge poorly. The
group.alpha_target argument (default 0.01)
searches for the smallest per-group group.dispersion.pwt
meeting an ellipsoid-violation target, floored at the user’s supplied
weight; results land in ps$group.pwt_calibration.
Runtime diagnostics live in lmebayes.
summary_sigma2() is a method on fitted lmerb
objects in lmebayes, not part of this package. See
vignette("Chapter-08", package = "lmebayes").
vignette("Chapter-B05") — the truncation window for
\(\tau^2\).vignette("Chapter-B06") — where \(\Omega\) enters the extended rate.vignette("Chapter-08", package = "lmebayes") — the
user-facing view.?dGamma_list, ?Prior_Setup_GLMMinst/DGAMMA_LIST_MARGINAL_AND_BOUNDS.md,
inst/P_MATRIX_MARGINAL_PRECISIONS.md,
inst/DERIVATION_sigma2_with_block2_v2.md,
inst/PREFLIGHT_dGamma_list.md in the package sources.