Chapter B00: Engine roadmap and notation

library(lmebayesCore)
#> Loading required package: MASS

Chapter B00: Engine roadmap and notation

lmebayesCore is the engine package behind lmebayes. It owns the formula-to-design pipeline (model_setup()), the Block 2 hyperprior calibration (Prior_Setup_GLMM()), the matrix-level two-block Gibbs samplers (rlmerb(), rglmerb(), and the rLMM_reg / rGLMM_reg route families), and the convergence theory that decides how many inner Gibbs sweeps each stored draw needs.

This B-series documents that engine. It is written for maintainers, method developers, and readers who want the derivations behind the numbers that lmebayes prints. If you want to fit models, start with the lmebayes curriculum (Chapter-00 through Chapter-13) instead.

Two vignette series, two audiences

Series Package Audience
Chapter-00 … Chapter-13 lmebayes Applied users fitting mixed models
Chapter-B00 … Chapter-B10 lmebayesCore Engine internals, theory, derivations
Chapter-00 … Chapter-18, Chapter-A01 … A12 glmbayes The iid GLM stack this builds on

The Chapter-A* files that ship in the glmbayes and glmbayesCore sources describe the iid envelope machinery — likelihood subgradient densities (Nygren and Nygren 2006), envelope construction, and the pfamily prior system. Every Block 1 and Block 2 draw made by this package ultimately calls into that machinery through glmbayesCore. Read glmbayes Chapter-A05 (simulation methods) and Chapter-A08 (envelope functions) if the accept-reject layer is unfamiliar.

The B-series

  • B00 — Engine roadmap and notation (this chapter). Symbol glossary and the map between the internal parameterization and lme4.
  • B01 — The two-block Gibbs architecture. Block 1 / Block 2, the sweep-outer driver, replicate chains, and where C++ takes over.
  • B02 — Ergodicity and total-variation convergence. The rate matrix \(A = P_{11}^{-1/2} P_{12} P_{22}^{-1} P_{21} P_{11}^{-1/2}\), the full-rank conditions, and why the chain is geometrically ergodic.
  • B03 — Calibrating m_convergence. From \(\lambda^*\) to a certified sweep count, and the pilot-versus-main cost trade-off.
  • B04 — The known-vcov Gaussian route and the exact iid sampler.
  • B05 — Estimated variance components (ING Block 2).
  • B06 — The extended local rate diagnostic.
  • B07 — Per-group measurement dispersion.
  • B08 — Design, identifiability and rank conditions.
  • B09 — Row-block engines.
  • B10 — Sweep history and convergence plots.

The model

For group \(j = 1, \dots, J\) with \(n_j\) observations and \(P\) group-varying coefficients, the engine works with a centered (Lindley–Smith) parameterization (Lindley and Smith 1972):

Stage 1 (likelihood). \[ y_j \mid \beta_j, \sigma^2_j \;\sim\; N\!\big(D_j \beta_j,\; \sigma^2_j I_{n_j}\big) \]

Stage 2 (hierarchical prior). \[ \beta_j \mid \gamma, \Psi \;\sim\; N\!\big(\mathcal{W}_j \gamma,\; \Psi\big) \]

The full, non-centered coefficient vector \(\beta_j\) — not a mean-zero deviation — enters the likelihood directly. This is the single most important structural difference from lme4, and it is what makes the two-block split natural: Block 1 draws every \(\beta_j\) given \(\gamma\), and Block 2 draws \(\gamma\) given all \(\beta_j\).

For non-Gaussian families the stage-1 normal is replaced by the corresponding log-concave GLM likelihood with link \(g\), and the Block 1 draw goes through the envelope sampler rather than a conjugate normal.

Symbol table

Symbol Name Dimension Role
\(y_j\) response, group \(j\) \(n_j \times 1\) data
\(D_j\) likelihood design, group \(j\) \(n_j \times P\) multiplies \(\beta_j\)
\(\beta_j\) group-level coefficients \(P \times 1\) not mean zero; Block 1 state
\(\sigma^2_j\) measurement (observation) dispersion scalar per-group or pooled
\(\mathcal{W}_j\) level-2 design, group \(j\) \(P \times q\) block-diagonal; \(\gamma \to E[\beta_j]\)
\(\gamma\) population / hyper-mean coefficients \(q \times 1\) Block 2 state
\(\Psi\) covariance of \(\beta_j\) about \(\mathcal{W}_j\gamma\) \(P \times P\) diagonal, \(\Psi = \mathrm{diag}(\tau^2_k)\)
\(\tau^2_k\) variance component \(k\) scalar known (dNormal) or sampled (ING)
\(u_j := \beta_j - \mathcal{W}_j\gamma\) derived deviation \(P \times 1\) mean zero; not sampled directly

\(\mathcal{W}_j\) is block-diagonal across the \(P\) coefficient dimensions, \[ \mathcal{W}_j = \mathrm{blockdiag}\big(W_{1j}, \dots, W_{Pj}\big), \qquad W_{kj} \in \mathbb{R}^{1 \times q_k}, \quad q = \sum_{k=1}^{P} q_k , \] so each random-effect coefficient may carry its own set of level-2 predictors. \(q_k = 1\) (intercept only) is the default when a coefficient has no group-level covariates. This raggedness is exactly why the level-2 designs are stored as a list X_hyper with one \(J \times q_k\) matrix per coefficient rather than as a single padded array.

Internal storage

Quantity Storage Shape
Stacked likelihood design Z (matrix argument x) \(n \times P\)
Group index group length \(n\), \(g_i = j\)
Level-2 designs X_hyper = list(X_1, …, X_P) X_hyper[[k]] is \(J \times q_k\)
Dimension bookkeeping q_k = ncol(X_hyper[[k]]), q = sum(q_k) —

Recovery rules: \(D_j\) is Z[group == j, ], and \(\mathcal{W}_j\) is blockdiag(X_1[j, ], …, X_P[j, ]), assembled on demand. The design is stored row-stacked, never as the block-diagonal \(n \times PJ\) object lme4 implicitly uses for \(Z\).

Full conditionals

Block 1 (\(\beta_j\), conjugate normal, one draw per group): \[ \beta_j \mid \cdot \;\sim\; N(m_j, V_j), \qquad V_j = \Big(\tfrac{1}{\sigma^2_j} D_j' D_j + \Psi^{-1}\Big)^{-1}, \qquad m_j = V_j\Big(\tfrac{1}{\sigma^2_j} D_j' y_j + \Psi^{-1}\mathcal{W}_j\gamma\Big). \]

Block 2 (\(\gamma\), generalized least squares across groups): \[ \hat\gamma = \Big(\sum_j \mathcal{W}_j' \Psi^{-1} \mathcal{W}_j\Big)^{-1} \sum_j \mathcal{W}_j' \Psi^{-1} \beta_j . \]

Because \(\Psi\) is diagonal and \(\mathcal{W}_j\) is block-diagonal, the Block 2 update factorizes across coefficients \(k\): each \(\gamma_k\) is drawn from its own \(q_k\)-dimensional normal regression of \(\beta_{\cdot k}\) on \(X_k\). That is why Block 2 is specified as a pfamily_list — one prior object per random-effect coefficient — rather than as one joint \(q\)-dimensional prior.

Correspondence to lme4

Substituting stage 2 into stage 1 recovers the familiar marginal form (Bates et al. 2015): \[ y_j = D_j(\mathcal{W}_j\gamma + u_j) + \varepsilon_j = \underbrace{(D_j\mathcal{W}_j)}_{X_j}\gamma + \underbrace{D_j}_{Z_j} u_j + \varepsilon_j . \]

lmebayesCore lme4 / Laird–Ware Relationship
\(D_j\) — ancestor of both \(X_j\) and \(Z_j\)
\(D_j \mathcal{W}_j\) \(X_j\) fixed-effects design block
\(D_j\) \(Z_j\) random-effects design block
\(\gamma\) \(\beta\) fixed effects
\(u_j = \beta_j - \mathcal{W}_j\gamma\) \(b_j\) random effects (mean zero)
\(\Psi\) per-group block of \(G\) \(G = I_J \otimes \Psi\)
\(\beta_j\) \(\hat\beta + \hat b_j\) lme4::coef(), when \(\mathcal{W}_j\) is intercept-only

Note the asymmetry: \(D_j\) has no single lme4 counterpart. When \(\mathcal{W}_j\) is intercept-only for every coefficient (\(q_k = 1\ \forall k\)), the model reduces to a standard random intercept / random slope lmer fit, and the posterior mean of \(\gamma\) should track lme4::fixef(), \(\hat\beta_j - \hat\gamma\) should track lme4::ranef(), and \(\hat\beta_j\) should track lme4::coef(). The package test suite uses exactly this equivalence as a correctness check.

Model restrictions

The engine deliberately covers a narrower model class than lme4:

  • A single grouping factor. No crossed or nested random effects.
  • Uncorrelated random-effect terms. Formulas use the || form, so \(\Psi\) is diagonal. This is what lets Block 2 factorize across coefficients and lets the variance components be handled one at a time.
  • Log-concave likelihoods. Gaussian, Poisson, binomial, and Gamma with the standard links. Log-concavity is what the envelope accept-reject sampler requires (Nygren and Nygren 2006).

model_setup() enforces these and reports the resulting design in a printable object. check_identifiability() adds the two rank conditions derived in Chapter-B08.

Naming conventions in the code

Prefix / suffix Meaning
rLMM* Gaussian (linear) mixed-model replicate-chain routes
rGLMM* Non-Gaussian replicate-chain routes
*_known_vcov Variance components \(\tau^2_k\) treated as fixed
*_estimated_vcov \(\tau^2_k\) sampled (ING Block 2)
*_iid Exact independent draws, no Gibbs sweeps
two_block_* Convergence theory, rates, TV bounds, pilot costing
.two_block_* Internal (unexported) helpers of the above
block_*, *_reg_group Row-block engines (Chapter-B09)

Names beginning with a dot, and topics marked \keyword{internal}, are not exported; reach them with lmebayesCore::: when reproducing the diagnostics in these chapters.

See also

  • Chapter B01 — the two-block architecture in detail.
  • Chapter B02 — ergodicity and total-variation convergence.
  • vignette("Chapter-A05", package = "glmbayes") — the iid simulation methods every block draw is built from.
  • inst/notation.md in the package sources, from which this chapter’s symbol table is derived.

References

Bates, Douglas, Martin Mächler, Ben Bolker, and Steve Walker. 2015. “Fitting Linear Mixed-Effects Models Using lme4.” Journal of Statistical Software 67 (1): 1–48. https://doi.org/10.18637/jss.v067.i01.
Lindley, D. V., and A. F. M. Smith. 1972. “Bayes Estimates for the Linear Model.” Journal of the Royal Statistical Society. Series B (Methodological) 34 (1): 1–41. https://doi.org/10.1111/j.2517-6161.1972.tb00899.x.
Nygren, K. N., and L. M. Nygren. 2006. “Likelihood Subgradient Densities.” Journal of the American Statistical Association 101 (475): 1144–56. https://doi.org/10.1198/016214506000000357.