library(lmebayes)

Chapter 12: Row-block models

lmbBlock() and glmbBlock() sit one level below the mixed-model interface. They partition the rows of a dataset into blocks and fit an independent Bayesian regression within each block — no hierarchy, no shared population level, no shrinkage between blocks.

That makes them the “no pooling” corner of Chapter 02, done properly: each block gets a full posterior rather than a point estimate, and the fits are independent by construction rather than by accident.

When you want this

As a baseline. Fitting the same data with and without a hierarchy is the clearest way to see what partial pooling actually did.

When the blocks are not exchangeable. The mixed model assumes group coefficients are draws from a common population. If your blocks are different experiments, different instruments, or different time periods with no reason to share information, that assumption is wrong and pooling them is worse than not.

As a conceptual sibling of Block 1. A Block 1 update in lmerb() is also a set of independent per-group regressions — the difference is that there the per-group priors come from the population level and are refreshed every sweep, whereas here they are fixed once. The two are not the same code path, though: lmbBlock() fits each block with glmbayes::lmb(), while Block 1 uses the compiled block engines in lmebayesCore. vignette("Chapter-B09", package = "lmebayesCore") traces both.

A worked example

data(sleepstudy, package = "lme4")
dat <- sleepstudy
dat$Days_c <- dat$Days - mean(dat$Days)

## Three subjects, to keep the output readable.
sub <- droplevels(subset(dat, Subject %in% c("308", "309", "310")))
table(sub$Subject)
#> 
#> 308 309 310 
#>  10  10  10

Priors come from Prior_SetupGroup(), the row-block counterpart of Prior_Setup_GLMM(). Note the argument is group:

psb <- Prior_SetupGroup(Reaction ~ Days_c, group = sub$Subject,
                        data = sub, pwt = 0.05)
psb
#> 
#> Call:
#> Prior_SetupGroup(formula = Reaction ~ Days_c, group = sub$Subject, 
#>     data = sub, pwt = 0.05)
#> 
#> --- Row-block prior setup ---
#>   formula : Reaction ~ Days_c
#>   blocks  : 3 (308, 309, 310)
#>   Each element is a Prior_Setup() result; build pfamilies with
#>   pfamily_list(). Use print(x, blocks = ...) for details.

The object is a list of ordinary Prior_Setup() results, one per block. Each block is calibrated against its own data — there is no borrowing, which is exactly the point.

pfb <- pfamily_list(psb)
length(pfb)
#> [1] 3
names(pfb)
#> [1] "308" "309" "310"

set.seed(1)
fitb <- lmbBlock(Reaction ~ Days_c, block = sub$Subject, data = sub,
                 pfamily_list = pfb, n = 200L)

class(fitb)
#> [1] "blmb"
names(fitb)
#> [1] "308" "309" "310"

The result is a blmb object: a named list of lmb fits, one per block, with call, formula, block, group_info, k, p, pred_names, and pfamily_lists attached as attributes. Each element is an ordinary glmbayes lmb fit — the block wrapper partitions rows and loops, and lmb() does the sampling — so everything you know about lmb() applies:

summary(fitb[["308"]])
#> Call
#> lmb(formula = Reaction ~ Days_c, subset = 1:10, n = 200L, data = sub) 
#> 
#> Expected Residuals:
#>         Min          1Q      Median          3Q         Max 
#> -104.739754  -10.239600    9.410135   24.935573   61.783097 
#> 
#> Prior and Maximum Likelihood Estimates with Standard Deviations
#> 
#>             Null Mode Prior Mean Prior.sd Max Like. Like.sd
#> (Intercept)    342.13     342.13    69.29    342.13   15.11
#> Days_c           0.00       0.00    24.12     21.76    5.26
#> 
#> Bayesian Estimates Based on 200 iid draws
#> 
#>             Post.Mode Post.Mean  Post.Sd MC Error Pr(Null_tail) SE(tail)
#> (Intercept)  342.1338  342.4117  15.1242   1.0694        0.4850    0.035
#> Days_c        20.6765   20.9907   5.2055   0.3681        0.0000    0.000
#>             Pr(Prior_tail)    
#> (Intercept)          0.485    
#> Days_c              <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>   Directional Tail Summaries:
#> 
#>                Metric vs Null vs Prior
#>  Mahalanobis Distance  4.0841   4.0841
#>      Tail Probability  0.0000   0.0000
#>   [Tail probabilities are P(delta^T * Z <= 0) in whitened space]
#> 
#> 
#> Distribution Percentiles
#> 
#>                1.0%    2.5%    5.0%  Median   95.0%   97.5%  99.0%
#> (Intercept) 306.647 309.088 315.223 341.464 366.589 368.801 374.51
#> Days_c        9.875  11.795  12.484  20.604  30.191  31.373  32.59
#> 
#> Effective Number of Parameters: 1.780797 
#> Expected Residual Deviance: 9.027449 
#> DIC: 107.5355 
#> 
#> Expected Mean dispersion: 2527.151 
#> Sq.root of Expected Mean dispersion: 50.27078 
#> 
#> Mean Likelihood Subgradient Candidates Per iid sample: 1

Comparing blocks

Because the elements are ordinary fits, comparison is a loop:

tab <- t(vapply(fitb, function(f) colMeans(f$coefficients), numeric(2)))
round(tab, 3)
#>     (Intercept) Days_c
#> 308     342.412 20.991
#> 309     214.982  2.099
#> 310     231.014  5.701

round(t(vapply(fitb, function(f) apply(f$coefficients, 2, sd), numeric(2))), 3)
#>     (Intercept) Days_c
#> 308      15.124  5.206
#> 309       3.045  1.032
#> 310       4.223  1.447

Versus the hierarchical fit

The same three subjects, fitted hierarchically:

form <- Reaction ~ Days_c + (1 + Days_c || Subject)
ps <- Prior_Setup_GLMM(form, data = sub, pop.pwt = 0.05)
set.seed(1)
fith <- lmerb(form, data = sub,
              pfamily_list     = pfamily_list(ps),
              dispersion_ranef = ps$group.dispersion,
              n = 200L, progbar = FALSE)

hier <- as.matrix(lme4::ranef(fith)$Subject)
round(cbind(
  block_intercept = tab[, "(Intercept)"],
  hier_intercept  = hier[rownames(tab), "(Intercept)"],
  block_slope     = tab[, "Days_c"],
  hier_slope      = hier[rownames(tab), "Days_c"]
), 3)
#>     block_intercept hier_intercept block_slope hier_slope
#> 308         342.412        340.125      20.991     20.421
#> 309         214.982        215.803       2.099      3.027
#> 310         231.014        231.926       5.701      6.372

The hierarchical estimates sit between the block estimates and their common average — that is shrinkage, visible directly. With only three subjects and ten observations each, the pull is substantial. It would be weaker with more data per subject and stronger with less, following the ratio in Chapter 02.

Which is right depends on a question the data cannot answer: are these three subjects draws from a population you care about? If yes, the hierarchical estimates are better — they borrow strength. If they are three unrelated experiments, the block estimates are the honest ones.

Generalized responses

glmbBlock() is the same idea with a family argument:

psg <- Prior_SetupGroup(y ~ x, group = dat$block, data = dat,
                        family = poisson(), pwt = 0.05)
fitg <- glmbBlock(y ~ x, block = dat$block, data = dat,
                  family = poisson(),
                  pfamily_list = pfamily_list(psg), n = 200L)

Each block runs glmbayes::glmb() independently, so the draws are iid within block and independent across blocks. There are no Gibbs sweeps anywhere: without a shared population level there is nothing to alternate between, and consequently no m_convergence, no tv_tol, and no pilot stage. Convergence is not a question that arises.

What you give up

	`lmbBlock()`	`lmerb()`
Blocks share information	No	Yes
Population-level parameters	None	Estimated
Small / noisy blocks	Estimated on their own	Shrunk toward the population
Group-level covariates	Not expressible	Level-2 design
Sampling	iid within block	iid draws via certified sweeps
Blocks with rank-deficient designs	Fail individually	Excluded, others unaffected

The most consequential row is the third. A block with two observations gets a posterior dominated by its prior, and that prior was calibrated from those same two observations. The hierarchical model handles this case properly by letting the other blocks inform it.

Chapter 12: Row-block models

Chapter 12: Row-block models

When you want this

A worked example

Comparing blocks

Versus the hierarchical fit

Generalized responses

What you give up

See also