Chapter B01: The two-block Gibbs architecture

Chapter-B00 introduced the two-stage model. This chapter follows one draw all the way down: from the user-facing call, through the replicate chain driver and the inner sweep loop, into the block updates and the glmbayesCore iid samplers where the randomness actually happens.

Why two blocks

The parameter vector splits naturally into

• Block 1 — the group-level coefficients \(\beta_1, \dots, \beta_J\) (dimension \(J \cdot P\)), and
• Block 2 — the population coefficients \(\gamma\) (dimension \(q\)), together with any sampled variance components \(\tau^2_k\).

Given \(\gamma\) and \(\Psi\), the \(\beta_j\) are conditionally independent across groups: group \(j\)’s conditional depends only on \(y_j\), \(D_j\), and \(\mathcal{W}_j\gamma\). Given all \(\beta_j\), the Block 2 conditional is a plain normal regression of the \(\beta\)’s on the level-2 designs, and because \(\Psi\) is diagonal it factorizes further across the \(P\) coefficients.

So a “two-block” sweep is really \(J + P\) small independent draws, never a \(J P + q\) dimensional one. The Block 1 precision matrix, which would be \(J p_{re} \times J p_{re}\), is never formed.

The price is autocorrelation between successive sweeps. The whole of Chapter-B02 and Chapter-B03 is about pricing that autocorrelation exactly, so the sampler can decide how many sweeps to discard between stored draws.

One sweep

for each inner sweep m = 1, ..., m_convergence:

  Block 1:  for j = 1, ..., J
              beta_j  ~  p(beta_j | y_j, gamma, Psi, sigma2_j)

  Block 2:  for k = 1, ..., P
              gamma_k ~  p(gamma_k | beta_.k, X_k, tau2_k)
              [ING only: also draw tau2_k, refresh Psi]

Two conventions matter downstream:

1. Block 2 is updated second. The total-variation theory of Chapter-B02 bounds the distance for the block updated second in each sweep. The stored Block 1 draw therefore lags by a half-step, which is why the sweep calibration adds one extra sweep.
2. The sampled \(\tau^2_k\) feeds back into Block 1. When a Block 2 component is dIndependent_Normal_Gamma, the drawn \(\tau^2_k\) replaces the corresponding diagonal entry of \(\Psi\) before the next Block 1 update.

Replicate chains, not one long chain

lmebayesCore does not produce one long correlated chain and thin it. It runs n independent replicate chains, each started at the exact joint posterior mean (Gaussian case) or the posterior mode (GLMM case), runs each for m_convergence inner sweeps, and stores one draw per chain.

This is the same output contract as the iid samplers in glmbayesCore: the n stored draws are mutually independent, so ordinary Monte Carlo standard errors apply and no effective-sample-size correction is needed. The cost is that every stored draw pays the full m_convergence burn-in; the benefit is that m_convergence is small (typically single digits) because the chain starts at the mean and the TV bound is explicit.

Starting every chain at the posterior mean also kills the mean term of the convergence bound: \(D_0 = 0\), so only the variance-convergence sum remains. See Chapter-B03.

The call stack

lmerb() / glmerb()                        [lmebayes: formula interface]
  |
  +-- model_setup()                       formula -> y, Z, group, X_hyper
  +-- Prior_Setup_GLMM()                  reference lmer/glmer -> prior scales
  +-- pfamily_list() / dGamma_list()      Block 2 and dispersion priors
  |
  +-- rlmerb() / rglmerb()                [lmebayesCore: matrix level]
        |
        +-- matrix_args_lmm()             design + priors -> route arguments
        +-- lmerb_posterior_mean()        exact Gaussian mean (ICM), or
        |   glmerb_posterior_mode()       GLMM mode -- also the chain start
        |
        +-- route table
              |
              +-- rLMMNormal_reg_known_vcov_iid()      exact iid, no sweeps
              +-- rLMMNormal_reg_known_vcov_two_bg()   two-block Gibbs
              +-- rLMMindepNormalGamma_reg_*()         ING sigma2 routes
              +-- rGLMM_reg_known_vcov()               non-Gaussian, fixed tau2
              +-- rGLMM_reg_estimated_vcov()           non-Gaussian, ING tau2
                    |
                    +-- rGLMM_sweep()     outer sweep loop (pilot and main)
                          |
                          +-- Block 1 across all chains
                          +-- Block 2 across all chains

rlmerb() and rglmerb() are the matrix-level entry points: they take a design object plus prior objects rather than a formula, and they are what lmebayes calls. Everything above them is formula parsing; everything below is sampling.

Sweep-outer versus chain-outer

Two loop orders exist in the codebase.

Sweep-outer (rGLMM_sweep(), the reference implementation): for each inner sweep m, all chains do Block 1, then all chains do Block 2.

for (m in seq_len(inner_sweeps)) {
  batch <- block1_all_chains(batch, ...)
  batch <- block2_all_chains(batch, ...)
  sweep_stats[[m]] <- snapshot_fixef_stats(batch)
}

Chain-outer (the older C++ driver): each chain runs all of its sweeps before the next chain starts.

Sweep-outer is preferred for three reasons. It exposes the natural parallelism — at fixed m the chains are independent within a block, so Block 1 across chains is embarrassingly parallel. It amortizes per-sweep setup (family functions, prior payloads) across chains rather than rebuilding it per chain. And it makes sweep_history meaningful: cross-chain means and variances at sweep m are a genuine snapshot of the sampler’s distribution after m sweeps, which is exactly what plot_mean_convergence() and plot_var_convergence() plot (Chapter-B10).

Some Gaussian fixed-\(\sigma^2\), fixed-\(\tau^2\) paths still use the chain-outer driver and therefore do not populate sweep_history; that gap is tracked in the package README.

Inside Block 1

Block 1 is split into a prep phase and a draw phase.

Prep assembles, per chain: the current Block 2 state \(\gamma\), the implied observation-level prior means \(\mathcal{W}_j\gamma\) (built by build_mu_all()), and the Block 1 prior precision \(\Psi^{-1}\) — refreshed from the current \(\tau^2\) draw when Block 2 is ING.

Draw then calls, once per chain:

• rNormal_reg_group() for Gaussian responses — a conjugate normal draw per group, and
• rNormalGLM_reg_group() for Poisson, binomial, and Gamma responses — the envelope accept-reject sampler per group, via glmbayesCore.

Both are row-block engines: they take the stacked design, a group index, and per-group priors, and return one draw per group. The same engines back lmebayes’s lmbBlock() / glmbBlock() interface, where the row partition is user-specified rather than a grouping factor (Chapter-B09).

Inside Block 2

Block 2 does three things per chain.

1. Align. The Block 1 draws are stored in group_levels order; the level-2 design X_hyper[[k]] has its own row order. The response passed to the Block 2 regression must follow rownames(X_hyper[[k]]). Skipping this alignment is silently correct on datasets where the two orders coincide and silently wrong everywhere else, so it is done unconditionally.

2. Draw \(\gamma_k\) for each random-effect coefficient \(k\), by calling glmbayesCore::rglmb() with pfamily_list[[k]]:

   • dNormal — \(\tau^2_k\) is the known regression dispersion, and the draw is a conjugate multivariate normal.
   • dIndependent_Normal_Gamma — a Gamma prior sits on the precision \(1/\tau^2_k\), and the sampler makes a joint \((\gamma_k, \tau^2_k)\) draw through the likelihood-subgradient envelope (Nygren and Nygren 2006; Griffin and Brown 2010). This is the route documented in Chapter-B05.

3. Write back. The new \(\gamma_k\) goes into the chain’s fixed-effect state; a new \(\tau^2_k\) (ING only) goes into \(\Psi\) for the next Block 1 prep.

The response for the Block 2 regression of coefficient \(k\) is the vector of \(J\) Block 1 draws for that coefficient. So Block 2 sees \(J\) “observations” — one per group — regardless of how many rows the original data has. This is why the prior-versus-data guard for ING components requires the effective prior sample size to be at most \(J\).

Where C++ takes over

The compiled entry points resolve their R-level callbacks through the glmbayesCore namespace rather than a local copy, so there is exactly one implementation of the envelope machinery in the ecosystem even though two packages link against it.

Optional OpenCL offload for envelope construction is inherited from the same layer; see vignette("Chapter-A10", package = "glmbayes").

What one sweep stores

With inner_sweeps = 1 the first sweep_history row equals colMeans(fixef_draws) exactly — a useful invariant when checking a new driver against an old one. For parity work between R and C++ kernels, compare cross-chain means at large n rather than expecting bit-identical RNG streams: the two implementations draw from different streams.

See also

• Chapter B02 — why the sweep loop converges, and how fast.
• Chapter B03 — how m_convergence is chosen.
• Chapter B09 — the row-block engines used inside Block 1.
• Chapter B10 — reading sweep_history.
• inst/ARCHITECTURE_glmerb.md in the package sources, the maintainer-level file-by-file map this chapter summarizes.