Every fit in this package goes through the same three steps.
formula + data
|
v
model_setup() what model am I actually fitting?
|
v
Prior_Setup_GLMM() calibrate priors against a reference lmer fit
|
v
pfamily_list() one prior object per random-effect coefficient
|
v
lmerb() / glmerb()You can skip the first step — lmerb() calls it
internally — but you should not, because model_setup() is
where you find out whether the model you wrote is the model you
meant.
model_setup()design <- model_setup(form, data = dat)
design
#> Call:
#> model_setup(formula = form, data = dat)
#>
#> --- Level 1 (Within-Group) ---
#> Reaction ~ (Intercept) + Days_c
#>
#> Observations : 180
#> Group-effect coefficients: 2
#> Group : Subject [18 levels]
#> Full-rank D_j: 18 of 18 groups
#> Full-rank & estimable: 18 of 18 full-rank (18 of 18 total Subject)
#>
#> --- Level 2 (Across-Group) ---
#> (Intercept) ~ 1
#> Days_c ~ 1
#>
#> --- Level 2 Rank ---
#> (Restricted to 18 estimable Subject)
#>
#> (Intercept) groups=18 predictors=1 full-rank
#> Days_c groups=18 predictors=1 full-rank
#>
#> Rank-deficient flags:
#> (Intercept) FALSE
#> Days_c FALSE
#>
#> popef.rank_ok: TRUE -- model rank looks OKThe printed report has three parts.
Level 1 lists the within-group design: which coefficients vary by group, how many observations there are, and how many groups.
Level 2 lists, for each random coefficient, which
group-level predictors enter its prior mean. A bare ~ 1
means “shrink toward a single common value”. Anything more means “shrink
toward a fitted line”.
Rank reports two separate checks.
Full-rank D_j: 18 of 18 groups is the within-group check:
does each group have enough distinct covariate values to identify its
own coefficients? Groups that fail stay in the fit but are set aside for
the next check. The Level 2 Rank block is the across-group check, run on
the groups that passed: do they span the level-2 predictor space?
popef.rank_ok summarizes it.
The asymmetry matters. A rank-deficient group is survivable — it
simply learns nothing from its data in that direction. A failed
level-2 check is not: the population parameter \(\gamma\) is then unidentified by every
data-bearing group, and the chain is null recurrent in that direction
even with a proper prior. Both conditions are derived in
vignette("Chapter-B08", package = "lmebayesCore").
The object itself carries the matrices:
names(design)
#> [1] "group_name" "group" "groupef.names"
#> [4] "y" "weights" "offset"
#> [7] "D" "W" "popef.moderation"
#> [10] "call" "formula" "family"
#> [13] "terms" "data" "dispformula"
#> [16] "lmer" "varcorr" "Psi"
#> [19] "dispersion" "groupef.rank" "groupef.estimable"
#> [22] "groupef.glm_check" "popef.rank" "popef.deficient"
#> [25] "popef.rank_ok"
dim(design$D) # stacked within-group design
#> [1] 180 2
str(design$groupef.names) # names of the random coefficients
#> chr [1:2] "(Intercept)" "Days_c"Prior_Setup_GLMM()ps <- Prior_Setup_GLMM(form, data = dat, pop.pwt = 0.01)
ps
#> Call: Prior_Setup_GLMM()
#>
#> --- Setup ---
#> Model family, prior-mean sources, and reference fit used for calibration.
#>
#> family : gaussian (identity link)
#> pop.intercept_source : null_model
#> pop.effects_source : null_effects
#> dispformula : ~1
#> calibration_source : lme4 (fixef / tau^2_k / sd_tau)
#>
#> --- Group dispersion (sigma^2) ---
#> Pooled Gamma prior on sigma^2 is stored in group.ing_prior.
#> (dGamma_list() requires a per-group dispformula.)
#> group.dispersion : 709.8 [pooled S_marg (A12 3.3.4)]
#> classical lmer sigma^2 : 654.9 (diagnostic)
#>
#> --- Design check (from model_setup) ---
#> Identifiability flags copied from the design object (not re-checked here).
#>
#> Full-rank groups (algebraic D_j): 18 of 18 Subject
#> Full-rank & estimable : 18 of 18 full-rank (18 of 18 total Subject)
#>
#> --- Group-effect covariance (group.Sigma) ---
#> Diagonal Psi with entries tau^2_k for each group-effect coefficient.
#>
#> (Intercept) Days_c
#> (Intercept) 1408.73 0.0000
#> Days_c 0.00 35.0717
#>
#> --- Population priors (pop.prior_list) ---
#> One prior list per group-effect coefficient (input to pfamily_list()).
#>
#> [(Intercept)]
#> pop.pwt : 0.01
#> pop.dispersion.pwt : 0.01 [derived from pwt]
#> pop.dispersion.nprior : 0.1818 (= J * pwt_disp / (1 - pwt_disp))
#>
#> mu:
#> (Intercept)
#> 298.5079
#>
#> Sigma:
#> (Intercept)
#> (Intercept) 8108.233
#>
#> dispersion : 1346.1456 (tau^2_k; A12 3.3.4 marginal)
#> pop.dispersion.ref : 1408.7301 (classical VarCorr / Psi)
#> ING tau^2 window : [632.3, 2998] (0.01/0.99 limiting-posterior quantiles; upper/tau2 = 2.23)
#> ING shape, rate : 1.091, 122.4 (Gamma on 1/tau^2_k; for ptypes = "dIndependent_Normal_Gamma")
#>
#> [Days_c]
#> pop.pwt : 0.01
#> pop.dispersion.pwt : 0.01 [derived from pwt]
#> pop.dispersion.nprior : 0.1818 (= J * pwt_disp / (1 - pwt_disp))
#>
#> mu:
#> (Intercept)
#> 0
#>
#> Sigma:
#> (Intercept)
#> (Intercept) 236.5569
#>
#> dispersion : 29.3722 (tau^2_k; A12 3.3.4 marginal)
#> pop.dispersion.ref : 35.0717 (classical VarCorr / Psi)
#> ING tau^2 window : [13.8, 65.42] (0.01/0.99 limiting-posterior quantiles; upper/tau2 = 2.23)
#> ING shape, rate : 1.091, 2.67 (Gamma on 1/tau^2_k; for ptypes = "dIndependent_Normal_Gamma")This is the calibration step. It fits a reference lmer()
model and uses its estimates to place priors on a scale that is
commensurate with the data, so you specify how much prior
information you want rather than guessing at numeric
hyperparameters.
The central argument is pop.pwt, the population
prior weight. Setting pop.pwt = 0.01 means: give
the population-level prior about as much weight as one percent of the
effective sample. The prior covariance is scaled accordingly.
| Setting | Effect |
|---|---|
pop.pwt = 0.01 (default-ish) |
Weak prior; posterior close to lmer() |
pop.pwt = 0.05 |
Noticeably regularizing |
pop.pwt = 0.25 |
Strong shrinkage toward the prior mean |
The equivalent alternative is pop.nprior, which states
the prior weight as a number of observations rather than a fraction, or
pop.sd, which sets prior standard deviations directly.
Two arguments control the prior mean.
pop.intercept_source — "null_model"
(default) centers the intercept at the intercept-only fit, so the prior
says “I know roughly the overall level” without pre-judging any
effect.pop.effects_source — "null_effects"
(default) centers every non-intercept coefficient at
zero. This is the conservative choice: the prior
expresses no opinion about direction, and any signal in the posterior
comes from the data.Setting either to "full_model" centers the prior on the
reference lmer estimates. That is legitimate for
sensitivity analysis but is not a prior in the usual sense — it uses the
data twice. Use pop.mu to supply your own means when you
have real external information.
names(ps)
#> [1] "formula" "family"
#> [3] "pop.pwt" "pop.nprior"
#> [5] "pop.sd" "pop.mu"
#> [7] "group.dispersion" "pop.dispersion.pwt"
#> [9] "pop.dispersion.nprior" "group.dispersion.pwt"
#> [11] "group.dispersion.nprior" "group.alpha_target"
#> [13] "group.pwt_calibration" "group.max_disp_perc"
#> [15] "pop.max_disp_perc" "dispformula"
#> [17] "pop.intercept_source" "pop.effects_source"
#> [19] "data" "block_formula"
#> [21] "group.tau_sd" "design"
#> [23] "mer_fit" "fit_ref"
#> [25] "group.dispersion.fit" "group.dispersion.ref"
#> [27] "pop.dispersion.ref" "calibration_source"
#> [29] "group.Sigma" "pop.prior_list"
#> [31] "pop.ing_prior" "group.ing_prior"The names are prefixed by which part of the model they configure.
| Prefix | Configures |
|---|---|
pop.* |
Block 2: the population coefficients \(\gamma\) and the variance components \(\tau^2_k\) |
group.* |
The observation-level dispersion \(\sigma^2\) |
The two you pass on to the sampler are:
ps$group.dispersion # sigma^2: observation variance
#> [1] 709.826
#> attr(,"source")
#> [1] "pooled S_marg (A12 3.3.4)"
#> attr(,"classical")
#> [1] 654.941
str(ps$pop.prior_list, max.level = 2)
#> List of 2
#> $ (Intercept):List of 3
#> ..$ mu : Named num 299
#> .. ..- attr(*, "names")= chr "(Intercept)"
#> ..$ Sigma : num [1, 1] 8108
#> .. ..- attr(*, "dimnames")=List of 2
#> ..$ dispersion: num 1346
#> $ Days_c :List of 3
#> ..$ mu : Named num 0
#> .. ..- attr(*, "names")= chr "(Intercept)"
#> ..$ Sigma : num [1, 1] 237
#> .. ..- attr(*, "dimnames")=List of 2
#> ..$ dispersion: num 29.4ps$group.dispersion is the fixed observation variance,
taken from the reference fit. ps$pop.prior_list holds, for
each random coefficient, the prior mean vector mu, prior
covariance Sigma, and the dispersion tau^2_k
that governs how tightly group coefficients cluster around \(\mathcal{W}_j\gamma\).
When you ask for sampled variance components (Chapter 07),
Prior_Setup_GLMM() also computes a truncation window:
ps_ing <- Prior_Setup_GLMM(form, data = dat, pop.pwt = 0.01,
pop.dispersion.pwt = 0.2)
ps_ing$pop.ing_prior
#> $`(Intercept)`
#> $`(Intercept)`$shape
#> [1] 3.25
#>
#> $`(Intercept)`$rate
#> [1] 3028.828
#>
#> $`(Intercept)`$disp_lower
#> [1] 632.3273
#>
#> $`(Intercept)`$disp_upper
#> [1] 2998.203
#>
#> $`(Intercept)`$max_disp_perc
#> [1] 0.99
#>
#>
#> $Days_c
#> $Days_c$shape
#> [1] 3.25
#>
#> $Days_c$rate
#> [1] 66.08751
#>
#> $Days_c$disp_lower
#> [1] 13.79707
#>
#> $Days_c$disp_upper
#> [1] 65.41931
#>
#> $Days_c$max_disp_perc
#> [1] 0.99disp_lower and disp_upper bound the support
of each \(\tau^2_k\). They are
quantiles of the limiting posterior for the variance component,
not of the prior — the reasoning is in
inst/ING_TRUNCATION_WINDOW.md in the package sources. The
lower bound matters for more than propriety: the sampler evaluates its
convergence rate at disp_lower, because smaller \(\tau^2\) means tighter coupling between
blocks and therefore a slower chain. A bound computed there is valid
across the whole window. See
vignette("Chapter-B03", package = "lmebayesCore").
pfamily_list()pf <- pfamily_list(ps)
pf
#>
#> [[(Intercept)]]
#>
#> Call: dNormal(mu = mu_k, Sigma = Sig_k, dispersion = d_k)
#>
#> Prior Family: dNormal
#>
#> Prior List:
#>
#> $mu
#> [,1]
#> (Intercept) 298.5079
#>
#> $Sigma
#> (Intercept)
#> (Intercept) 8108.233
#>
#> $dispersion
#> [1] 1346.146
#>
#> $ddef
#> [1] FALSE
#>
#> attr(,"Prior Type")
#> [1] "dNormal"
#>
#> [[Days_c]]
#>
#> Call: dNormal(mu = mu_k, Sigma = Sig_k, dispersion = d_k)
#>
#> Prior Family: dNormal
#>
#> Prior List:
#>
#> $mu
#> [,1]
#> (Intercept) 0
#>
#> $Sigma
#> (Intercept)
#> (Intercept) 236.5569
#>
#> $dispersion
#> [1] 29.37223
#>
#> $ddef
#> [1] FALSE
#>
#> attr(,"Prior Type")
#> [1] "dNormal"One pfamily object per random-effect coefficient. The
prior family determines what happens to that coefficient’s variance
component.
| Family | \(\tau^2_k\) | Chapter |
|---|---|---|
dNormal |
Fixed at the reference estimate | 01, 04 |
dIndependent_Normal_Gamma |
Sampled, with a Gamma prior on the precision | 07 |
Ask for the second with ptypes:
pf_ing <- pfamily_list(ps_ing, ptypes = "dIndependent_Normal_Gamma")
class(pf_ing[["Days_c"]])
#> [1] "pfamily"ptypes accepts one value applied to all coefficients, or
a named vector if you want to sample some variance components and fix
others.
Nothing forces you to use the calibrated object. A
pfamily_list is a plain named list, so you can replace an
element:
pf_manual <- pf
pf_manual[["Days_c"]] <- dNormal(
mu = matrix(5, 1, 1), # prior mean slope: 5 ms per day
Sigma = matrix(4, 1, 1), # prior SD 2 on that mean
dispersion = 25 # tau^2: subject slopes vary with SD 5
)
pf_manual[["Days_c"]]$prior_list[c("mu", "Sigma", "dispersion")]
#> $mu
#> [,1]
#> [1,] 5
#>
#> $Sigma
#> [,1]
#> [1,] 4
#>
#> $dispersion
#> [1] 25Three numbers, three different jobs, and they are easy to confuse:
mu — the prior mean of the
population slope \(\gamma_k\).Sigma — how uncertain you are about
that population slope.dispersion — \(\tau^2_k\), how much individual
subjects scatter around it. This is the variance component, not
prior uncertainty.Setting Sigma small says “I know the average slope”.
Setting dispersion small says “everyone’s slope is nearly
the same”. They are independent statements.
set.seed(1)
fit <- lmerb(
form,
data = dat,
pfamily_list = pf,
dispersion_ranef = ps$group.dispersion,
n = 200L,
progbar = FALSE
)fit
#> Call:
#> lmerb(formula = form, data = dat, pfamily_list = pf, dispersion_ranef = ps$group.dispersion,
#> n = 200L, progbar = FALSE)
#>
#> Bayesian linear mixed model [200 draws, exact iid]
#> Formula: Reaction ~ Days_c + (1 + Days_c || Subject)
#>
#> Random effects (variance components fixed at lmer estimates):
#> Groups Name Std.Dev.
#> Subject (Intercept) 37.533
#> Subject.1 Days_c 5.922
#> Residual 25.592
#> Number of obs: 180, groups: Subject, 18
#>
#> --- Posterior means (ICM exact, under fixed variance components) ---
#>
#> RE component parameter fixef.mode fixef.means draws SD
#> ------------ ----------- ------------ ------------ ----------
#> (Intercept) (Intercept) 298.5079 297.8895 9.4520
#> Days_c (Intercept) 10.3748 10.2699 1.3633Because the whole prior is one argument, sensitivity analysis is a
loop over pop.pwt:
pwts <- c(0.01, 0.10, 0.50)
res <- lapply(pwts, function(w) {
ps_w <- Prior_Setup_GLMM(form, data = dat, pop.pwt = w)
f <- lmerb(form, data = dat,
pfamily_list = pfamily_list(ps_w),
dispersion_ranef = ps_w$group.dispersion,
n = 200L, progbar = FALSE)
unlist(f$popef.means)
})tab <- do.call(rbind, res)
rownames(tab) <- paste0("pop.pwt = ", pwts)
round(tab, 3)
#> (Intercept).(Intercept) Days_c.(Intercept)
#> pop.pwt = 0.01 296.872 10.451
#> pop.pwt = 0.1 299.947 9.133
#> pop.pwt = 0.5 298.466 3.609The intercept barely moves — it is centered at the null-model
intercept, so prior and data already agree. The Days_c
slope is centered at zero, so heavier prior weight pulls it toward zero.
That contrast is the whole point of
pop.effects_source = "null_effects": the prior is
deliberately skeptical about effects, and you can see exactly how much
skepticism costs.
~ 1.?Prior_Setup_GLMM, ?pfamily_list,
?model_setupvignette("Chapter-A12", package = "glmbayes") — the
derivations behind the Gaussian prior calibration this inherits.