Chapter B09: Row-block engines

A row-block engine partitions the rows of a design into groups and runs one regression per group. That description fits two quite different things in this ecosystem, and conflating them is the main hazard this chapter exists to prevent:

• rNormal_reg_group() / rNormalGLM_reg_group() — the compiled block engines in this package. They are what Block 1 of the mixed-model sampler calls, every sweep.
• lmbBlock() / glmbBlock() in lmebayes — independent “by group” fits with no population level. Despite the name, they do not use the engines above; they loop over blocks calling glmbayes::lmb().

Same shape, different code, different statistics.

The engines

args(rNormal_reg_group)
#> function (n, y, x, group, prior_list = NULL, prior_lists = NULL, 
#>     offset = NULL, weights = 1, Gridtype = 2L) 
#> NULL
args(rNormalGLM_reg_group)
#> function (n, y, x, group, prior_list = NULL, prior_lists = NULL, 
#>     offset = NULL, weights = 1, family = gaussian(), Gridtype = 2L, 
#>     n_envopt = NULL, use_parallel = TRUE, use_opencl = FALSE, 
#>     verbose = FALSE, progbar = FALSE) 
#> NULL

The Gaussian engine draws each block’s coefficients from its normal full conditional; the GLM engine uses the envelope sampler of Nygren and Nygren (2006) per block. Both accept either a single prior_list shared across blocks or a prior_lists giving one per block.

The compiled path partitions once and loops in C++:

rNormalGLM_reg_group()  →  .block_rNormalGLM_cpp()
                        →  block_rNormalGLM_cpp_export()   [src/block_utils.cpp]
                        →  rNormalGLMBlocks()              [src/rNormalGLMBlocks.cpp]
                        →  rNormalGLM()  per block         [glmbayesCore]

rNormal_reg_group()     →  .block_rNormalReg_cpp()
                        →  block_rNormalReg_cpp_export()
                        →  rNormalRegBlocks()
                        →  rNormalReg()  per block         [glmbayesCore]

The per-block draw itself is glmbayesCore’s. What this package adds is the partition, the per-block prior plumbing, and the loop — which is precisely the part that has to be fast, since Block 1 executes it once per sweep per chain.

A small worked example:

set.seed(1)
J <- 6; n_j <- 20
g <- factor(rep(seq_len(J), each = n_j))
x <- cbind("(Intercept)" = 1, slope = rnorm(J * n_j))
beta_true <- cbind(rnorm(J, 10, 2), rnorm(J, 2, 0.5))
y <- rowSums(x * beta_true[as.integer(g), ]) + rnorm(J * n_j)

fit <- rNormal_reg_group(
  n = 100, y = y, x = x, group = g,
  prior_list = list(mu = c(0, 0), Sigma = diag(c(100, 100)), dispersion = 1)
)
class(fit)
#> [1] "rNormal_reg_group" "list"
dim(fit$coefficients)
#> [1] 6 2
round(head(fit$coefficients, 3), 3)
#>   (Intercept) slope
#> 1       8.564 1.757
#> 2      12.589 1.808
#> 3      10.182 1.446

coefficients and coef.mode are \(k \times l_1\) matrices — one row per block, one column per predictor. This is worth internalizing, because it is the opposite of the draws-in-rows convention used elsewhere: here the block index is the row.

rownames(fit$coefficients)
#> [1] "1" "2" "3" "4" "5" "6"
colnames(fit$coefficients)
#> [1] "(Intercept)" "slope"

str(fit$group_info, max.level = 1)
#> List of 5
#>  $ k        : int 6
#>  $ ids      : chr [1:6] "1" "2" "3" "4" ...
#>  $ l2_blocks: int [1:6] 20 20 20 20 20 20
#>  $ starts   : int [1:6] 1 21 41 61 81 101
#>  $ rows     :List of 6

Block 1 of the mixed-model sampler

This is where the engines earn their keep. Inside rGLMM_sweep():

rGLMM_sweep
  └── .two_block_block1_all_chains()
        ├── .two_block_block1_prep_all_chains()
        │     └── build_mu_all(design, fixef_i)        # γ → per-group prior means
        │     └── .two_block_block1_prior_with_tau2()  # τ² → prior precision
        └── .two_block_block1_draw_all_chains()
              └── .two_block_block1_draw_block()
                    └── rNormal_reg_group() / rNormalGLM_reg_group()

The prep step is the whole difference between a Gibbs sweep and a set of independent fits. Each sweep, for each chain:

fixef_i    <- .two_block_batch_fixef_chain(batch, i)
mu_all     <- build_mu_all(design, fixef_i, batch$group_levels)$mu_all
prior_list <- .two_block_block1_prior_with_tau2(block1_prior, tau2_i, ...)

build_mu_all() evaluates \(\mathcal{W}_j\gamma\) for every group at the current Block 2 draw and returns a \(p_{re} \times J\) matrix. That matrix becomes prior_list$mu, in which column \(j\) is block \(j\)’s prior mean, matched by position. Separately, .two_block_block1_prior_with_tau2() refreshes the diagonal of the prior precision, P[k, k] <- 1/tau2_k, from the current \(\tau^2\).

So both halves of each block’s prior — mean and precision — are rebuilt from the population state before every draw. That is the coupling; the engine itself is stateless.

Comparing the two uses side by side:

Two of these are easy to trip over when reading the source. n = 1L is right because the sweep loop supplies the iteration; asking for more draws per sweep would be sampling from a conditional that is about to change. And the Block 1 wrapper deliberately does not forward design$offset or design$weights — those are already absorbed upstream, so passing them again would double-count.

The reordering step matters because normalize_group() orders blocks by factor level while the sweep batch carries its own group_levels ordering; .two_block_block1_draw_reorder() reconciles them so that downstream code can index by position.

lmbBlock() and glmbBlock() are not these engines

Here is the part most worth stating explicitly. lmbBlock() has exactly one runtime dependency on this package:

lmebayesCore:::normalize_group(block_vec, l2)

called from .blmb_formula_block_meta(). That is all. The flow is:

1. .blmb_formula_block_meta() — build the model frame, resolve the block specification, call normalize_group().
2. .mrglmb_normalize_pfamily_lists() / .validate_pfamily_for_rlmb() — prior validation, from glmbayesCore.
3. For b = 1..k: subset the rows, call glmbayes::lmb() or glmbayes::glmb().
4. .blmb_assemble() / .bglmb_assemble() — collect the results.

No call reaches rNormal_reg_group(), rNormalGLM_reg_group(), or any compiled code in this package. The block wrappers borrow a partition helper and otherwise delegate wholesale to glmbayes.

The result is a blmb (or bglmb) object: a named list of lmb (or glmb) fits with call, formula, block, group_info, k, p, pred_names and pfamily_lists as attributes — plus family for the GLM case.

Why keep them separate at all? Because they answer different questions. The block wrappers give genuinely independent per-block posteriors, which is what you want when the blocks are not exchangeable or when you want a no-pooling baseline to compare a mixed model against. Block 1’s draws are conditional on the population state and are not marginally independent; they are one component of a Gibbs sweep, not a standalone posterior.

The user-facing treatment is vignette("Chapter-12", package = "lmebayes").

Per-block prior calibration

Prior_SetupGroup() is the row-block analogue of Prior_Setup_GLMM(): it calls glmbayesCore::Prior_Setup() once per block on that block’s rows and returns a named list of per-block results, classed "Prior_SetupGroup", with call, formula, group and group_info attached.

There is no second level here, so there is nothing to calibrate across blocks — each block’s prior is derived from its own data and nothing else. That independence is exactly what makes the identifiability requirement stricter, which is the next section.

Identifiability, without a second level

Chapter B08 established that a rank-deficient group in the hierarchical model is survivable: the population conditional supplies what the data cannot. Independent block fitting has no population conditional, so that repair is unavailable.

Consequently only level 1 applies, and it is not advisory. A block whose design matrix is rank deficient cannot identify its coefficients from its data, and under a weak prior the posterior is improper in the deficient direction. Such blocks should be dropped before fitting, not passed in and flagged.

block_check_identifiability() in lmebayes does the screening. It returns level1_table, level1_keep, level1_drop, level2_rank, level2_q, level2_ok and action, and takes on_failure = c("warn", "stop").

The level-2 fields are present for the coupled case rather than the independent one, and its failure message is the sharpest statement of why that condition is the fatal one:

block_check_identifiability: Level 2 FAILED -- rank of X_nbhd restricted to
Level-1-identified blocks is {l2_rank} (need {q}). The population parameter
mu is not identified by the data-bearing blocks; the coupled Gibbs chain
will be null recurrent in {q - l2_rank} direction(s) of mu even with a
proper prior.

“Even with a proper prior” is the operative phrase, and the contrast with level 1 — where a proper prior is the fix — is the whole content of Chapter B08’s asymmetry.

Who owns what

Since three packages are involved, the division of labor in one table:

The dependency worth noting: glmbayes has no dependency on lmebayesCore. The arrow runs the other way, and lmbBlock() reaching in for normalize_group() is the single exception in the whole stack.

See also

• vignette("Chapter-B01") — the sweep loop that drives Block 1.
• vignette("Chapter-B08") — the identifiability conditions.
• vignette("Chapter-12", package = "lmebayes") — lmbBlock() for users.
• ?rNormal_reg_group, ?rNormalGLM_reg_group, ?Prior_SetupGroup
• inst/BLOCK1_draw_partition_and_labeling.md, inst/LMBBLOCK_LMEBAYESCORE_DEPENDENCIES.md in the package sources.