library(lmebayesCore)

.c_extdata_rds <- function(name) {
  p <- system.file("extdata", name, package = "lmebayesCore")
  if (nzchar(p) && file.exists(p)) return(p)
  for (root in c(
    file.path("..", "inst", "extdata", name),
    file.path("inst", "extdata", name)
  )) {
    if (file.exists(root)) return(normalizePath(root))
  }
  stop(
    "Cannot find ", name, " in inst/extdata/. ",
    "Run devtools::load_all() or reinstall the package.",
    call. = FALSE
  )
}

C04 <- readRDS(.c_extdata_rds("Chapter-C04-bwc-twoblock.rds"))
meta <- C04$meta
conv <- C04$convergence

This chapter fits the same model and prior specification as Chapter-C02, but forces the two-block Gibbs route instead of the exact iid sampler. Both chapters target the same multivariate normal posterior; here the draws come from independent chains, each run through 25 inner sweeps before being stored.

Prior specification

Two quantities are held constant throughout this chapter: the observation dispersion $\sigma^2$ and the between-group variance matrix $\Psi$. Both are fixed at the defaults Prior_Setup_GLMM() calibrates from the data, so neither carries a prior and neither is sampled:

C04$summary$tau2_prior_overview[, c("Prior", "E[tau2]")]
#>                   Prior   E[tau2]
#> (Intercept)     dNormal 159.13613
#> distracted_ppvt dNormal   6.73747
#> distracted_a1   dNormal  11.47851
#> Residual          fixed 388.80365

The three dNormal rows are the components of $\Psi$, declared known. The Residual row is $\sigma^2$, likewise fixed.

That leaves the population effects $\gamma$ as the only quantity with a prior of its own. They receive the default multivariate normal prior, specified separately for each group-effect component:

round(C04$summary$popef_prior_overview[, c("Prior Mean", "Prior.sd")], 3)
#>                                   Prior Mean Prior.sd
#> (Intercept)::(Intercept)              88.376   78.450
#> (Intercept)::private_school            0.000   96.853
#> (Intercept)::title1                    0.000  100.278
#> (Intercept)::free_reduced_lunch        0.000    1.571
#> distracted_ppvt::(Intercept)           0.000   17.506
#> distracted_a1::(Intercept)             0.000   41.419
#> distracted_a1::free_reduced_lunch      0.000    0.584

Every coefficient is centered at zero except the grand intercept, so the prior mean vector describes a null, intercept-only model — the same default as in Chapter-C02 and the glmbayes family of packages. The setup is built by two calls:

ps <- Prior_Setup_GLMM(form, data = dat, pop.pwt = 0.01)
pf <- pfamily_list(ps)

pop.pwt = 0.01 gives the prior on $\gamma$ the weight of about one percent of the data.

Estimating the model using the default settings

The call matches Chapter-C02 except for the simulation method:

fit <- rlmerb(
  n              = 10000,
  design         = model_setup(form, data = dat),
  pfamily_list   = pf,
  dispprior_list = list(dispersion = ps$group.dispersion),
  sim_method     = "TWO_BLOCK_GIBBS"
)

The design, pfamily_list and dispprior_list arguments carry the same meaning as in the previous chapter. sim_method = "TWO_BLOCK_GIBBS" forces the iterative route on this known-variance-components specification; with sim_method = "DEFAULT" the package would draw directly from the closed-form posterior instead.

Under this specification the joint posterior is exactly multivariate normal, as in Chapter-C02. Each stored draw is one independent chain after 25 inner Gibbs sweeps. That budget is not a guess: it is certified from the convergence rate (lambda_star = 0.9308) against a total variation tolerance, by the theory of vignette("Chapter-C03"). The derivation is in the section below.

Properties of the posterior and the two-block Gibbs procedure

The log posterior splits into the same three components as in Chapter-C02: the log-likelihood, the prior on the group effects, and the prior on the population effects. Because $\sigma^2$, $\Psi$ and $V$ are all known, the posterior is multivariate normal and factorizes into the same conditional distributions.

Block 1: group effects given the population effects

Given $\gamma$, the group effects are conditionally independent across groups:

\[ \beta_j \mid \gamma, y \sim N\!\left(b_j,\; \Sigma_j\right), \qquad \begin{aligned} \Sigma_j &= \left(\sigma^{-2} D_j' D_j + \Psi^{-1}\right)^{-1}, \\ b_j &= \Sigma_j\left(\sigma^{-2} D_j' y_j + \Psi^{-1}\mathcal{W}_j\gamma\right). \end{aligned} \]

Block 1 of the Gibbs sampler draws each $\beta_j$ from this conditional, holding the current $\gamma$ fixed.

Block 2: population effects given the group effects

Given the group effects, the population effects receive a conjugate normal update that combines the prior on $\gamma$ with the information each group’s current $\beta_j$ contributes through $\mathcal{W}_j$. The result is another multivariate normal conditional, $\gamma \mid \beta, y$, with precision that adds the prior precision to a sum of group-level contributions. Block 2 draws $\gamma$ from this conditional, holding the current $\beta$ fixed.

One sweep and one stored draw

One inner sweep alternates the two blocks: update all $\beta_j$, then update $\gamma$. Repeating the pair moves the chain closer to the target. After 25 sweeps the chain’s distribution is close enough that the final state can be treated as a draw from the posterior. Replicating that procedure across 10000 independent chains produces the stored sample.

Unlike the iid route of Chapter-C02, which draws $\gamma$ from its marginal and then each $\beta_j$ from its conditional in one pass, the Gibbs route needs the sweeps because it never integrates the group effects out analytically. Both routes target the same posterior; the Gibbs route trades the closed-form factorization for an iterative scheme whose convergence can be diagnosed from the sweep history.

Where the sweep count comes from

vignette("Chapter-C03") develops Nygren (2020)’s total variation bounds for a two-block Gibbs sampler on a multivariate normal target in the abstract. This model is exactly that setting, so the results apply verbatim and fix the sweep count that produced the draws above. This section instantiates them.

The rate matrix for this design

Take $\gamma$ — the 7-vector of population effects — as the first block, updated second in each sweep, and the group effects $\beta$ as the second block, updated first. Claim 2 says the mean bias contracts by the matrix $P_{11}^{-1}P_{12}P_{22}^{-1}P_{21}$ each sweep, and Remark 8 says the relevant eigenvalues are those of the symmetrized form $A = P_{11}^{-1/2}P_{12}P_{22}^{-1}P_{21}P_{11}^{-1/2}$ — which is why the rate can be read off the precision blocks without ever forming $\Sigma$.

Here $P_{11}$ is the prior precision of $\gamma$, $P_{22}$ is block diagonal with $j$th block $\sigma^{-2}D_j'D_j + \Psi^{-1}$, and $P_{12}$ collects the $-\Psi^{-1}\mathcal{W}_j$ coupling terms. two_block_rate() assembles these and solves the generalized eigenproblem:

C04$rate
#> Two-block Gibbs convergence rate (Nygren 2020, Remark 8)
#>   groups J = 34, p_re = 3, q = 7, n_obs = 400  [family: gaussian]
#> 
#> Eigenvalues of A (q x q):
#> [1] 0.930794 0.910041 0.865935 0.196770 0.191864 0.176490 0.172073
#> 
#> lambda* = 0.930794
#> 
#> Implied sweeps per tolerance:
#>   tol          (lambda*)^m proxy     Theorem 3 bound   Corollary 1 bound
#>   0.01                        65                  24                  45
#>   0.001                       97                  39                  61
#>   1e-06                      193                  87                 108
#>   (TV bounds assume start at the posterior mean: D0 = 0)

The spectrum splits sharply. Three eigenvalues sit near 0.87–0.93 and four near 0.17–0.20. The large ones belong to the intercept and the two slope directions, where the prior on $\gamma$ is weak relative to the information the group effects carry, so $\gamma$ and $\beta$ are strongly coupled and the sampler moves through that subspace slowly. The rate is $\lambda^{*} =$ 0.9308, the largest of them.

From the rate to a certified sweep count

Every chain starts at the exact posterior mean, so $D_0 = (x_1^{(0)}-\mu_1)'\Sigma_{11}^{-1}(x_1^{(0)}-\mu_1) = 0$ and the mean term of Theorem 3 vanishes identically. What remains is the variance sum $\sum_i d_i^{(l)}$, evaluated at the Remark 11 ratios $r_i^{(l)} = (1-a_{i-1}^{2l})/(1-a_i^{2l})$:

l <- c(1, 5, 10, 15, 20, 24, 30)
data.frame(
  sweeps     = l,
  theorem3   = signif(lmebayesCore:::two_block_tv_bound(C04$rate, l, "theorem3"), 3),
  corollary1 = signif(lmebayesCore:::two_block_tv_bound(C04$rate, l, "corollary1"), 3)
)
#>   sweeps theorem3 corollary1
#> 1      1  0.78500     1.0000
#> 2      5  0.24900     1.0000
#> 3     10  0.09200     1.0000
#> 4     15  0.03890     0.8850
#> 5     20  0.01740     0.4010
#> 6     24  0.00933     0.2150
#> 7     30  0.00375     0.0861

Both columns bound the same total variation distance, and the ordering between them is the price of Corollary 1’s closed form. Theorem 3 evaluates the error function differences exactly; Corollary 1 replaces each by its Remark 5 and Remark 17 envelope, which is loose at small $l$ — it exceeds 1, and so says nothing, until sweep 12 or so — but decays at the same geometric rate afterwards. The (lambda*)^m column of the printout above is neither: it is the naive proxy $(\lambda^{*})^{m}$, which ignores that the variance converges at $a_i^{2l}$ rather than $a_i^{l}$ and so is conservative by roughly a factor of two in sweeps.

At the package default tolerance of 0.01, Theorem 3 certifies 24 sweeps. The sampler adds one, giving 25:

c(
  theorem3      = lmebayesCore:::two_block_l_for_tv(C04$rate, 0.01, "theorem3"),
  stored        = C04$convergence$m_convergence
)
#> theorem3   stored 
#>       24       25

The extra sweep is a half-step correction, not slack. The bound is a statement about the block updated second in each sweep — $\gamma$ — and the stored Block 1 draw $\beta$ lags it by half a sweep, so one more sweep brings both blocks inside the tolerance.

Two things are worth noting about the cost. Because the bound is geometric, tightening the tolerance is cheap: the printed table above shows that going from 0.01 to $10^{-6}$ — four orders of magnitude — costs roughly three and a half times the sweeps, not four orders of magnitude more. And because $D_0 = 0$ removes the mean term outright, the entire budget here is paid for variance convergence; routes that cannot start at the posterior mean (non-Gaussian families, or sampled variance components) must certify $D_0$ from a pilot stage and carry the Lemma 1 mean term as well.

Summarizing the output

summary() collects the draws into the same tables as in Chapter-C02:

summary(fit, groups = c("30", "15"))

The output is organized by block: reference random effects, Block 2 variance components (fixed here, so posterior tables are suppressed), Block 2 population effects in three tables, and Block 1 group effects with optional per-group detail.

Two features of the population effect tables follow from the default pop.pwt = 0.01. The posterior means land close to the lmer estimates, displaced slightly toward the null prior means. Read on the scale that matters, the two fits agree throughout: every Bayesian estimate falls within a fifth of a standard error of its lmer counterpart.

The Pr(Prior_tail) column measures the posterior mass on the opposite side of the prior mean from the estimate — a one-sided statement about how far the data have moved each coefficient off its prior value. Because the prior means are the null model, that reference point is zero for every coefficient except the intercept, so for those rows the column tracks half the two-sided p-value the corresponding lmer fit would report.

The output below was produced by data-raw/make_Chapter-C04.R with n = 10000 and set.seed(1), then stored, so that building this vignette does not rerun the sampler.

C04$summary
#> Call
#> rlmerb(n = N_DRAWS, design = design, pfamily_list = pf, dispprior_list = list(dispersion = ps$group.dispersion), 
#>     progbar = TRUE, verbose = FALSE, print_icm_table = FALSE, 
#>     sim_method = "TWO_BLOCK_GIBBS")
#> 
#> Bayesian mixed model fit based on 10000 draws [two-block Gibbs, 25 sweeps/draw]
#> Formula: score_ppvt ~ private_school + title1 + free_reduced_lunch + distracted_ppvt + distracted_a1 + free_reduced_lunch:distracted_a1 + (1 + distracted_ppvt + distracted_a1 || school_id) 
#> Number of obs: 400,  groups: school_id, 34
#> 
#> Random effects (lmer reference):
#>           grp            var1 var2   vcov  sdcor
#> 1   school_id     (Intercept) <NA> 191.79 13.849
#> 2 school_id.1 distracted_ppvt <NA>  26.03  5.102
#> 3 school_id.2   distracted_a1 <NA>  41.05  6.407
#> 4    Residual            <NA> <NA> 366.41 19.142
#> 
#> === Block 2 variance components (tau^2_k) ===
#> 
#> Prior and lmer reference
#> 
#>                   Prior  E[tau2] sqrt(E[tau2]) disp_lower disp_upper     lmer
#> (Intercept)     dNormal 159.1361       12.6149         NA         NA 191.7850
#> distracted_ppvt dNormal   6.7375        2.5957         NA         NA  26.0342
#> distracted_a1   dNormal  11.4785        3.3880         NA         NA  41.0508
#> Residual          fixed 388.8036       19.7181         NA         NA 366.4058
#>                 lmer SD
#> (Intercept)     13.8486
#> distracted_ppvt  5.1024
#> distracted_a1    6.4071
#> Residual        19.1417
#> 
#>   (All variance components held fixed; nothing sampled.)
#> 
#> === Block 2 population effects (level-2 hyperparameters) ===
#> 
#> Prior Estimates with Standard Deviations
#> 
#>                                   Prior Mean  Prior.sd      lmer lmer.se
#> (Intercept)::(Intercept)            88.37564  78.45025  99.16637   7.885
#> (Intercept)::private_school          0.00000  96.85269  -0.16630   9.734
#> (Intercept)::title1                  0.00000 100.27761  23.24038  10.078
#> (Intercept)::free_reduced_lunch      0.00000   1.57052  -0.38866   0.158
#> distracted_ppvt::(Intercept)         0.00000  17.50644  -1.61211   1.759
#> distracted_a1::(Intercept)           0.00000  41.41908   1.46279   4.163
#> distracted_a1::free_reduced_lunch    0.00000   0.58411  -0.09333   0.059
#> 
#> Bayesian Estimates Based on 10000 draws
#> 
#>                                   Post.Mode Post.Mean  Post.Sd MC Error
#> (Intercept)::(Intercept)           98.70159  98.72679  7.25423    0.073
#> (Intercept)::private_school         0.71693   0.71448  8.88178    0.089
#> (Intercept)::title1                22.93297  22.96431  9.31385    0.093
#> (Intercept)::free_reduced_lunch    -0.38190  -0.38312  0.14556    0.001
#> distracted_ppvt::(Intercept)       -1.60229  -1.57857  1.53076    0.015
#> distracted_a1::(Intercept)          0.78681   0.75771  3.63296    0.036
#> distracted_a1::free_reduced_lunch  -0.09218  -0.09231  0.05012    0.001
#>                                   Pr(Prior_tail)
#> (Intercept)::(Intercept)                  0.0756
#> (Intercept)::private_school               0.4674
#> (Intercept)::title1                       0.0063
#> (Intercept)::free_reduced_lunch           0.0041
#> distracted_ppvt::(Intercept)              0.1492
#> distracted_a1::(Intercept)                0.4153
#> distracted_a1::free_reduced_lunch         0.0334
#> 
#> Distribution Percentiles
#> 
#>                                         1.0%       2.5%       5.0%     Median
#> (Intercept)::(Intercept)           81.868168  84.734978  86.900613  98.728120
#> (Intercept)::private_school       -19.987908 -16.607816 -13.886773   0.724683
#> (Intercept)::title1                 1.359595   4.815026   7.853031  22.986348
#> (Intercept)::free_reduced_lunch    -0.719311  -0.663217  -0.618687  -0.383932
#> distracted_ppvt::(Intercept)       -5.126314  -4.562351  -4.098339  -1.581182
#> distracted_a1::(Intercept)         -7.749593  -6.274689  -5.195309   0.778453
#> distracted_a1::free_reduced_lunch  -0.208035  -0.189575  -0.173844  -0.092025
#>                                        95.0%      97.5%   99.0%
#> (Intercept)::(Intercept)          110.645543 113.112455 116.037
#> (Intercept)::private_school        15.235681  18.197290  21.400
#> (Intercept)::title1                38.350818  41.179280  44.366
#> (Intercept)::free_reduced_lunch    -0.142926  -0.094848  -0.044
#> distracted_ppvt::(Intercept)        0.904222   1.395001   2.009
#> distracted_a1::(Intercept)          6.749852   7.887948   9.190
#> distracted_a1::free_reduced_lunch  -0.010015   0.005335   0.022
#> 
#> === Block 1 group effects ===
#> 
#> Distribution of posterior modes across groups
#> 
#>                     Mean       SD      Min       Q1   Median       Q3      Max
#> (Intercept)      91.2640  14.0905  61.4926  82.4736  92.3806 102.7179 114.1792
#> distracted_ppvt  -1.6033   0.9026  -3.4764  -1.8115  -1.4829  -1.1302   0.7325
#> distracted_a1    -4.3086   3.7974 -11.3164  -7.0521  -5.7084  -2.0022   1.9539
#>                 MCMC.mean
#> (Intercept)        91.265
#> distracted_ppvt    -1.587
#> distracted_a1      -4.341
#> 
#> Per-group detail (requested levels)
#> 
#>  group (Intercept).mode (Intercept).mean (Intercept).sd distracted_ppvt.mode
#>     30            80.59            80.67          6.133               -2.773
#>     15            96.69            96.76          6.790               -1.717
#>  distracted_ppvt.mean distracted_ppvt.sd distracted_a1.mode distracted_a1.mean
#>                -2.774              2.814            -11.316            -11.367
#>                -1.709              2.825              1.396              1.352
#>  distracted_a1.sd
#>             3.283
#>             4.532

Variance convergence during the main phase

The exact iid route of Chapter-C02 produces no sweep history because it never runs a Markov chain. The Gibbs route does, and the history is stored on the fit as fit$sweep_history. The variance chart measures Claim 3 of vignette("Chapter-C03") empirically. That claim is an equality, $\Sigma_{11}^{-1/2}\Sigma^{(l)}_{11}\Sigma_{11}^{-1/2} = I - A^{2l}$, so after $l$ sweeps the cross-chain covariance of the Block 2 fixed effects should sit below the exact posterior covariance by exactly $A^{2l}$ on the whitened scale, and rise toward it as $l$ grows.

Because $\sigma^2$ and $\Psi$ are known here, that target covariance is known as well. plot_var_convergence() resolves it automatically when both the design and the measurement prior list are supplied — it calls lmerb_posterior_covariance() internally and uses the result as $\mathrm{Var}_{\mathrm{ref}}$, not the empirical variance at the final sweep. The artifact stores a stripped design and a measurement prior list built the same way the sampler does: Block~1 precision from pfamily_list() (the fixed $\tau^2_k$ plug-ins), not the REML group.Sigma row of Prior_Setup_GLMM() (which lists the reference lmer variances in its lmer column). Using the wrong $\Psi$ in the reference would make the slope traces sit far above 1 even when the chain has converged.

plot_var_convergence(
  C04$sweep_history,
  design                 = C04$design,
  measurement_prior_list = C04$measurement_prior_list,
  whitened               = FALSE,
  n_chains               = C04$meta$n_chains
)
#> 
#> === main -- Intercept predictors sweep history (Var / Var_exact ratio; base) ===

#> 
#> === main -- Slope predictors sweep history (Var / Var_exact ratio; base) ===

Cross-chain variance ratio by inner sweep, unwhitened coefficients

With whitened = FALSE, each trace is a named population effect and the y-axis is $\mathrm{Var}(l) / \mathrm{Var}_{\mathrm{ref}}$ for that coefficient, where $\mathrm{Var}(l)$ is the cross-chain variance at inner sweep $l$ and $\mathrm{Var}_{\mathrm{ref}}$ is the corresponding diagonal entry of the exact posterior covariance. Under Claim 3 the ratio is bounded above by 1 and rises toward it as $l$ grows, so a well-calibrated run shows every trace flattening just below 1 by sweep 25. The dotted band, when n_chains is supplied, is a naive Monte Carlo ruler for how much departure from 1 is sampling noise in the numerator alone — appropriate here because the denominator is exact.

The chart uses the main-phase history only. This route has no pilot stage; every chain starts at the exact posterior mean and the sweeps carry the cross-chain spread up to the target variance.

References

Nygren, Kjell N. 2020. “On the Total Variation Distance Between Multivariate Normal Densities with Applications to Two-Block Gibbs Samplers.” Unpublished manuscript.

Chapter C04: The two-block Gibbs sampler