C02 <- readRDS(system.file(
"extdata", "Chapter-C02-bwc-iid.rds", package = "lmebayesCore"
))
meta <- C02$metaThis chapter fits the model of Chapter-C01 under the
simplest prior specification the package offers, and shows that the
resulting posterior can be sampled exactly rather than iteratively.
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:
C02$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.80365The 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, so the intercept
block, the distracted_ppvt block and the
distracted_a1 block each get their own mean vector and
covariance matrix:
round(C02$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.584The pattern in the mean column is worth noting. Every coefficient is centered at zero except the grand intercept, so the prior mean vector describes a null, intercept-only model: no covariate has any effect, and no group effect varies with anything at level two. The prior therefore carries the same interpretation the null model does in classical statistics, and the analysis can be read as asking how far the data move each coefficient away from that null. What differs is that here the null is a starting point that the data revise rather than a hypothesis that is accepted or rejected.
This is the same default the glmbayes family of packages uses for ordinary iid regression: a multivariate normal prior on the regression coefficients, centered at the null model and calibrated from the data rather than specified by hand. What the mixed model adds is that the specification is repeated once per group effect.
The setup is built by two calls. Prior_Setup_GLMM()
performs the calibration, and pfamily_list() wraps the
result as the Block 2 prior families the sampler consumes:
pop.pwt = 0.01 gives the prior on \(\gamma\) the weight of about one percent of
the data. Because every component of pf comes back as
dNormal, the variance components are known, which is what
the first table above reports.
The sampler takes the design and the two prior objects directly:
fit <- rlmerb(
n = 1000,
design = model_setup(form, data = dat),
pfamily_list = pf,
dispprior_list = list(dispersion = ps$group.dispersion),
sim_method = "DEFAULT"
)The design argument carries the data and both design
matrices, pfamily_list supplies the Block 2 priors from the
previous section, and dispprior_list supplies the fixed
\(\sigma^2\). Passing a bare
list(dispersion = ...) rather than a dGamma()
prior is what declares that dispersion known.
Under this specification the joint posterior over the group effects
and the population effects is exactly multivariate normal, which admits
a direct sampling scheme with no Markov chain.
sim_method = "DEFAULT" selects it.
The same model is fitted again later by a two-block Gibbs sampler. Running both against one posterior is deliberate: the answer here is exact, so it serves as the standard against which the iterative sampler is judged.
Collect the group effects as \(\beta = (\beta_1', \dots, \beta_J')'\). Up to an additive constant, the log posterior splits into the three pieces of the model:
\[ -2\log \pi(\gamma, \beta \mid y) = \underbrace{\frac{1}{\sigma^2}\sum_{j} \left\|y_j - D_j\beta_j\right\|^2}_{\text{(i) log-likelihood}} + \underbrace{\sum_{j} (\beta_j - \mathcal{W}_j\gamma)' \Psi^{-1} (\beta_j - \mathcal{W}_j\gamma)}_{\text{(ii) prior on the group effects}} + \underbrace{(\gamma - \mu_0)' V^{-1} (\gamma - \mu_0)}_{\text{(iii) prior on the population effects}}. \]
Because \(\sigma^2\), \(\Psi\) and \(V\) are all known, every term is quadratic in \((\gamma, \beta)\), and a quadratic log density is a normal density. The posterior is therefore multivariate normal, and it can be factorized into a marginal for \(\gamma\) and a conditional for \(\beta\).
Integrating the group effects out of (i) and (ii) is a linear-Gaussian calculation. Substituting \(\beta_j = \mathcal{W}_j\gamma + u_j\) with \(u_j \sim N(0, \Psi)\) gives
\[ y_j \mid \gamma \sim N\!\left(D_j\mathcal{W}_j\,\gamma,\; V_j\right), \qquad V_j = \sigma^2 I_{n_j} + D_j \Psi D_j', \]
so each group contributes a Gaussian likelihood in \(\gamma\) with design \(X_j = D_j\mathcal{W}_j\). Combining these with (iii) gives the marginal posterior
\[ \gamma \mid y \sim N\!\left(m_\gamma,\; P_\gamma^{-1}\right), \qquad \begin{aligned} P_\gamma &= V^{-1} + \sum_j X_j' V_j^{-1} X_j, \\ m_\gamma &= P_\gamma^{-1}\Big(V^{-1}\mu_0 + \sum_j X_j' V_j^{-1} y_j\Big). \end{aligned} \]
The term \(D_j \Psi D_j'\) in \(V_j\) is what distinguishes this from ordinary regression: it inflates the variance of the groups that would otherwise be most influential, which is how the model discounts a group that is unusual rather than merely large.
Given \(\gamma\), terms (i) and (ii) separate across groups, so the group effects are conditionally independent — each is an ordinary conjugate normal update, combining that group’s own data with the population prediction \(\mathcal{W}_j\gamma\):
\[ \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} \]
The precision weights in \(b_j\) are the shrinkage: a group with plenty of data has a large \(D_j' D_j\) and stays near its own fit, while a small group is pulled toward \(\mathcal{W}_j\gamma\).
The two results above factorize the posterior exactly,
\[ \pi(\gamma, \beta \mid y) = \pi(\gamma \mid y)\;\prod_j \pi(\beta_j \mid \gamma, y), \]
so one draw from the joint posterior is obtained in two steps: draw \(\gamma\) from its marginal, then draw each \(\beta_j\) from its conditional given that value. Repeating the pair produces independent draws from the target itself — no burn-in to discard, no autocorrelation, and no convergence rate to bound. It also makes the ICM point estimate the exact posterior mean, since the mode and mean of a normal distribution coincide.
summary() collects the draws into tables:
The output is organized by block, following the structure of the model. After the call and a header giving the number of draws, the engine used and the data dimensions, it reports in turn:
lmer
variance components, for orientation.lmer reference, then the posterior
mode, mean, standard deviation and Monte Carlo error, then distribution
percentiles. Rows are named component::parameter, following
the indexing described in Chapter-C01.groups.The groups argument is optional; without it the
per-group detail is replaced by a note pointing at
fit$groupef.mode. The two levels requested here are the
largest and smallest schools in the data.
Two features of the population effect tables follow from the default
pop.pwt = 0.01. Because the prior is given the weight of
only about one percent of the data, it barely competes with the
likelihood: the posterior means land close to the lmer
estimates, displaced slightly toward the prior means of the previous
section. Since those prior means are the null model, that displacement
is a mild shrinkage toward zero for every coefficient but the
intercept.
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, and most within a twentieth, so no
coefficient in the table supports a different reading under one fit than
the other. Exact agreement was never expected: the prior carries a small
amount of weight, and this chapter fixes the variance components where
lmer estimates them, so the two are answering slightly
different questions. The differences that remain are too small to be
worth interpreting individually.
The Pr(Prior_tail) column has a related reading. For
every row it measures the same thing: the posterior mass falling on the
opposite side of the prior mean from the posterior estimate. It is thus
a one-sided statement about how far the data have moved a coefficient
away from its prior value, and how much posterior belief remains on the
other side of it.
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 answers the same question a one-sided p-value does. Across them
it stays in the same range as half the two-sided p-value the
corresponding lmer fit would report, typically running
somewhat smaller because the posterior standard deviations are not
identical to the lmer standard errors, and the ordering it
induces matches the classical one. The intercept is read the same way,
but against the calibrated prior mean rather than the origin, so it is
not comparable to a test against zero.
The output below was produced by
data-raw/make_Chapter-C02.R with n = 10000 and
set.seed(1), then stored, so that building this vignette
does not rerun the sampler.
C02$summary
#> Call
#> rlmerb(n = N_DRAWS, design = design, pfamily_list = pf, dispprior_list = list(dispersion = ps$group.dispersion),
#> progbar = FALSE, verbose = FALSE, print_icm_table = FALSE,
#> sim_method = "DEFAULT")
#>
#> Bayesian mixed model fit based on 10000 draws [exact iid draws]
#> 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.65977 7.24785 0.072
#> (Intercept)::private_school 0.71693 0.75439 8.96814 0.090
#> (Intercept)::title1 22.93297 22.95896 9.24689 0.092
#> (Intercept)::free_reduced_lunch -0.38190 -0.38139 0.14460 0.001
#> distracted_ppvt::(Intercept) -1.60229 -1.59775 1.52220 0.015
#> distracted_a1::(Intercept) 0.78681 0.74485 3.73106 0.037
#> distracted_a1::free_reduced_lunch -0.09218 -0.09169 0.05139 0.001
#> Pr(Prior_tail)
#> (Intercept)::(Intercept) 0.0771
#> (Intercept)::private_school 0.4631
#> (Intercept)::title1 0.0070
#> (Intercept)::free_reduced_lunch 0.0039
#> distracted_ppvt::(Intercept) 0.1489
#> distracted_a1::(Intercept) 0.4199
#> distracted_a1::free_reduced_lunch 0.0372
#>
#> Distribution Percentiles
#>
#> 1.0% 2.5% 5.0% Median
#> (Intercept)::(Intercept) 82.008226 84.527706 86.791378 98.657050
#> (Intercept)::private_school -20.316336 -17.281343 -14.340589 0.835984
#> (Intercept)::title1 1.699511 4.935409 8.031938 22.942436
#> (Intercept)::free_reduced_lunch -0.715026 -0.669400 -0.621988 -0.378938
#> distracted_ppvt::(Intercept) -5.133680 -4.571532 -4.109684 -1.598932
#> distracted_a1::(Intercept) -8.022068 -6.525714 -5.410193 0.763723
#> distracted_a1::free_reduced_lunch -0.209467 -0.190931 -0.176635 -0.091491
#> 95.0% 97.5% 99.0%
#> (Intercept)::(Intercept) 110.662854 112.981720 115.985
#> (Intercept)::private_school 15.372904 18.123172 20.765
#> (Intercept)::title1 37.948480 41.225127 44.715
#> (Intercept)::free_reduced_lunch -0.144732 -0.099333 -0.048
#> distracted_ppvt::(Intercept) 0.888930 1.389696 1.900
#> distracted_a1::(Intercept) 6.827513 8.043016 9.468
#> distracted_a1::free_reduced_lunch -0.006376 0.008700 0.026
#>
#> === 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.275
#> distracted_ppvt -1.609
#> distracted_a1 -4.313
#>
#> Per-group detail (requested levels)
#>
#> group (Intercept).mode (Intercept).mean (Intercept).sd distracted_ppvt.mode
#> 30 80.59 80.67 6.175 -2.773
#> 15 96.69 96.76 6.693 -1.717
#> distracted_ppvt.mean distracted_ppvt.sd distracted_a1.mode distracted_a1.mean
#> -2.782 2.826 -11.316 -11.383
#> -1.719 2.824 1.396 1.349
#> distracted_a1.sd
#> 3.312
#> 4.556lmer fitThe lmer columns used throughout this chapter come from
a single REML fit that model_setup() performs on all
observations. Prior_Setup_GLMM() returns it as
mer_fit, and the same object is stored on the design, so
either of these reaches it:
It serves two purposes: it calibrates the default priors of the first section, and it populates the reference columns in the summary tables. Its own summary is:
cat(C02$lmer_summary, sep = "\n")
#> Linear mixed model fit by REML ['lmerMod']
#> Formula:
#> score_ppvt ~ private_school + title1 + free_reduced_lunch + distracted_ppvt +
#> distracted_a1 + free_reduced_lunch:distracted_a1 + ((1 |
#> school_id) + (0 + distracted_ppvt | school_id) + (0 + distracted_a1 |
#> school_id))
#> Data: bwc
#> Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
#>
#> REML criterion at convergence: 3565.2
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -2.59677 -0.68699 -0.05662 0.70471 2.90372
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> school_id (Intercept) 191.78 13.849
#> school_id.1 distracted_ppvt 26.03 5.102
#> school_id.2 distracted_a1 41.05 6.407
#> Residual 366.41 19.142
#> Number of obs: 400, groups: school_id, 34
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) 99.16637 7.88455 12.577
#> private_school -0.16629 9.73406 -0.017
#> title1 23.24038 10.07828 2.306
#> free_reduced_lunch -0.38866 0.15784 -2.462
#> distracted_ppvt -1.61210 1.75946 -0.916
#> distracted_a1 1.46278 4.16277 0.351
#> free_reduced_lunch:distracted_a1 -0.09333 0.05870 -1.590
#>
#> Correlation of Fixed Effects:
#> (Intr) prvt_s title1 fr_rd_ dstrc_ dstr_1
#> privat_schl -0.800
#> title1 0.093 -0.076
#> fr_rdcd_lnc -0.682 0.547 -0.747
#> dstrctd_ppv -0.051 0.009 0.013 -0.002
#> distractd_1 -0.094 -0.003 0.000 0.062 -0.114
#> fr_rdcd_:_1 0.091 -0.003 -0.007 -0.077 0.017 -0.859Two points connect this back to the chapter. First, the fixed effects
are reported in lmer’s flat parameterization, with the
cross-level term written as
free_reduced_lunch:distracted_a1. The summary tables carry
the same seven coefficients but regroup them as
component::parameter, so that interaction appears as
distracted_a1::free_reduced_lunch — a level-two predictor
of the distracted_a1 slope. Chapter-C01 covers
the mapping.
Second, the Random effects block is the quantity this
chapter holds fixed, and it is worth seeing that it is not fixed at
these values:
round(C02$summary$tau2_prior_overview[, c("E[tau2]", "lmer")], 2)
#> E[tau2] lmer
#> (Intercept) 159.14 191.79
#> distracted_ppvt 6.74 26.03
#> distracted_a1 11.48 41.05
#> Residual 388.80 366.41The residual variance is calibrated close to the REML value, and the
intercept variance is moderately tighter, but the two slope components
are held near a quarter of what lmer estimates. This is the
one place where the two fits are far apart, and it is a difference in
specification rather than in estimation: a smaller \(\Psi\) shrinks the group effects harder
toward their population predictions, so the group effects reported in
the previous section are more tightly pooled than lmer
would make them.
That the population effects nonetheless agree so closely is worth
noting. Holding \(\Psi\) at these
values costs little there, which is what makes this specification a
reasonable starting point even though it is not the one
lmer fits. The later chapters estimate the variance
components rather than fixing them, and are the place to look for what
changes when the Random effects block above is treated as
unknown.
vignette("Chapter-C01") — the data, model and design
used here.vignette("Chapter-B04") — route selection and the
compiled kernel behind the exact engine.?rlmerb, ?summary.rlmerb,
?model_setup, ?Prior_Setup_GLMM,
?pfamily_listdata-raw/make_Chapter-C02.R in the package sources for
the script that produced the stored results.