Chapter B04 held every \(\tau^2_k\)
fixed. This appendix documents what happens when a component’s prior is
dIndependent_Normal_Gamma (ING) instead of
dNormal, so that \(\tau^2_k\) is sampled alongside the
population coefficients \(\gamma_k\).
The route becomes lmm_fixed_estimated, implemented by
rLMMNormal_reg_estimated_vcov() with
needs_pilot = TRUE.
The key re-framing is that, conditional on the current Block 1 state, each random-effect component \(k\) is an ordinary Gaussian regression whose response is the vector of group coefficients.
Fix component \(k\). Let \(b_{jk}\) be group \(j\)’s coefficient, \(X_k\) the \(J
\times p_k\) level-2 design (design$W[[k]]), and
\(\gamma_k\) the level-2 coefficients.
Then
\[ b_{jk} \sim N\!\bigl(X_k[j,\cdot]^\top \gamma_k,\ \tau^2_k\bigr), \qquad j = 1, \dots, J, \]
or in stacked form \(\mathbf{b} = X_k\gamma_k + \varepsilon\) with \(\varepsilon \sim N(0, \tau^2_k I_J)\).
That is a linear model with \(J\) observations, and \(\tau^2_k\) is its residual variance. Once you see it this way, the ING prior is just the standard independent normal–gamma prior for a regression:
\[ \gamma_k \sim N(\mu_\gamma, \Sigma_\gamma), \qquad \lambda_k = 1/\tau^2_k \sim \mathrm{Gamma}(\mathrm{shape}, \mathrm{rate}), \]
with the two independent a priori — which is what distinguishes ING from the conjugate normal-gamma, where \(\Sigma_\gamma\) would be scaled by \(\tau^2_k\).
Alignment matters: the order of b_k must match the rows
of X_k, which
two_block_align_b_to_xhyper(b_k, X_k, group_levels)
enforces.
Note what does not enter. mu_all is a
Block 1 object; the \(\tau^2\) update
needs only the residual \(\mathbf{b} -
X_k\gamma_k\), which is the same as subtracting row \(k\) of mu_all. Nor does
another component’s state appear — a diagonal \(\Psi\) is exactly what makes Block 2
factorize across \(k\).
Freezing \(\gamma_k\) reduces the problem to conjugate Gamma–Normal algebra. With \(\eta_j = X_k[j,\cdot]^\top\gamma_k\) and \(\mathrm{RSS} = \sum_j (b_{jk} - \eta_j)^2\), \(J\) independent normals of variance \(\tau^2_k\) contribute \(J/2\) to the shape and \(\mathrm{RSS}/2\) to the rate:
\[ \mathrm{shape}_{\mathrm{post}} = \mathrm{shape} + \frac{J}{2}, \qquad \mathrm{rate}_{\mathrm{post}} = \mathrm{rate} + \frac{\mathrm{RSS}}{2}, \]
and the implementation reports the mode
\[ \widehat{\tau^2_k} = \frac{\mathrm{rate}_{\mathrm{post}}} {\mathrm{shape}_{\mathrm{post}} - 1}. \]
Neither \(\mu_\gamma\) nor \(\Sigma_\gamma\) appears, because \(\gamma_k\) is held fixed. This path is not the production sampler; it is the algebra that Method B collapses to when \(\gamma_k\) is frozen, and it is what the ICM start uses.
rglmb(), the production pathThe Gibbs Block 2 update in two_block_block2_one_chain()
fits the full ING model rather than freezing \(\gamma_k\), by handing the regression to
glmbayesCore:
fit_k <- rglmb(
n = 1L,
y = b_al, # aligned group coefficients
x = X_k, # level-2 design
family = gaussian(),
pfamily = pf_ing # dIndependent_Normal_Gamma(mu, Sigma, shape, rate)
)
tau2_k <- fit_k$dispersion[1L]
gamma_k_new <- fit_k$coefficients[1L, ]So \(\gamma_k\) and \(\tau^2_k\) are updated jointly,
from their exact joint conditional, using the same envelope sampler
documented in
vignette("Chapter-A05", package = "glmbayes"). That joint
update is what keeps the two-block structure intact: Block 2 is still
one draw, not a sub-Gibbs alternation between \(\gamma\) and \(\tau^2\).
The cost shows up in Cand/draw in
summary(): the ING conditional is not conjugate, so each
draw needs several accept-reject candidates.
Before any iteration, \(\tau^2\) has to be set from the prior alone — the ICM start and the \(\Sigma_{\text{ranef}}\) used by Block 1 both need a number before there is any \(b\) to condition on:
\[ \tilde\tau^2_k = \frac{\mathrm{rate}}{\mathrm{shape} - 1} \qquad (\text{requires } \mathrm{shape} > 1), \]
computed by .two_block_tau2_plug_in_from_pfamily() and
.two_block_tau2_plug_in_vector(). The conservative rate
bounds use a different reference,
.two_block_tau2_ref_from_pfamily(), which returns \(\mathrm{rate}/\mathrm{shape}\) —
deliberately smaller, since smaller \(\tau^2\) means a slower chain.
Prior_Setup_GLMM() calibrates the prior so that the
plug-in reproduces the marginal hyper-regression center:
\[ \mathrm{shape}_k = \frac{n_0 + 1 + p_k}{2}, \qquad \mathrm{rate}_k = \hat\tau^2_k\,\frac{n_0 + p_k - 1}{2}, \]
where \(\hat\tau^2_k\) is the
marginal center from compute_gaussian_prior applied to
\(b_{\cdot k} \sim N(W_k\gamma_k, \tau^2_k
I_J)\) with \(\gamma_k \sim N(\mu_k,
\Sigma_k)\) — stored in
pop.prior_list[[k]]$dispersion — and \(n_0\) is the effective prior sample size
for dispersion, set by pop.dispersion.pwt. The classical
VarCorr plug-in is kept separately in
pop.dispersion.ref.
By construction \(\mathrm{rate}_k /
(\mathrm{shape}_k - 1) = \hat\tau^2_k\), so the ING prior and the
dNormal prior agree on the center and differ only
in whether the component is allowed to move. That is what makes the
ptypes switch a clean comparison: changing it does not
silently change the location of the prior.
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,
pop.dispersion.pwt = 0.2)
str(ps$pop.ing_prior)
#> List of 2
#> $ (Intercept):List of 5
#> ..$ shape : num 3.25
#> ..$ rate : num 3029
#> ..$ disp_lower : num 632
#> ..$ disp_upper : num 2998
#> ..$ max_disp_perc: num 0.99
#> $ Days_c :List of 5
#> ..$ shape : num 3.25
#> ..$ rate : num 66.1
#> ..$ disp_lower : num 13.8
#> ..$ disp_upper : num 65.4
#> ..$ max_disp_perc: num 0.99pop.ing_prior is one list per component, each carrying
the Gamma shape / rate on the precision plus
the truncation window.
ing <- ps$pop.ing_prior
round(data.frame(
plug_in = vapply(ing, function(p) p$rate / (p$shape - 1), numeric(1)),
ref_bound = vapply(ing, function(p) p$rate / p$shape, numeric(1)),
dNormal = vapply(ps$pop.prior_list, function(p) p$dispersion, numeric(1))
), 4)
#> plug_in ref_bound dNormal
#> (Intercept) 1346.1456 931.9470 1346.1456
#> Days_c 29.3722 20.3346 29.3722The plug-in matches the dNormal dispersion exactly; the
bound reference sits below it.
Each ING component carries disp_lower and
disp_upper, and the sampled \(\tau^2_k\) is confined to that interval.
The bounds are quantiles of the limiting posterior of
the variance component — not of the prior.
The distinction is the whole point. A Gamma prior on a precision puts mass arbitrarily close to \(\tau^2 = 0\), and \(\tau^2 \to 0\) is the degenerate limit in which every group coefficient is pinned to \(\mathcal{W}_j\gamma\): the hierarchy collapses, Block 1 and Block 2 become perfectly coupled, and the rate \(\lambda^*\) tends to 1. Truncating that region away is what keeps the chain geometrically ergodic, and the argument mirrors the drift-condition reasoning of Jones and Hobert (2001).
The lower bound does double duty. Because \(\lambda^*\) is monotone decreasing in \(\tau^2\) for this structure, evaluating the
rate at disp_lower gives a bound valid over the entire
window:
\[ \sup_{\tau^2 \in [\,\text{disp\_lower},\, \text{disp\_upper}\,]} \lambda^*(\tau^2) \;=\; \lambda^*(\text{disp\_lower}). \]
That is conservative rather than correct-on-average, which is the right direction for a certificate: the sweep count it produces is at least as large as any \(\tau^2\) actually visited would require.
The upper bound is less load-bearing but keeps the envelope
construction well conditioned and prevents a component from wandering
into a region the data cannot support.
inst/ING_TRUNCATION_WINDOW.md gives the quantile choices
and the derivation.
pf_ing <- pfamily_list(ps, ptypes = "dIndependent_Normal_Gamma")
design <- model_setup(form, data = dat)
set.seed(1)
fit <- rlmerb(n = 200L, design = design,
pfamily_list = pf_ing,
dispprior_list = list(dispersion = ps$group.dispersion),
progbar = FALSE, verbose = FALSE, print_icm_table = FALSE)c(engine = fit$convergence_info$sim_method_used,
m = fit$convergence_info$m_convergence)
#> engine m
#> "TWO_BLOCK_GIBBS" "3"
names(fit$sweep_history)
#> [1] "stage" "n_sweeps" "re_names" "fixef_mode" "table"
#> [6] "cov_by_sweep" "coef_index" "disp_table"
round(colMeans(fit$popef.dispersion), 3)
#> (Intercept) Days_c
#> 1400.245 30.481Two structural changes relative to Chapter B04.
sweep_history now has a pilot element as well
as main, and popef.dispersion holds draws
rather than a repeated constant.
The Theorem 3 bound of Chapter B02 has two terms: one for the distance of the starting mean from the target mean, and one for the variance. Chapter B04’s route sets the first to exactly zero by starting every chain at the exact posterior mean.
With an ING component there is no exact mean to start from. The joint
target over \((\gamma, \tau^2, b)\) is
not Gaussian, so the ICM start is a conditional mode, and mode \(\neq\) mean. The pilot stage produces a
start point whose distance from the truth can be bounded
instead: averaging n_pilot independent pilot chains gives a
start whose sampling covariance is \(\Sigma /
n_{\text{pilot}}\), so at confidence level \(\alpha\)
\[ D_0 = \frac{\chi^2_{p,\alpha}}{n_{\text{pilot}}}. \]
vignette("Chapter-B03") works through the optimization
that trades pilot sweeps against main sweeps.
ptypes accepts a single value or a named vector, so
components can be mixed:
pf_mixed <- pfamily_list(ps, ptypes = c(
"(Intercept)" = "dNormal",
"Days_c" = "dIndependent_Normal_Gamma"
))
vapply(pf_mixed, function(p) class(p)[1L], character(1))
#> (Intercept) Days_c
#> "pfamily" "pfamily"Because Block 2 factorizes across components, this is coherent: the
intercept’s \(\tau^2\) stays fixed
while the slope’s is sampled. The route key is still
lmm_fixed_estimated — pf_summary$all_dNormal
is FALSE as soon as any component is ING — so the
pilot stage and the conservative rate apply to the whole model. There is
no partial-exactness discount.
vignette("Chapter-B03") — pilot sizing and the cost
optimization.vignette("Chapter-B04") — the exact route this one
generalizes.vignette("Chapter-B06") — the extended rate diagnostic
for ING models.vignette("Chapter-07", package = "lmebayes") — the
user-facing view.inst/TAU2_ING_FORMULAS.md,
inst/ING_TRUNCATION_WINDOW.md,
inst/BLOCK_ING_RINDEPNORMALGAMMA_REG.md in the package
sources.