Package 'lmebayesCore'

lmebayesCore: Core C++ Sampling Engine for lmebayes

Description

Core C++ engine for envelope-based iid samplers, two-block Gibbs mixed-model engines, and optional OpenCL acceleration. Full-featured developer backend for lmebayes and related extensions. End users should install lmebayes for lmer/glmer-style mixed-model workflows.

Details

Low-level entry points include envelope construction (EnvelopeBuild, EnvelopeOrchestrator), registered simulation pipelines (rNormalGLM_std, rIndepNormalGammaReg_std), matrix-input samplers (rglmb, rlmb), the two-block Gibbs mixed-model engines (rlmerb, rglmerb), and OpenCL kernel loaders. Formula interfaces lmerb and glmerb live in lmebayes; lmebayesCore supplies the sampling engine.

Vignette index (the B series)

The Chapter-B* vignettes document the engine itself: the theory it implements, the routes it selects between, and the diagnostics it emits. They assume the reader has met the models already, for which see the Chapter-00 to Chapter-13 curriculum in lmebayes.

vignette("Chapter-B00")

Roadmap, notation, and the correspondence to lme4.

vignette("Chapter-B01")

The two-block architecture: replicate chains, the sweep loop, and the call stack.

vignette("Chapter-B02")

Ergodicity and total-variation convergence; the rate $\lambda^*$ and the rank conditions.

vignette("Chapter-B03")

Calibrating m_convergence from tv_tol, and the pilot-versus-main cost trade-off.

vignette("Chapter-B04")

The known-vcov Gaussian route and the exact iid sampler.

vignette("Chapter-B05")

Estimated variance components (dIndependent_Normal_Gamma Block 2).

vignette("Chapter-B06")

The extended local rate diagnostic.

vignette("Chapter-B07")

Per-group measurement dispersion.

vignette("Chapter-B08")

Design, identifiability and rank conditions.

vignette("Chapter-B09")

Row-block engines.

vignette("Chapter-B10")

Sweep history and convergence plots.

Vignette index (the C series)

The Chapter-C* vignettes take one model class each and work through its posterior algebra, the eigenvalue calculation behind the convergence rate $\lambda^*$, and simulation results with convergence plots. Where the B series describes what the engine does, the C series derives why it is correct and then shows the numbers.

vignette("Chapter-C01")

The Gaussian running example: data, model, design, rank restriction and calibrated priors, shared by every Gaussian chapter in the series.

vignette("Chapter-C02")

The exact Gaussian iid sampler: known $\Psi$ and $\sigma^2$, closed-form posterior, no Markov chain.

vignette("Chapter-C03")

Total variation distances between multivariate normal densities: Nygren (2020) reproduced in full — the generalized error functions, Lemmas 1–2, Theorems 1–3, Claims 1–4, Corollary 1 and their proofs. Model-free.

vignette("Chapter-C04")

The same Gaussian model via the two-block Gibbs sampler: the Chapter-C03 rate and sweep budget instantiated on the running design, sweep history, and variance-convergence diagnostics.

The Chapter-A* vignettes shipped in the sources are inherited iid GLM reference material from glmbayes and are not built into the package tarball.

OpenCL startup checks

In interactive sessions, attaching the package with library(lmebayesCore) may emit a short packageStartupMessage when glmbayesCore_has_opencl() is FALSE but a GPU or OpenCL stack appears available on the host. Set options(glmbayes.quiet_opencl_startup = TRUE) to suppress attach notes (recommended for CI and R CMD check).

Author(s)

Maintainer: Kjell Nygren [email protected]

Other contributors:

• The R Core Team (R Mathlib sources, R stats modeling code, and derived/adapted routines) [contributor, copyright holder]

• The R Foundation (Portions of R Mathlib and R source code) [copyright holder]

• Ross Ihaka (R Mathlib and original R modeling infrastructure) [contributor, copyright holder]

• Robert Gentleman (Portions of R Mathlib) [contributor, copyright holder]

• Simon Davies (Original R glm implementation) [contributor]

• Morten Welinder (Portions of R Mathlib) [contributor, copyright holder]

• Martin Maechler (Portions of R Mathlib) [contributor]

References

There are no references for Rd macro ⁠\insertAllCites⁠ on this help page.

See Also

lmebayes for the end-user modelling package.

---

Test whether group random effects lie in widetilde B(delta_2).

Description

Test whether group random effects lie in widetilde B(delta_2).

Usage

beta_in_marginal_safe_set(b, beta_set)

Arguments

b	Group effects as a J x p_re matrix (rows = groups) or a stacked numeric vector in group_precision_floor order.
beta_set	A beta_marginal_safe_set object.

Value

Logical scalar: TRUE if $\Xi(\beta) \le r_{\mathrm{Gauss}}$.

See Also

rNormal_reg_group_safe, rNormalGLM_reg_group_safe

---

Marginal-mode beta posterior mode (gamma integrated out)

Description

Finds the Newton mode $\beta^\dagger$ of the gamma-integrated marginal posterior on stacked group coefficients. Reusable input for beta_marginal_safe_set and group_precision_floor.

Usage

beta_marginal_mode(
  design,
  pfamily_list,
  family,
  dispprior_list = NULL,
  offset = NULL,
  weights = NULL,
  tol = 1e-10,
  maxit = 200L
)

Arguments

design	A model_setup list.
pfamily_list	Block~2 prior list from pfamily_list().
family	A family object.
dispprior_list	Optional Block~1 dispersion prior for gaussian().
offset, weights	Observation offset and weights; default from design.
tol, maxit	Newton controls (unused beyond initial solve).

Value

An object of class "beta_marginal_mode" with beta, hessian, f_mode, engine, prior_int, group_data, P_b, and family.

See Also

beta_marginal_safe_set, group_precision_floor

---

Marginal-mode beta safe set for Rosenthal / TV certificates

Description

Build the beta-safe set $\widetilde B(\delta_2) = \{\Xi(\beta) \le r_{\mathrm{Gauss}}(n,\delta_2)\}$ at fixed $\beta^\dagger$ from beta_marginal_mode. Does not compute data-precision floors; use group_precision_floor for $\Gamma_j^{\mathrm{LB}}$.

Usage

beta_marginal_safe_set(
  design = NULL,
  pfamily_list = NULL,
  family = NULL,
  delta_2 = 0.01,
  beta_mode = NULL,
  dispprior_list = NULL,
  offset = NULL,
  weights = NULL,
  tol = 1e-10,
  maxit = 200L,
  verbose = FALSE
)

Arguments

design	A model_setup list (used if beta_mode is NULL).
pfamily_list	Block~2 prior list (used if beta_mode is NULL).
family	A family object (used if beta_mode is NULL).
delta_2	Tail mass budget in $(0,1)$ for the Laplace reference at $r_{\mathrm{Gauss}}$.
beta_mode	Optional result of beta_marginal_mode; if NULL, marginal mode is computed from design.
dispprior_list, offset, weights	Passed to beta_marginal_mode when beta_mode is NULL.
tol, maxit	Passed to beta_marginal_mode when beta_mode is NULL.
verbose	If TRUE, emit message notes about the calibrated level.

Value

An object of class "beta_marginal_safe_set" with level, beta_mode, delta_2, P_b, and certification metadata.

See Also

beta_marginal_mode, group_precision_floor, gamma_beta_tv_certificate

---

Restricted two-block Gibbs minorization certificate (Chapter C05).

Description

Computes the EM fixed point $\gamma^\star$, the mode profile $\varepsilon(\gamma^\star)$, and the spectrum-calibrated restricted-chain multiplier $\varepsilon = \varepsilon_d Q(\widetilde C_d)$.

Usage

certificate(
  design,
  pfamily_list,
  family = gaussian(),
  delta = 0.01,
  dispprior_list = NULL,
  estep = c("exact", "aghq", "mc"),
  acceleration = c("none", "squarem"),
  icm_init = TRUE,
  icm_tol = 1e-08,
  icm_maxit = 200L,
  n = 10000L,
  mc_seed = NULL,
  tol = NULL,
  em_tol = 1e-10,
  maxit = 200L
)

Arguments

design	A model_setup list.
pfamily_list	Block~2 prior list from pfamily_list().
family	A family object.
delta	Tail / escape probability budget in (0, 1).
dispprior_list	Optional Block~1 dispersion prior for gaussian().
estep	E-step tier: "exact" (default Gaussian closed form when unrestricted), "mc" (simulated conditional means; used automatically when beta_set is supplied), or "aghq" (not yet implemented).
acceleration	"none" or "squarem" (not yet implemented).
icm_init	If TRUE, run iterated conditional modes (Block~1 modes via rglmb, Block~2 via the C05 mean map) before EM. Skipped on the Gaussian exact closed-form path. Gaussian Block~1 uses exact conditional modes.
icm_tol	ICM convergence tolerance on $\|\Delta\gamma\|_{P_{11}}$.
icm_maxit	Maximum ICM iterations.
n	Number of $\beta_j$ draws per group when estep = "mc". The sample mean $\bar b_j$ has elementwise MC standard error $\mathrm{sd}(\beta_j \mid \gamma, y)/\sqrt{n}$; relative error on means is therefore $O(1/\sqrt{n})$ (e.g. n = 10000 targets $\approx 1\%$ of the posterior scale).
mc_seed	Optional seed for the MC E-step (first iteration only).
tol	If supplied, compute a sweep count n such that (1 - eps)^n + delta <= tol (Theorem 2).
em_tol	EM convergence tolerance passed to population_mode as tol (distinct from sweep-count tol below).
maxit	Maximum EM iterations.

Details

The certificate applies to the restricted Gibbs chain on $\widetilde C_d$, not necessarily the production sampler; see inst/CHAPTER_C05_IMPLEMENTATION.md section 8. Set sizing inverts the weighted-$\chi^2$ escape tail on the coupling spectrum at $\gamma^\star$ (section 4A.3). Outside Gaussian closure the exact escape mass may differ from delta.

Value

An object of class "minorization_certificate".

See Also

population_mode, epsilon, deficiency_spectrum, gamma_beta_tv_certificate

---

Check identifiability and estimability of a single-factor mixed-model design

Description

Runs the "two-step identifiability assessment" used internally by model_setup: a Level~1 (within-group) rank/estimability check on the level-1 design D, followed by a Level~2 (across-group) rank check on the level-2 design matrices W, restricted to the Level-1-estimable groups. Takes the same D/group/W matrix inputs (and the same conventions – see Details) as rLMM_reg/rGLMM_reg, so it can be run directly on hand-built matrices from those matrix-level workflows, or on model_setup()'s design$D/design$group/ design$W directly, without requiring a full lmer/glmer reference fit.

Usage

check_identifiability(
  y,
  D,
  group,
  W,
  family = gaussian(),
  group_name = attr(group, "group_name"),
  verbose = FALSE
)

Arguments

y	Response vector of length nrow(D).
D	Level-1 design matrix (nrow(D) x p_re); must have unique, non-empty colnames(D) (random-effect coefficient names).
group	Grouping factor of length nrow(D) (must be a factor); levels(group) fixes the row order assumed by W and names groupef.rank/groupef.estimable in the return value – there is no separate group_levels argument. Unlike rLMM_reg, group_name is a separate argument here (see below) rather than derived from group via substitute().
W	Named list of group-level design matrices (J x q_k), one per column of D; names must match colnames(D) exactly.
family	A family object (or a string/function accepted by family); normalized internally via the same logic model_setup uses. Defaults to gaussian().
group_name	Display name of the grouping factor (e.g. "school_id"), used only in error and verbose messages (not in any computation). Defaults to attr(group, "group_name") (the convention already used by rLMM_reg); must be resolvable to a non-empty string, either this way or by passing it explicitly. Unlike rLMM_reg, there is no fallback to substitute(group) – this function may be called on matrices that are not bare variables.
verbose	If TRUE, emit message()s summarizing the Level 1 and Level 2 outcomes. Default FALSE (silent), so calling this from model_setup introduces no new console output.

Details

Level 1 (within-group). For each group $j$, the within-group submatrix $D_j$ (D[group == $j$, ]) is checked for full column rank (groupef.rank). A rank-deficient group has too few distinct observations to estimate all random slopes independently. Algebraically full-rank groups are then additionally checked for estimability (groupef.estimable, with per-group detail in groupef.glm_check), with the specific check depending on family:

• binomial(), poisson(), Gamma(): a classical glm(y ~ D_j - 1, family) fit is attempted; a group can be full column rank yet have no finite MLE under complete or quasi-complete separation (binomial), an all-zero count response (poisson), or a non-positive response (Gamma) – Gamma groups are additionally required to have n_j > p_j so the dispersion/shape parameter is estimable.

• gaussian(): no fit is attempted (the OLS coefficient MLE always exists given full column rank), but the group must still have n_j > p_j so the residual variance/dispersion is estimable (otherwise the group achieves a perfect/saturated fit with zero residual degrees of freedom).

• other families: groupef.estimable is set equal to groupef.rank (no additional check yet).

Level 2 (across-group). Restricting to the estimable groups from Level 1 (groupef.estimable), each level-2 design matrix W[[k]] is checked for full column rank (popef.rank), using only the rows at the estimable positions (row $j$ of W[[k]] corresponds to levels(group)[j], exactly as in rLMM_reg). Rank deficiency at this level means the level-2 hyperparameters are not identified by the data, even as the number of estimable groups grows. The scalar popef.rank_ok is TRUE only when every popef.rank entry is TRUE.

Inputs must already be consistently keyed. This function does not infer or coerce missing names – it validates them up front and errors immediately on a mismatch, rather than risking a silently wrong match deeper in the rank/estimability computation:

• D must have unique, non-empty colnames(D): these are the random-effect coefficient names used to key W (there is no separate groupef.names argument to override them, matching rLMM_reg).

• group must be a factor (not coerced from a character vector); levels(group) fixes the row order assumed by W (there is no separate group_levels argument).

• W must be a named list, with names matching colnames(D) exactly (as sets; no missing, empty, or duplicate names), and each W[[k]] must have exactly length(levels(group)) rows.

• group_name, used only for error/verbose messages, must be resolvable from either the group_name argument or attr(group, "group_name").

Value

A list with:

groupef.rank

Named logical, one entry per group (names levels(group)): TRUE if $D_j$ is full column rank.

groupef.estimable

Named logical, same names as groupef.rank: TRUE if the group is additionally estimable (see Details).

groupef.glm_check

Per-group diagnostic data frame for binomial()/poisson()/Gamma()/gaussian(), or NULL for other families.

popef.rank

Named logical, one entry per RE coefficient (names colnames(D)): TRUE if W[[k]] is full column rank when restricted to groupef.estimable groups.

popef.deficient

Negation of popef.rank.

popef.rank_ok

Scalar TRUE only when every popef.rank entry is TRUE.

group_levels

levels(group), i.e. the exact names used on groupef.rank/groupef.estimable.

group_name

The resolved grouping-factor display name.

See Also

model_setup, rLMM_reg, rGLMM_reg, extract_re_hyper_matrices

---

Calibrate reference radius r(delta) and deficiency d = r^2/2.

Description

Inverts the Gaussian-reference escape tail on $B$-weights,

$\Pr\Bigl(\sum_i w_i Z_i^2 > r^2\Bigr) = \delta, \qquad w_i = \frac{\kappa_i}{1-\kappa_i},$

then sets d = r^2/2 (equivalently 2d = r^2) and set_cut = exp(-d).

Usage

deficiency_calibrate(delta, spectrum)

Arguments

delta	Tail / escape probability budget in (0, 1).
spectrum	A list from deficiency_spectrum with weights.

Details

Weighted-$\chi^2$ tails for general spectra use CompQuadForm::imhof() when the suggested package CompQuadForm is installed, otherwise mgcv::psum.chisq() if mgcv is available.

Value

A list with r2, r, d, delta, escape_mass, set_cut, and spectrum.

See Also

epsilon, deficiency_spectrum

---

Coupling eigenvalue spectrum at the population mode.

Description

Returns symmetrized coupling eigenvalues $\kappa_i = \mathrm{eig}(A)$ and deficiency weights $w_i = \kappa_i/(1-\kappa_i) = \mathrm{eig}(B)$ with $A = P_{11}^{-1/2} S P_{11}^{-1/2}$, $S = P_{11}\tilde J$, $B = A(I-A)^{-1}$ (spectrally; matrix product only when $A$ and $P_{11}$ share eigenvectors). Computed at the terminal conditional-mean E-step of population_mode.

Usage

deficiency_spectrum(mode)

Arguments

mode	A population_mode result with p11 and tilde_J.

Details

Limits: $\kappa_i \to 0 \Rightarrow w_i \to 0 \Rightarrow d \to 0$; $\kappa_i \to 1 \Rightarrow w_i \to \infty \Rightarrow d \to \infty$.

Value

A list with kappa, weights, kappa_max, rho, q, and method.

See Also

deficiency_calibrate, epsilon

---

Build a named list of dGamma measurement-dispersion priors

Description

Generic for constructing a named list of dGamma prior objects from a prior-specification object. For mixed models, one dGamma per group level supplies observation-level $\sigma^2_j$ priors for lmerb(..., group.dispersion = ...).

Usage

dGamma_list(object, ...)

Arguments

object	A prior-specification object.
...	Additional arguments passed to methods (e.g. max_disp_perc_measurement). For print.dGamma_list(), passed to print.pfamily.

Details

Where pfamily_list configures the between-group variance components $\tau^2_k$, dGamma_list() configures the within-group observation variance. It returns one dGamma prior per group level, so each group gets its own $\sigma^2_j$ rather than sharing a pooled value.

Each element places a Gamma prior on the precision $1/\sigma^2_j$ (hence Inv_Dispersion = TRUE) and carries a truncation window disp_lower/disp_upper calibrated from that group's own data. Because groups differ in size and in residual scatter, the windows differ across groups – often by an order of magnitude – which is the point of the per-group treatment.

For Prior_Setup_GLMM with a per-group dispformula, each group's dGamma() uses shape_ING/rate from object$group.ing_prior (Part VI calibration; see inst/DGAMMA_LIST_MARGINAL_AND_BOUNDS.md). Truncation bounds are posterior-shape quantiles at that group's max_disp_perc.

Value

For dGamma_list(), an object of class "dGamma_list" (also inherits from "list"): a named list of "pfamily" objects (one per group). Attributes such as window_diagnostics and measurement_prior_group are retained for downstream use but are not shown by print.dGamma_list(). For print.dGamma_list(), x invisibly.

See Also

dGamma, Prior_Setup_GLMM, pfamily_list, rLMM_reg

---

Restricted-chain Doeblin multiplier epsilon from spectrum-calibrated sizing.

Description

Combines the mode profile eps_star with the reference radius r(delta) and deficiency d = r^2/2 from the $B$-spectrum weights $w_i = \kappa_i/(1-\kappa_i)$ at mode. Produces eps_d = exp(-d) eps_star and eps = eps_d Q(C_tilde_d) with Q(C_tilde_d) at Pr(sum kappa_i Z_i^2 <= r^2).

Usage

epsilon(eps_star, delta, mode)

Arguments

eps_star	Mode profile $\varepsilon(\gamma^\star)$ from epsilon_star.
delta	Tail budget used to calibrate d.
mode	A population_mode result with p11 and tilde_J. The spectrum uses the terminal conditional-mean E-step.

Details

Weighted-$\chi^2$ tails for general coupling spectra require suggested package CompQuadForm (preferred) or mgcv.

Value

A list with eps_star, d, r, r2, eps_d, Q_mass_lb, eps, spectrum, kappa, cal, and route.

---

Numerical epsilon(gamma*) from the minorization profile.

Description

Solves the convex program of Chapter C05 Definition 5 at the EM fixed point: minimize $g(\gamma' \mid \gamma^\star) = \log q(\gamma' \mid \gamma^\star) - \log q_Q(\gamma')$ over destination $\gamma'$, with the kernel integral over $\beta$ uses an exact Gaussian convolution when family = gaussian(), otherwise Monte Carlo draws from $\pi(\beta \mid \gamma^\star, y)$.

Usage

epsilon_optimize(
  mode,
  n = NULL,
  mc_seed = NULL,
  gamma_prime_start = NULL,
  control = list(maxit = 200L),
  ...
)

Arguments

mode	Object returned by population_mode.
n	Number of $\beta$ draws for the MC inner integral when the likelihood is not Gaussian. Defaults to mode$n from population_mode, or 10000L.
mc_seed	Optional seed for the cached $\beta$ draws (fixed for the whole optimization).
gamma_prime_start	Initial destination $\gamma'$; defaults to mode$gamma_star.
control	A list passed to optim (method = "BFGS").
...	Ignored.

Value

A list with eps_star, method, certified, g_opt, gamma_prime, optim, and n.

See Also

population_mode, epsilon_star

---

Mode profile epsilon(gamma*) for the restricted Gibbs certificate.

Description

Mode profile epsilon(gamma*) for the restricted Gibbs certificate.

Usage

epsilon_star(mode, method = c("closure", "optimize"), ...)

Arguments

mode	Object returned by population_mode.
method	"closure" uses $\det(I+\tilde J)^{-1/2}$ when available (Gaussian reference); "optimize" calls epsilon_optimize.
...	Passed to epsilon_optimize when method = "optimize".

Value

A list with eps_star, method, and certified.

---

Floor coupling eigenvalues from data-precision lower bounds.

Description

Forms

$S^{\mathrm{LB}} = \sum_j H_j^\top P_b (P_{22,j}^{\mathrm{LB}})^{-1} P_b H_j$

and returns $\kappa_i^{\mathrm{LB}} = \mathrm{eig}(P_{11}^{-1/2} S^{\mathrm{LB}} P_{11}^{-1/2})$. Rosenthal drift constants are computed separately by rosenthal_drift_constants.

Usage

floor_coupling_eigenvalues(mode, floor)

Arguments

mode	A population_mode result with p11.
floor	A group_precision_floor result with Gamma_lb.

Value

A list with kappa_lb, kappa_max_lb, S_lb, and q.

See Also

rosenthal_drift_constants, rosenthal_tv_bound

---

Backward-compatible wrapper combining eigenvalues and drift constants.

Description

Backward-compatible wrapper combining eigenvalues and drift constants.

Usage

floor_coupling_spectrum(mode, beta_set)

Arguments

mode	A population_mode result.
beta_set	A group_precision_floor or legacy object with Gamma_lb.

Value

List from floor_coupling_eigenvalues plus Rosenthal drift fields from rosenthal_drift_constants.

---

Sharpest displayed gamma-beta total-variation certificate.

Description

Orchestrates categories 1–8: population_mode, epsilon_star, beta_marginal_mode, beta_marginal_safe_set, group_precision_floor, floor_coupling_eigenvalues, and optimal_rosenthal_tv_bound.

Usage

gamma_beta_tv_certificate(
  design,
  pfamily_list,
  family = gaussian(),
  delta_2 = 0.01,
  dispprior_list = NULL,
  k = NULL,
  tol = NULL,
  inner_tol = NULL,
  total_tol = NULL,
  alpha = NULL,
  display_mode = c("sharp", "general"),
  include_full_pi_gamma = TRUE,
  estep = c("exact", "aghq", "mc"),
  kappa_method = c("laplace", "crude", "none"),
  acceleration = c("none", "squarem"),
  n = 10000L,
  mc_seed = NULL,
  icm_init = TRUE,
  icm_tol = 1e-08,
  icm_maxit = 200L,
  em_tol = 1e-10,
  maxit = 200L,
  delta = NULL,
  optimize_alpha = TRUE,
  verbose = FALSE
)

Arguments

design	A model_setup list.
pfamily_list	Block~2 prior list from pfamily_list().
family	A family object.
delta_2	Asymptotic tail allowance for the full marginal.
dispprior_list	Optional Block~1 dispersion prior for gaussian().
k	If set, evaluate at this sweep count only.
tol	Alias for inner_tol when inverting for sweep count.
inner_tol	Target for the inner Rosenthal bound.
total_tol	Target for full $\pi_\gamma$ bound when include_full_pi_gamma = TRUE (uses total_tol - delta_2).
alpha	If set, skip alpha optimization.
display_mode	"sharp" (default) or "general".
include_full_pi_gamma	If TRUE, report full $\pi_\gamma$ bound as inner bound plus delta_2.
estep	E-step tier: "exact" (default Gaussian closed form when unrestricted), "mc" (simulated conditional means; used automatically when beta_set is supplied), or "aghq" (not yet implemented).
kappa_method	Passed to group_precision_floor.
acceleration	"none" or "squarem" (not yet implemented).
n	Number of $\beta_j$ draws per group when estep = "mc". The sample mean $\bar b_j$ has elementwise MC standard error $\mathrm{sd}(\beta_j \mid \gamma, y)/\sqrt{n}$; relative error on means is therefore $O(1/\sqrt{n})$ (e.g. n = 10000 targets $\approx 1\%$ of the posterior scale).
mc_seed	Optional seed for the MC E-step (first iteration only).
icm_init	If TRUE, run iterated conditional modes (Block~1 modes via rglmb, Block~2 via the C05 mean map) before EM. Skipped on the Gaussian exact closed-form path. Gaussian Block~1 uses exact conditional modes.
icm_tol	ICM convergence tolerance on $\|\Delta\gamma\|_{P_{11}}$.
icm_maxit	Maximum ICM iterations.
em_tol	EM convergence tolerance passed to population_mode as tol.
maxit	Maximum EM iterations.
delta	Tail budget for display_mode = "general" minorization via epsilon (not used in the default sharp display).
optimize_alpha	If TRUE, minimize over alpha.
verbose	If TRUE, pass verbose to beta_marginal_safe_set.

Value

An object of class "gamma_beta_tv_certificate".

See Also

certificate, optimal_rosenthal_tv_bound

---

Conditional means of group random effects given population effects

Description

Computes the C05 EM E-step: per-group conditional means $E[\beta_j \mid \gamma, y]$ (and posterior covariances $V_j$) at a supplied population location fixef. For gaussian() with estep = "exact" and no beta_set, the posterior mean is closed form. Otherwise means are Monte Carlo averages.

Usage

group_effects_conditional_mean(
  design = NULL,
  pfamily_list = NULL,
  family = NULL,
  fixef = NULL,
  mode = NULL,
  dispprior_list = NULL,
  beta_set = NULL,
  estep = c("exact", "aghq", "mc"),
  n = 10000L,
  mc_seed = NULL
)

Arguments

design	A model_setup list (required if mode is NULL).
pfamily_list	Block~2 prior list (required if mode is NULL).
family	A family object (required if mode is NULL).
fixef	Named list of population fixed effects ($\gamma$ in C05 notation). Ignored when mode is supplied.
mode	Optional result of population_mode; supplies fixef, design, p11, and model context.
dispprior_list	Optional Block~1 dispersion prior for gaussian().
beta_set	Optional beta_marginal_safe_set object ($\widetilde B(\delta_2)$); enables joint-draw rejection (see Details).
estep	"exact" (Gaussian closed form when unrestricted), "mc" (simulated means), or "aghq" (not yet implemented). Ignored when beta_set is supplied (always uses rejection MC).
n	Target number of accepted joint draws when beta_set is set; otherwise number of $\beta_j$ draws per group when estep = "mc" (required for non-Gaussian families).
mc_seed	Optional RNG seed for the MC E-step.

Details

When beta_set ($\widetilde B(\delta_2)$) is supplied, means are Monte Carlo averages of accepted joint draws from rNormal_reg_group_safe or rNormalGLM_reg_group_safe (by family).

Value

An object of class "group_effects_conditional_mean" with b_mean, V_list, b_mc_se (when MC is used), fixef, estep, n, n_target, n_tried, accept_rate (when beta_set is set), beta_set, restricted, and model context (design, family, p11).

See Also

population_mode, beta_marginal_safe_set, rNormal_reg_group_safe, rNormalGLM_reg_group_safe

---

Certified group-level data precision lower bounds

Description

Compute per-group lower-bound data precision matrices $\Gamma_j$ on a certified beta level set. Typically used after beta_marginal_mode and beta_marginal_safe_set.

Usage

group_precision_floor(
  design = NULL,
  pfamily_list = NULL,
  family = NULL,
  beta_mode = NULL,
  beta_set = NULL,
  epsilon = 0.01,
  dispprior_list = NULL,
  offset = NULL,
  weights = NULL,
  kappa_method = c("crude", "laplace", "none"),
  union_over = c("groups", "functionals"),
  mode_method = c("marginal_newton", "icm"),
  level_method = c("prop2_groups", "prop2_functionals", "r_gauss_joint"),
  tol = 1e-10,
  maxit = 200L,
  verbose = FALSE
)

Arguments

design	Used when beta_mode and beta_set are both NULL (legacy standalone path).
pfamily_list	Block~2 prior list.
family	A family object.
beta_mode	Result of beta_marginal_mode (preferred).
beta_set	Result of beta_marginal_safe_set supplying the joint level $r_{\mathrm{Gauss}}$; if NULL, epsilon and level_method define the set.
epsilon	Escape budget when building levels without beta_set.
dispprior_list	Optional Block~1 dispersion prior for gaussian().
offset, weights	Observation offset and weights; default from design.
kappa_method	"crude" (certified), "laplace", or "none".
union_over	"groups" or "functionals" (Prop 2 schemes).
mode_method	Default "marginal_newton"; ignored when beta_mode is supplied.
level_method	Default "r_gauss_joint" when beta_set is supplied; otherwise "prop2_groups" for legacy.
tol, maxit	Newton controls (unused beyond initial solve).
verbose	If TRUE, emit message notes (legacy Prop 2 path).

Value

An object of class "group_precision_floor".

See Also

beta_marginal_safe_set, floor_coupling_eigenvalues

---

Model setup for generalized linear mixed models

Description

Wrapper around lmer or glmer for models with exactly one grouping factor. Design matrices come from formula (including cross-level RE moderation terms). Random-effect and residual variance components (Psi, dispersion) come from the same reference lmer/glmer fit on formula, so printed summaries match summary(lmer(...)) on that formula.

Usage

model_setup(
  formula,
  data = NULL,
  family = gaussian(),
  REML = TRUE,
  control = NULL,
  start = NULL,
  verbose = 0L,
  subset,
  weights,
  na.action,
  offset,
  contrasts = NULL,
  devFunOnly = FALSE,
  fit_mer = TRUE,
  dispformula = ~1,
  ...
)

Arguments

formula	A two-sided linear formula, as in lmer: response on the left of ~, fixed and random-effects terms ((terms | group)) on the right. Used for both design extraction and the reference lmer/glmer fit (fixed effects, hyper calibration, and variance components); see Details for the single-grouping-factor / || constraints.
data	An optional data frame containing the variables named in formula (as in lmer); if omitted, variables are taken from the environment of formula.
family	As in glm: a character string naming a family function, a family function, or the result of a call to one. Defaults to gaussian(), in which case lmer is used. Non-Gaussian families use glmer.
REML	Logical scalar, as in lmer's own REML argument: should the reference fit's estimates optimize the REML criterion rather than the log-likelihood? Passed to lmer when family = gaussian(); ignored otherwise (glmer has no REML criterion).
control	lmerControl when family = gaussian(), otherwise glmerControl – including the nonlinear optimizer and its settings (as in lmer/glmer's own control argument); passed through to the reference fit when fit_mer = TRUE. NULL (default) resolves to lmerControl()'s defaults for Gaussian families, or glmer's own default control otherwise.
start	Optional starting values, passed to the reference lmer/glmer fit's own start argument (numeric vector or named list with par/theta, and for glmer also fixef/beta; see their documentation). NULL (default) uses lme4's own defaults.
verbose	Integer scalar passed to the reference lmer (or glmer for non-Gaussian family); see their verbose argument (> 0 prints optimizer iterations; > 1 also prints inner PIRLS steps).
subset	An optional expression indicating the subset of rows of data to use, as in lmer's subset argument. Passed to both extract_re_hyper_matrices (via lFormula) and the reference lmer/glmer (and glmmTMB when used) so the sampler design matrices and the classical reference fit share the same rows.
weights	An optional vector of prior weights, as in lmer's weights argument. Passed to design extraction (lFormula), the reference lmer/glmer fit, and glmmTMB when used, so dropped rows match subset/ na.action. Stored on the returned design as length-n weights (default all ones). Mixed-model sampler engines that still hard-code unit weights do not consume this vector yet; see Prior_Setup_GLMM.
na.action	A function indicating what should happen when data contain NAs (default na.omit, as in lmer). Passed to both design extraction (lFormula) and the reference fit so dropped rows match.
offset	An optional a priori known component to include in the linear predictor, as in lmer's offset argument (plus any offset() term in formula). Passed to design extraction and the reference fit; stored on the returned design as length-n offset (default all zeros). Mixed-model sampler engines that still hard-code a zero offset do not consume this vector yet.
contrasts	An optional list of contrasts (see the contrasts.arg of model.matrix.default, as in lmer). Passed to both design extraction (lFormula) and the reference fit so factor coding in W/D matches fixef/vcov.
devFunOnly	Logical, as in lmer's devFunOnly argument: if TRUE, return only the deviance evaluation function instead of a fitted model. Gaussian lmer fits only; ignored for glmer.
fit_mer	If TRUE (default), fit the reference lmer/glmer model and extract variance components – analogous to glmmTMB's doFit argument. If FALSE, skip the reference fit and return design matrices and rank diagnostics only (used by glmerb() in lmebayes, which fits its own glmer separately).
dispformula	As in glmmTMB's dispformula argument: a one-sided formula selecting the measurement-dispersion structure: ~1 (default, pooled) or ~<group_name> (matching the random-effects grouping factor exactly, requesting per-group dispersion). Only affects design$glmmTMB_fit (see @return); design$lmer/design$glmer are always the plain pooled-dispersion lmer/glmer fit, regardless of dispformula. ~<group_name> additionally requires family = gaussian() (no observation-level dispersion parameter to model per group otherwise). Mirrors the dispformula argument on Prior_Setup_GLMM and on lmerb()/glmerb() in lmebayes; all three are independent arguments that must be kept consistent by the caller.
...	Ignored.

Details

Uncorrelated random effects (||). The sampler treats Sigma_ranef as diagonal (no off-diagonal covariance). Multi-coefficient random terms must use ||, e.g. (1 + x || group) rather than (1 + x | group). A single random intercept may use (1 | group); (1 || group) is not supported by lme4.

Fixed-effect constraints. model_setup accepts the same formula language as lmer, subject to one structural rule: every fixed effect that does not correspond to a random-slope term must be a group-constant (level-2) covariate—a predictor whose value is the same for every observation within a given group. School-level attributes such as private_school or title1 satisfy this constraint. Student-level covariates that vary within groups may appear as fixed main effects only when they also appear as random slopes (they then represent the population mean slope $\gamma_{10}$, e.g., distracted_ppvt). Cross-level interactions of the form level2_var:random_slope (e.g., free_reduced_lunch:distracted_a1) are additionally permitted; they moderate the prior mean of the corresponding random slope across groups (see extract_re_hyper_matrices). Fixed terms that are none of these three types—level-2 covariate, population mean slope, or cross-level moderation interaction—are rejected with an informative error.

Two-step identifiability assessment. After fitting lmer/glmer, model_setup calls check_identifiability to assess whether the model is empirically identified at both the within-group (groupef.rank, groupef.estimable, groupef.glm_check) and across-group (popef.rank, popef.deficient, popef.rank_ok) levels, and copies its fields onto the returned design object – see check_identifiability for the full algorithm description (including the family-specific estimability checks and how popef.rank_ok is derived). Non-estimable groups are flagged but retained in the lmer/glmer fit; Prior_Setup_GLMM excludes them when calibrating priors, and requires popef.rank_ok = TRUE to derive default hyperpriors automatically.

The example uses big_word_club from bayesrules (see ?bayesrules::big_word_club) and the same formula as the full lmerb() demo in lmebayes (demo("Ex_12_lmerb_BigWordClub", package = "lmebayes")).

Value

Object of class "model_setup" with:

y

Response vector, length nrow(D).

weights

Numeric prior-weight vector of length length(y) (ones when unspecified), aligned with y/ D/group after subset/na.action.

offset

Numeric offset vector of length length(y) (zeros when unspecified and no formula offset()), aligned with y/D/group.

D

Group-effects model matrix (n_obs x p_re): per-observation loadings on the within-group group-effect coefficient vector (columns groupef.names).

group

Factor of length nrow(D) for block subsetting; levels(group) fixes group order used elsewhere. Also carries attr(group, "group_name") (same value as group_name) so group alone can be passed to sampler engine functions without a separate group_name argument.

group_name

Grouping factor name.

groupef.names

Group-effect coefficient names from lme4; always identical to names(W) (and to colnames(D)).

W

Named list of group-level hyper-design matrices (one per groupef.names entry, one row per group level).

popef.moderation

Cross-level moderation metadata; see extract_re_hyper_matrices.

formula

The formula argument, as supplied.

family

The normalized family object.

call

The matched call.

terms

The terms object for formula, as in lm/glm (and glmbayes::glmb/lmb).

data

The data argument, as supplied (NULL when omitted), as in lm/glm (and glmbayes::glmb/lmb).

dispformula

The dispformula argument, as supplied.

lmer / glmer

The reference lmer (Gaussian) or glmer (otherwise) fit on formula; NULL when fit_mer = FALSE.

varcorr

The raw VarCorr object from the reference fit.

Psi

Named numeric vector (names groupef.names) of group-effect variances.

dispersion

Residual (measurement) variance from the reference fit.

glmmTMB_fit

NULL when dispformula = ~1 or fit_mer = FALSE; otherwise a glmmTMB fit on formula with the same dispformula (Gaussian models only), fit via the same data/REML/ control-derived arguments as lmer. This is additive: lmer/glmer are never replaced by a glmmTMB fit. Prior_Setup_GLMM and lmerb() (in lmebayes) reuse glmmTMB_fit as their per-group-dispersion calibration reference instead of fitting glmmTMB a second time.

groupef.rank

Named logical (names levels(group)): TRUE if that group's $D_j$ is full column rank; see check_identifiability.

groupef.estimable

Named logical, same names as groupef.rank: TRUE if the group is additionally estimable (family-specific MLE/dispersion check); see check_identifiability.

groupef.glm_check

Per-group diagnostic data frame from check_identifiability (family-dependent; may be NULL).

popef.rank

Named logical (names groupef.names): TRUE if W[[k]] is full column rank when restricted to groupef.estimable groups.

popef.deficient

Negation of popef.rank.

popef.rank_ok

Scalar TRUE only when every popef.rank entry is TRUE; required by Prior_Setup_GLMM to derive default hyperpriors automatically.

See Also

check_identifiability, extract_re_hyper_matrices

---

Legacy MVN-calibrated reference radius (equal-weight chi-square shortcut).

Description

Computes d = 0.5 * qchisq(1 - delta, q). This equal-weight route is retained for reference only; epsilon and certificate use deficiency_calibrate on the coupling spectrum from deficiency_spectrum instead.

Usage

mvn_calibrate(q, delta)

Arguments

q	Hyperparameter dimension.
delta	Tail / escape probability budget in (0, 1).

Value

A list with r2, r, d, and set_cut = exp(-d).

See Also

deficiency_calibrate, epsilon

---

Optimize Rosenthal TV bound over alpha and/or sweep count k.

Description

Optimize Rosenthal TV bound over alpha and/or sweep count k.

Usage

optimal_rosenthal_tv_bound(
  mode,
  eigenvalues,
  eps,
  k = NULL,
  alpha = NULL,
  inner_tol = NULL,
  total_tol = NULL,
  delta_2 = 0,
  include_full_pi_gamma = TRUE,
  optimize_alpha = TRUE,
  display_mode = c("sharp", "general"),
  d = NULL,
  V_sup = NULL
)

Arguments

mode	A population_mode result.
eigenvalues	Output of floor_coupling_eigenvalues.
eps	Minorization constant for the bound.
k	If set, evaluate at this sweep count only.
alpha	If set, skip alpha optimization.
inner_tol	Target for the inner Rosenthal bound.
total_tol	Target for full $\pi_\gamma$ bound when include_full_pi_gamma = TRUE (uses total_tol - delta_2).
delta_2	Asymptotic tail allowance for the full marginal.
include_full_pi_gamma	If TRUE, set full_bound to inner bound plus delta_2.
optimize_alpha	If TRUE, minimize over alpha.
display_mode	"sharp" (default) or "general".
d	Deficiency level for $V_{\sup}(d)$ when V_sup is omitted.
V_sup	Override for $V_{\sup}$ in general mode.

Value

A list of class "optimal_rosenthal_tv_bound".

See Also

rosenthal_tv_bound, gamma_beta_tv_certificate

---

Build a named list of pfamily objects from a prior-setup object

Description

Generic for converting a prior-setup object into a named list of pfamily objects (class "pfamily_list") for block or mixed-model samplers.

Usage

pfamily_list(object, ptypes = NULL, ...)

Arguments

object	A prior-setup object. Dispatch is method-specific; see Methods and linked help pages for required classes and arguments.
ptypes	Prior-family name(s) per component; meaning depends on the active method (for row-block setup see Prior_SetupGroup).
...	Passed to methods (unused for Prior_SetupGroup).

Details

For a Prior_Setup_GLMM object, each group-effect coefficient $k$ uses parameters from object$pop.prior_list[[k]]:

• "dNormal": dNormal(mu, Sigma, dispersion) with group-effect variance $\tau^2_k$ from pop.prior_list[[k]]$dispersion (A12 §3.3.4 marginal hyper-regression center for Gaussian models; classical VarCorr in pop.dispersion.ref).

• "dIndependent_Normal_Gamma": the same mu and Sigma, plus a Gamma prior on precision $1/\tau^2_k$ from object$pop.ing_prior[[k]] (calibrated by Prior_Setup_GLMM via pop.dispersion.pwt / pop.dispersion.nprior, mean-matched at that same marginal $\tau^2_k$). With effective prior sample size $n_0$ and $p_k$ population predictors for component $k$,

  $shape = (n_0 + 1 + p_k)/2, \qquad rate = \tau^2_k \, (n_0 + p_k - 1)/2.$

  Requires $n_0 \le J$ (equivalently pop.dispersion.pwt $\le 0.5$). Truncation bounds disp_lower / disp_upper come from the limiting posterior window stored on pop.ing_prior; see inst/ING_TRUNCATION_WINDOW.md.

This is separate from observation-level $\sigma^2$ priors: those use dGamma_list (per-group dispformula) or a single dGamma built from object$group.ing_prior when dispersion is pooled.

Value

For pfamily_list(), an object of class "pfamily_list" (also inherits from "list"): a named list of pfamily objects. Attribute "ptypes" records the resolved prior-family name per component when the method sets it. For print.pfamily_list(), x invisibly.

See Also

Prior_SetupGroup, Prior_Setup_GLMM, pfamily, dNormal, dNormal_Gamma, dIndependent_Normal_Gamma

---

Combined mean-convergence plot for Block~2 fixed effects (Claim 1)

Description

Tracks convergence of the Block~2 fixed-effects (hyperparameter) posterior mean across inner Gibbs sweeps as a single combined chart (all traces overlaid on one axis), the mean-side counterpart of plot_var_convergence's variance-side Claim~3 chart. Per the two-block Gibbs ergodicity result's Claim~1, the cross-chain mean bias $E[\gamma^{(l)}] - \gamma^*$ decays like $A^l$ (linear in l, vs. Claim~3's $A^{2l}$ for variance). Each coefficient's per-sweep deviation from its reference mean is standardized by that same sweep's cross-chain SD, $Z(l) = (\mathrm{mean}(l) - \mathrm{mean}_{ref}) / \mathrm{sd}(l)$, so coefficients on very different scales can be overlaid on one chart and Z(l) is expected to shrink toward 0 as l grows.

Usage

plot_mean_convergence(hist, ...)

## Default S3 method:
plot_mean_convergence(
  hist,
  coef_focus = NULL,
  component = c("fixef", "precision"),
  design = NULL,
  measurement_prior_list = NULL,
  whitened = FALSE,
  engine = c("base", "ggplot"),
  n_chains = NULL,
  conf_level = 0.95,
  stage_label = hist$stage,
  split = c("auto", "none"),
  max_whitened = 4L,
  ...
)

## S3 method for class 'rLMMNormal_reg'
plot_mean_convergence(
  hist,
  coef_focus = NULL,
  component = c("fixef", "precision"),
  whitened = FALSE,
  engine = c("base", "ggplot"),
  n_chains,
  conf_level = 0.95,
  stage = c("main", "pilot"),
  stage_label = NULL,
  ...
)

## S3 method for class 'rLMMindepNormalGamma_reg'
plot_mean_convergence(
  hist,
  coef_focus = NULL,
  component = c("fixef", "precision"),
  whitened = FALSE,
  engine = c("base", "ggplot"),
  n_chains,
  conf_level = 0.95,
  stage = c("main", "pilot"),
  stage_label = NULL,
  ...
)

## S3 method for class 'rlmerb'
plot_mean_convergence(
  hist,
  coef_focus = NULL,
  component = c("fixef", "precision"),
  whitened = FALSE,
  engine = c("base", "ggplot"),
  n_chains,
  conf_level = 0.95,
  stage = c("main", "pilot"),
  stage_label = NULL,
  ...
)

## S3 method for class 'rGLMM_reg'
plot_mean_convergence(
  hist,
  coef_focus = NULL,
  component = c("fixef", "precision"),
  whitened = FALSE,
  engine = c("base", "ggplot"),
  n_chains,
  conf_level = 0.95,
  stage = c("main", "pilot"),
  stage_label = NULL,
  ...
)

## S3 method for class 'rglmerb'
plot_mean_convergence(
  hist,
  coef_focus = NULL,
  component = c("fixef", "precision"),
  whitened = FALSE,
  engine = c("base", "ggplot"),
  n_chains,
  conf_level = 0.95,
  stage = c("main", "pilot"),
  stage_label = NULL,
  ...
)

Arguments

hist, coef_focus, component, engine, n_chains, conf_level, split, max_whitened	Same meaning as in plot_var_convergence.
...	Passed to methods; see plot_mean_convergence.default.
design, measurement_prior_list	When both are supplied, the exact reference mean is resolved automatically via lmerb_posterior_mean and, for whitened = TRUE, the exact reference covariance (via lmerb_posterior_covariance) needed to whiten the mean deviation vector; see plot_var_convergence for when these are available.
whitened	When FALSE (default), plots Z(l) per named coefficient in coef_focus. When TRUE, whitens each sweep's cross-chain mean deviation vector by the reference covariance ($\Sigma_{ref}^{-1/2} (\mathrm{mean}(l) - \mathrm{mean}_{ref})$) and plots its components – basis-invariant – labeled z1, z2, ... instead of named coefficients.
stage_label	Defaults to the resolved hist$stage; pass explicitly to override.
stage	"main" (default) or "pilot"; see plot_var_convergence.rLMMNormal_reg.

Details

What is plotted

The x-axis is the inner Gibbs sweep index $l$, and each point is a cross-chain statistic computed over the n_chains independent replicate chains, all of which sit at sweep $l$ simultaneously. Nothing here is a running statistic along a single chain.

The plotted quantity is the standardized deviation

$Z(l) = \frac{\mathrm{mean}(l) - \mathrm{mean}_{ref}}{\mathrm{sd}(l)},$

which divides out the scale of each coefficient so that an intercept measured in hundreds and a slope measured in tenths can share one axis. A converged run shows every trace scattered around 0 with no trend; a systematic drift toward or away from 0 means the sweep budget is too small. Supplying n_chains adds a naive band that separates Monte Carlo noise from genuine bias.

mean_ref is the exact posterior mean when design and measurement_prior_list are supplied (only possible when $\sigma^2$ and the RE variance components are all fixed), and otherwise the empirical last-sweep mean. As with plot_var_convergence, the empirical fallback cannot detect a run that never converged, because the reference is drawn from the same unconverged distribution.

Reading it alongside the variance chart

This chart is the location-side counterpart of plot_var_convergence, which tracks spread. The variance converges at twice the rate of the mean, so this is the slower — and therefore the more informative — of the two when deciding whether m_convergence is large enough. A fit whose variance ratios have settled at 1 while these traces still drift has converged in spread but not in location.

Because every replicate chain starts at the posterior mean, the bias begins at zero on Gaussian fits with fixed variance components: expect a flat chart from sweep 1 there. Read a non-flat one as evidence that something upstream — an estimated variance component, or a start at the mode rather than the mean — has moved the starting point off the target mean. vignette("Chapter-B02", package = "lmebayesCore") derives the two rates.

When there is nothing to plot

As for plot_var_convergence: the exact route runs no sweeps to record, some Gaussian paths do not export per-sweep cross-chain statistics, and whitened = TRUE needs slots that only the sweep-recording engines populate. See that topic, which also explains how to obtain n_chains from a fitted object.

Value

hist invisibly.

See Also

plot_var_convergence, lmerb_posterior_mean

---

Combined variance-convergence plot for Block~2 fixed effects (Claim 3)

Description

Tracks convergence of the Block~2 fixed-effects (hyperparameter) posterior variance across inner Gibbs sweeps as a single combined chart (all traces overlaid on one axis, following lmebayes::Ex_16's style), per the two-block Gibbs ergodicity result (Claim~3 of the package's ergodicity reference, inst/BLOCK_GIBBS_ERGODICITY.md): after l sweeps, the cross-chain covariance $\Sigma^{(l)}_{11}$ satisfies

$\Sigma_{11}^{-1/2} \Sigma^{(l)}_{11} \Sigma_{11}^{-1/2} = I - A^{2l}$

for a convergent (symmetric, eigenvalues in $[0, 1)$) matrix $A$. Var(l) / Var_ref ratios are therefore bounded above by 1 and increase toward it as l grows.

Usage

plot_var_convergence(hist, ...)

## Default S3 method:
plot_var_convergence(
  hist,
  coef_focus = NULL,
  component = c("fixef", "precision"),
  design = NULL,
  measurement_prior_list = NULL,
  whitened = FALSE,
  engine = c("base", "ggplot"),
  n_chains = NULL,
  conf_level = 0.95,
  stage_label = hist$stage,
  split = c("auto", "none"),
  max_whitened = 4L,
  ...
)

## S3 method for class 'rLMMNormal_reg'
plot_var_convergence(
  hist,
  coef_focus = NULL,
  component = c("fixef", "precision"),
  whitened = FALSE,
  engine = c("base", "ggplot"),
  n_chains,
  conf_level = 0.95,
  stage = c("main", "pilot"),
  stage_label = NULL,
  ...
)

## S3 method for class 'rLMMindepNormalGamma_reg'
plot_var_convergence(
  hist,
  coef_focus = NULL,
  component = c("fixef", "precision"),
  whitened = FALSE,
  engine = c("base", "ggplot"),
  n_chains,
  conf_level = 0.95,
  stage = c("main", "pilot"),
  stage_label = NULL,
  ...
)

## S3 method for class 'rlmerb'
plot_var_convergence(
  hist,
  coef_focus = NULL,
  component = c("fixef", "precision"),
  whitened = FALSE,
  engine = c("base", "ggplot"),
  n_chains,
  conf_level = 0.95,
  stage = c("main", "pilot"),
  stage_label = NULL,
  ...
)

## S3 method for class 'rGLMM_reg'
plot_var_convergence(
  hist,
  coef_focus = NULL,
  component = c("fixef", "precision"),
  whitened = FALSE,
  engine = c("base", "ggplot"),
  n_chains,
  conf_level = 0.95,
  stage = c("main", "pilot"),
  stage_label = NULL,
  ...
)

## S3 method for class 'rglmerb'
plot_var_convergence(
  hist,
  coef_focus = NULL,
  component = c("fixef", "precision"),
  whitened = FALSE,
  engine = c("base", "ggplot"),
  n_chains,
  conf_level = 0.95,
  stage = c("main", "pilot"),
  stage_label = NULL,
  ...
)

Arguments

hist	Object of class "two_block_sweep_history" (see print.two_block_sweep_history). For whitened = TRUE, hist must also carry cov_by_sweep/coef_index, which are only populated by sweeps-outer/chains-inner engines – i.e. every engine driven by rGLMM_sweep: the ⁠_run_with_pilot()⁠ family (rLMMNormal_reg_estimated_vcov(), rLMMindepNormalGamma_reg_known_vcov(), rLMMindepNormalGamma_reg_estimated_vcov()) and, as of the rGLMM_sweep()-based rewrite of its engine, rLMMNormal_reg_known_vcov(sim_method = "TWO_BLOCK_GIBBS") as well.
...	Passed to methods (e.g. coef_focus, design, whitened, n_chains, ...; see below). Ignored by plot_var_convergence.default beyond the documented formals.
coef_focus	Optional list of length-2 character vectors c(re_component, covariate) to restrict the plot to; default (NULL) uses every coefficient tracked in hist.
component	"fixef" (default) plots the Block~2 fixed-effects/hyperparameter series described above. "precision" instead plots the estimated random-effects precision ($1/\tau^2_k$, i.e. the diagonal of $\Sigma_{ranef}^{-1}$) for every ING (dIndependent_Normal_Gamma) Block~2 component – requires hist$disp_table (populated whenever at least one Block~2 component has an estimated, not fixed, variance – see fit$any_non_normal); errors if that is empty (e.g. a known-vcov fit) or if whitened = TRUE (only the cross-chain variance of each component's precision is captured so far, not its covariance with other components/fixed effects, so whitening is not yet possible – see inst/BLOCK_GIBBS_ERGODICITY.md for the deferred extension). Precision is plotted (not the variance $\tau^2_k$ itself) because $E[1/\tau^2_k]$ stays well-defined under a weak/vague prior where $E[\tau^2_k]$ need not be. This is a distinct quantity from the observation-level dispersion estimated by rLMMindepNormalGamma_reg_*() (despite both being called "dispersion" internally) – component = "precision" is always about the RE variance-covariance, never the observation noise variance. design/measurement_prior_list are ignored for component = "precision": no exact reference exists yet for an estimated $\Sigma_{ranef}$, so the denominator is always the empirical last-sweep value.
design, measurement_prior_list	When both are supplied, the exact reference covariance is resolved automatically via lmerb_posterior_covariance and used as the denominator – this is not an opt-in toggle, it is used whenever the ingredients for it exist. Only valid when dispersion and the RE variance-covariance are both fixed (not sampled). When omitted (the only option for estimated dispersion or estimated-vcov models), the denominator falls back to the empirical last-sweep cross-chain variance/covariance ("Var_final").
whitened	When FALSE (default), plots Var(l) / Var_ref per named coefficient in coef_focus. When TRUE, whitens each sweep's full cross-chain covariance by the reference covariance ($\Sigma_{ref}^{-1/2} \mathrm{Cov}(l) \Sigma_{ref}^{-1/2}$) and plots its eigenvalues – basis-invariant, so no eigenvectors of the theoretical $A$ are ever needed – labeled var1, var2, ... instead of named coefficients (per Claim~3, these eigen-directions mix the original coefficients together, so named labels no longer apply).
engine	"base" (default) or "ggplot" (requires ggplot2); both draw one combined chart with all traces overlaid.
n_chains	Optional integer, the number of independent chains each plotted Var(l) (and Var_ref, when it is the empirical Var_final rather than an exact reference) was computed across – i.e. the same n_chains/n_pilot/n passed to the sweep-generating call (see Details). When supplied, draws a single pair of horizontal dotted reference lines: a naive conf_level band around 1 for Ratio(l), under the null that the true ratio at that sweep already equals 1 (unrealistic at small l, but useful as a ruler for how much of the departure from 1 is just Monte Carlo noise). The band's distribution depends on whether Var_ref is exact or empirical (see design/measurement_prior_list above): ref_source = "exact": only Var(l) (the numerator) has sampling error, so $\mathrm{Ratio}(l) \sim \chi^2_{n_{chains}-1} / (n_{chains}-1)$. ref_source = "empirical (Var_final)": Var_ref itself is a sample variance/covariance across the same n_chains chains (just at the final sweep), so both numerator and denominator carry sampling error and the classic two-sample variance-ratio result applies instead: $\mathrm{Ratio}(l) \sim F(n_{chains}-1,\, n_{chains}-1)$. This is wider than the chi-squared band, and reduces to it exactly as Var_final's own degrees of freedom go to infinity (stats::qf(p, df1, Inf) == stats::qchisq(p, df1) / df1), so both cases are computed via a single stats::qf() call with df2 = Inf for the exact case. Not series-specific (a single band, not one per coefficient/eigenvalue) since it depends only on n_chains/ref_source. For whitened = TRUE this is an additional approximation regardless of ref_source: whitened eigenvalues are eigenvalues of a Wishart-distributed matrix, not marginally chi-squared/F, and tend to spread out more than this band implies (eigenvalue repulsion), so treat it as an illustrative ruler, not an exact whitened confidence band. Also note that Var(l) and Var_final come from the same n_chains chains at two different sweep indices, not independent samples the way a textbook two-sample F-test assumes – another reason this is a naive/illustrative band rather than an exact one. Omit (the default, NULL) to skip the band entirely.
conf_level	Confidence level for the n_chains band (default 0.95); ignored when n_chains is NULL.
stage_label	Defaults to the resolved hist$stage; pass explicitly to override.
split	"auto" (default) or "none". "auto" draws one separate chart – its own plot page/window, never stacked via layout/par(mfrow = ...) on top of another – per group: for whitened = FALSE, one chart for every coefficient whose re_component is "(Intercept)" ("Intercept predictors") and one for every other re_component combined ("Slope predictors"), whichever of the two is non-empty; for whitened = TRUE, consecutive batches of at most max_whitened eigenvalues each. "none" restores the single-combined-chart behaviour (every tracked coefficient/eigenvalue overlaid on one panel), same as before split existed.
max_whitened	Maximum number of eigenvalue series per chart when whitened = TRUE and split = "auto" (default 4); ignored otherwise.
stage	"main" (default) or "pilot" – which sweep history to plot (fit$sweep_history or fit$pilot$draws$sweep_history). Ignored by plot_var_convergence.default, which always takes hist literally.

Details

What is plotted

The x-axis is the inner Gibbs sweep index $l$, not the stored-draw index: these charts diagnose convergence within the burn-in that precedes each stored draw. Each point is a cross-chain statistic at sweep $l$ — the variance across the n_chains independent replicate chains, all of which are at sweep $l$ simultaneously — not a running statistic along a single chain. This is why the sweep-outer loop order matters: it is what makes "the sampler's distribution after $l$ sweeps" an observable quantity.

The y-axis is the ratio Var(l) / Var_ref. Under Claim 3 that ratio is bounded above by 1 and rises toward it, so a converged run shows every trace flattening just below 1. Traces that are still climbing at the last sweep indicate the sweep budget is too small; traces that overshoot 1 substantially are Monte Carlo noise, and n_chains draws the band that says how much overshoot is expected.

Var_ref is the exact posterior variance when design and measurement_prior_list are both supplied (only possible when $\sigma^2$ and the RE variance components are all fixed), and otherwise the empirical last-sweep variance. The empirical fallback cannot detect a run that has not converged at all, since the denominator is then drawn from the same unconverged distribution; prefer the exact reference whenever the model allows it.

Reading it alongside the mean chart

Location and spread converge at different speeds, and the two plotting functions separate them: plot_mean_convergence tracks the mean, plot_var_convergence the variance. The variance converges at twice the rate, so the mean chart settles later and is usually the one that decides whether the sweep budget is adequate. A run whose variance ratios sit at 1 while its standardized mean deviations still drift has converged in spread but not in location.

Because every chain starts at the posterior mean, the mean chart is flat from sweep 1 on well-behaved Gaussian fits, leaving the variance chart to carry the information. vignette("Chapter-B02", package = "lmebayesCore") derives the two rates.

When there is nothing to plot

sweep_history is recorded only when all chains advance together, one sweep at a time, so that the cross-chain spread at sweep $l$ is meaningful. Two situations produce no history:

• The exact route (fixed $\sigma^2$, all components dNormal, sim_method = "DEFAULT") runs no sweeps at all, because the draws are already exact. There is nothing to diagnose. Passing sim_method = "TWO_BLOCK_GIBBS" forces the iterative path on the same target if you want a history anyway.

• Some Gaussian paths run sweeps but do not export per-sweep cross-chain statistics; fit$sweep_history is then NULL.

whitened = TRUE additionally requires cov_by_sweep and coef_index, which only the sweep-outer engines above populate, and component = "precision" additionally requires disp_table, which requires at least one ING Block 2 component.

Supplying n_chains

n_chains is not stored on hist itself, so it must be supplied from the fitted object that produced it: for the main-stage history (fit$sweep_history) use nrow(fit$popef[[k]]) for any RE component k (all have the same number of rows, the n passed to the ⁠_run_with_pilot()⁠ call); for the pilot-stage history (fit$pilot$draws$sweep_history) use fit$pilot$n.

Value

hist invisibly.

See Also

plot_mean_convergence, lmerb_posterior_covariance

---

EM fixed point for the C05 population mode gamma*.

Description

Finds the fixed point of the conditional-mean map (Chapter C05 Stage 1). Uses conditional means for Block~1, not modes.

Usage

population_mode(
  design,
  pfamily_list,
  family = gaussian(),
  dispprior_list = NULL,
  beta_set = NULL,
  estep = c("exact", "aghq", "mc"),
  acceleration = c("none", "squarem"),
  n = 10000L,
  mc_seed = NULL,
  icm_init = TRUE,
  icm_tol = 1e-08,
  icm_maxit = 200L,
  icm_screen = NULL,
  mc_alpha = 0.05,
  mc_stable = 2L,
  tol = 1e-10,
  maxit = 200L,
  verbose = FALSE
)

Arguments

design	A model_setup list.
pfamily_list	Block~2 prior list from pfamily_list().
family	A family object.
dispprior_list	Optional Block~1 dispersion prior for gaussian().
beta_set	Optional beta_marginal_safe_set; enables widetilde-B rejection MC (see Details).
estep	E-step tier: "exact" (default Gaussian closed form when unrestricted), "mc" (simulated conditional means; used automatically when beta_set is supplied), or "aghq" (not yet implemented).
acceleration	"none" or "squarem" (not yet implemented).
n	Number of $\beta_j$ draws per group when estep = "mc". The sample mean $\bar b_j$ has elementwise MC standard error $\mathrm{sd}(\beta_j \mid \gamma, y)/\sqrt{n}$; relative error on means is therefore $O(1/\sqrt{n})$ (e.g. n = 10000 targets $\approx 1\%$ of the posterior scale).
mc_seed	Optional seed for the MC E-step (first iteration only).
icm_init	If TRUE, run iterated conditional modes (Block~1 modes via rglmb, Block~2 via the C05 mean map) before EM. Skipped on the Gaussian exact closed-form path. Gaussian Block~1 uses exact conditional modes.
icm_tol	ICM convergence tolerance on $\|\Delta\gamma\|_{P_{11}}$.
icm_maxit	Maximum ICM iterations.
icm_screen	If TRUE (default when the E-step is Monte Carlo), after ICM run mc_stable MC mean-map fixed-point checks at the ICM fixef; skip the EM loop when all pass. If NULL, enabled whenever icm_init is TRUE and the E-step uses MC.
mc_alpha	Two-sided level for the MC fixed-point screen and EM stop (tol_eff = max(tol, z_{1-alpha/2} * mc_delta_floor)).
mc_stable	Number of consecutive MC fixed-point passes required at the same fixef (ICM screen and EM stopping rule).
tol	Convergence tolerance on the Mahalanobis change in fixef under the C05 metric $\|\Delta\gamma\|_{P_{11}}$. When the E-step is Monte Carlo, the effective tolerance is max(tol, z_{1-alpha/2} * mc_delta_floor).
maxit	Maximum EM iterations.
verbose	If TRUE, emit message() diagnostics for ICM initialization and each EM record (gradient $g$, $\|g\|$, MC $\log \pi(\gamma \mid y)$ when joint safe draws are available).

Details

In the restricted certificate, the EM E-step uses conditional means under the $\beta$-truncated target on $\widetilde B(\delta_2)$. Pass Pass beta_set from beta_marginal_safe_set; routes to rNormal_reg_group_safe (Gaussian) or rNormalGLM_reg_group_safe (non-Gaussian).

Value

A list with fixef ($\gamma^\star$), gamma_star, b_mean, b_mc_se (terminal MC standard errors when estep = "mc"), V_list, n, model context (design, family, measurement_prior_list, p11), refresh objects (P11, Sigma_star), tilde_J, kappa_spectrum, weights, kappa, rho, Sigma_pi, eps_star_closure, icm and icm_screen diagnostics, EM diagnostics (tol_eff, mc_delta_floor, em_route, em_trace), beta_set, restricted, and related fields.

See Also

group_effects_conditional_mean, rNormal_reg_group_safe, rNormalGLM_reg_group_safe

---

Print simulated group coefficients from an LMM/GLMM result

Description

S3 generic for distribution-sampler style display of groupef draws (long data frame), parallel to print() for population coefficients. Filters by grouping levels and optional RE columns. Methods are registered for rLMMNormal_reg, rLMMindepNormalGamma_reg, rGLMM_reg, rlmerb, and rglmerb; the default method accepts any list-like object with a usable groupef data frame.

Usage

print_groupef(
  x,
  groups = NULL,
  components = NULL,
  draws = NULL,
  digits = max(3, getOption("digits") - 3),
  ...
)

## Default S3 method:
print_groupef(
  x,
  groups = NULL,
  components = NULL,
  draws = NULL,
  digits = max(3, getOption("digits") - 3),
  ...
)

## S3 method for class 'rLMMNormal_reg'
print_groupef(
  x,
  groups = NULL,
  components = NULL,
  draws = NULL,
  digits = max(3, getOption("digits") - 3),
  ...
)

## S3 method for class 'rLMMindepNormalGamma_reg'
print_groupef(
  x,
  groups = NULL,
  components = NULL,
  draws = NULL,
  digits = max(3, getOption("digits") - 3),
  ...
)

## S3 method for class 'rlmerb'
print_groupef(
  x,
  groups = NULL,
  components = NULL,
  draws = NULL,
  digits = max(3, getOption("digits") - 3),
  ...
)

## S3 method for class 'rGLMM_reg'
print_groupef(
  x,
  groups = NULL,
  components = NULL,
  draws = NULL,
  digits = max(3, getOption("digits") - 3),
  ...
)

## S3 method for class 'rglmerb'
print_groupef(
  x,
  groups = NULL,
  components = NULL,
  draws = NULL,
  digits = max(3, getOption("digits") - 3),
  ...
)

Arguments

x	An object with a groupef data frame (e.g. an rLMM_reg / rGLMM_reg return, or rlmerb / rglmerb).
groups	NULL (default; all group levels), a character vector of levels in the grouping column, or integer indices into the unique levels present in groupef.
components	NULL (all RE columns) or a character vector / integer indices selecting RE coefficient columns (excluding the draw and grouping columns).
draws	NULL (all draws) or integer indices / logical of length n selecting which draw values to display (print-only).
digits	Number of significant digits for printed numeric values.
...	further arguments passed to or from other methods.

Details

The groupef layout

Population coefficients ($\gamma$) and group coefficients ($\beta_j$) are stored in different shapes, which is why they have different printers. popef is a named list of matrices, one per random-effect component, each n draws by $q_k$ level-2 predictors — printed by the print() methods. groupef is a single long-format data frame with one row per (draw, group level) pair:

draw	integer draw index, 1:n
grouping column	the group level, named after the grouping factor ("Subject", "state", ...)
one column per RE coefficient	"(Intercept)" and any random slopes

So a fit with n = 1000 draws over 20 groups produces 20,000 rows. Printing it whole is rarely what anyone wants, which is the reason this function exists.

Filtering

The three filters are independent and compose; all of them are print-only, and x is never modified.

• draws — integer indices into 1:n, or a logical mask of length n. Filters rows by the draw column.

• groups — level names as characters, or integer indices into the unique levels present in groupef (which is the ordering the fit used, not necessarily levels() of the original factor). Filters rows by the grouping column.

• components — names or integer indices selecting which random-effect columns to show. The draw and grouping columns are always kept and are never counted in the indexing.

print_groupef(fit, draws = 1:5, groups = 1:3) is the usual idiom: the first few draws for the first few groups, enough to confirm the sampler is producing sensible numbers without flooding the console. An empty result after filtering prints a placeholder rather than erroring.

Numeric columns are formatted to digits significant figures for display only; the stored values keep full precision.

Dispatch

All registered methods share the default implementation — they exist so that the generic dispatches cleanly on every sampler return class rather than falling through, and so that class-specific layouts can be added later without changing call sites. The default method works on any list-like object carrying a usable groupef data frame, which includes intermediate objects produced during route dispatch.

Point-estimate fits (simulate = FALSE) have no groupef at all; use groupef.mode directly for those.

Value

No return value, called for side effects.

See Also

rLMM_reg, rGLMM_reg, rlmerb, rglmerb

---

Print a population_mode_em_trace object.

Description

Print a population_mode_em_trace object.

Usage

## S3 method for class 'population_mode_em_trace'
print(x, ...)

Arguments

x	A population_mode_em_trace object (from population_mode$em_trace).
...	Ignored.

Value

x, invisibly.

---

Prior setup for generalized linear mixed models

Description

Calibrates priors for the level-2 fixed effects (fixef) of a hierarchical mixed model using the reference lmer/glmer fits on all groups (from model_setup). Per-group design rank (groupef.rank) and estimability (groupef.estimable – classical-glm MLE existence for family = binomial()/poisson()/Gamma(), residual degrees of freedom for family = gaussian()) are diagnostic checks computed by model_setup and do not subset the data or the reference glmer fit. Random-effect variances are treated as fixed at their mixed-model estimates. The returned object organizes calibrated priors as two independent blocks:

Usage

Prior_Setup_GLMM(
  formula,
  data,
  family = gaussian(),
  REML = TRUE,
  control = NULL,
  start = NULL,
  verbose = 0L,
  subset,
  weights,
  na.action,
  offset,
  contrasts = NULL,
  dispformula = ~1,
  pop.pwt = NULL,
  pop.pwt_default_low = 0.01,
  pop.pwt_default_high = 0.05,
  pop.nprior = NULL,
  pop.sd = NULL,
  pop.dispersion.pwt = NULL,
  pop.dispersion.nprior = NULL,
  pop.max_disp_perc = 0.99,
  pop.intercept_source = c("null_model", "full_model"),
  pop.effects_source = c("null_effects", "full_model"),
  pop.mu = NULL,
  group.dispersion.pwt = NULL,
  group.dispersion.nprior = NULL,
  group.alpha_target = 0.01,
  group.max_disp_perc = 0.99
)

Arguments

formula	Mixed-model formula passed to model_setup, whose reference lmer/glmer fits (all groups) supply the calibration quantities.
data	Data frame containing all variables in formula.
family	Model family. Default gaussian(). Non-Gaussian families use glmer for calibration; Block~1 dispersion Gamma calibration is omitted (analogous to Prior_Setup for flat GLMs).
REML	Logical scalar passed to model_setup (and thus to the Gaussian reference lmer / glmmTMB fits). Default TRUE. Ignored for non-Gaussian families (glmer has no REML criterion). Affects the classical estimates used for prior calibration only.
control	Optional lmerControl / glmerControl object passed to model_setup. When NULL (default), Gaussian models use lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5)) and non-Gaussian models leave glmer's own control default.
start	Optional starting values passed to model_setup for the reference lmer/glmer fit. Default NULL.
verbose	Integer scalar passed to model_setup for the reference lmer/glmer optimizer. Default 0L.
subset, na.action	Optional arguments passed through to model_setup (same meanings as in lmer / glmer). model_setup() applies them to both design extraction (lFormula) and the reference lmer/glmer/glmmTMB fits, so sampler matrices (design$y/D/W) and classical calibration use the same rows.
weights, offset	Optional arguments passed through to model_setup (same meanings as in lmer/glmer). Stored on design$weights / design$offset as length-n vectors aligned with design$y (after subset/ na.action), for Phase~1 availability to samplers. Not fully active for mixed-model sampling yet. They affect the classical reference fit used for prior calibration. Low-level two_block_rNormal_reg already accepts and echoes prior.weights/offset2 when passed explicitly; high-level mixed routes that still hard-code unit weights / zero offset (and ICM / prior-weight algebra that assume unit weights) do not consume design$weights/design$offset yet. Prefer incorporating an offset in the formula or waiting until those paths are wired before relying on these for the Bayesian sampler.
contrasts	Optional contrast list passed to model_setup (as in lmer). Applied to both design extraction and the reference fit so factor coding in W matches fixef/vcov.
dispformula	One-sided formula selecting the Block~1 measurement- dispersion structure: ~1 (default, pooled) or ~<group_name> (matching the random-effects grouping factor exactly, requesting per-group dispersion). Passed through to model_setup, which fits and stores the glmmTMB reference (design$glmmTMB_fit) reused below as fit_ref. Gaussian models only compute the per-group group.ing_prior calibration (a within-group regression fit for every group level, used only by dGamma_list()) when dispformula is the group formula; dispformula = ~1 skips it entirely, so group.dispersion.pwt / group.dispersion.nprior must then be scalar (pooled). Mirrors the dispformula argument on lmerb()/glmerb() in lmebayes, which gates the analogous choice of sampler route; the two are independent arguments that must be kept consistent by the caller.
pop.pwt	Prior weight(s) in $(0, 1)$ for the population-level (Block~2, fixef) coefficients. Either a scalar (applied to every random-effect component and every Block~2 predictor), or a list with one element per random-effect component (named with the RE coefficient names in any order, or unnamed positional). Each list element is a scalar (recycled over that component's Block~2 predictors) or a vector of length $p_k$ (optionally named with the predictor column names of X_hyper[[k]], reordered to match). The prior covariance for each fixef_k block is scaled relative to vcov(fit_ref) following the Prior_Setup convention: $(1-\mathrm{pop.pwt})/\mathrm{pop.pwt}$ for a scalar, and elementwise $\sqrt{(1-\mathrm{pop.pwt}_i)/\mathrm{pop.pwt}_i}\, \sqrt{(1-\mathrm{pop.pwt}_j)/\mathrm{pop.pwt}_j}$ for vectors. At most one of pop.pwt (explicitly supplied), pop.nprior, and pop.sd may be supplied. When pop.pwt is NULL (the default), dimension-adaptive defaults are chosen separately for each random-effect component (see pop.pwt_default_low/ pop.pwt_default_high).
pop.pwt_default_low	Scalar in $(0, 1)$, default 0.01. Used when pop.pwt is NULL and a random-effect component's Block~2 predictor count $p_k = ncol(X\_hyper[[k]])$ is strictly less than 14 (mirrors Prior_Setup pwt_default_low).
pop.pwt_default_high	Scalar in $(0, 1)$, default 0.05. Used when pop.pwt is NULL and $p_k \ge 14$ for that component (mirrors Prior_Setup pwt_default_high).
pop.nprior	Optional absolute effective prior sample size(s), in group units ($J =$nlevels(group)), for the population-level (Block~2, fixef) coefficients – an alternative to pop.pwt, on the same $n_i = J\,w_i/(1-w_i)$ duality as pop.dispersion.pwt/pop.dispersion.nprior. Same shape as pop.pwt: a scalar, or a list with one element per random-effect component (named or unnamed positional), each element NULL (keep the pop.pwt-derived default for that component), a scalar, or a vector of length $p_k$ (optionally named with the predictor column names of X_hyper[[k]], reordered to match). Converted internally to an effective weight $w_i = \mathrm{pop.nprior}_i / (\mathrm{pop.nprior}_i + J)$, used in place of pop.pwt for that component's Sigma scaling. At most one of pop.pwt (explicitly supplied), pop.nprior, and pop.sd may be supplied. The returned pop.nprior field (see Value below) is always populated with the effective prior sample size actually implied (derived from whichever of pop.pwt/ pop.sd was used when pop.nprior is not supplied), in the same shape this argument expects – call once without pop.nprior to inspect the defaults, then feed a modified copy of ps$pop.nprior back in to override.
pop.sd	Optional prior standard deviation(s), on the coefficient scale of vcov(fit_ref), for the population-level (Block~2, fixef) coefficients – an alternative to pop.pwt. Same shape as pop.pwt: a scalar, or a list with one element per random-effect component (named or unnamed positional), each element NULL (keep the pop.pwt-derived default for that component – a partial override across components), a scalar, or a vector of length $p_k$ (optionally named with the predictor column names of X_hyper[[k]], reordered to match). Converted internally to an effective weight $w_i = V_{ii} / (V_{ii} + \mathrm{pop.sd}_i^2)$, where $V_{ii}$ is the corresponding diagonal entry of vcov(fit_ref), and used in place of pop.pwt for that component's Sigma scaling. At most one of pop.pwt (explicitly supplied), pop.nprior, and pop.sd may be supplied. The returned pop.sd field (see Value below) is always populated with the effective standard deviations actually used ($\sqrt{\mathrm{diag}(\Sigma_{\mathrm{fixef},k})}$, derived from whichever of pop.pwt/pop.nprior was used when pop.sd is not supplied), in the same shape this argument expects – call once without pop.sd to inspect the defaults, then feed a modified copy of ps$pop.sd back in to override.
pop.dispersion.pwt	Optional relative prior weight(s) in $(0, 1)$ for the Block~2 dispersion (precision) prior, decoupled from pop.pwt. A scalar, or a list / numeric vector with one value per random-effect component (named or positional). Converted internally to an effective prior sample size $n_k = J\,w_k/(1-w_k)$ where $J$ is the number of groups. At most one of pop.dispersion.pwt and pop.dispersion.nprior may be supplied; when neither is, the value is derived from pop.pwt (the mean across a component's predictors), keeping the Block~2 $\tau^2_k$ dispersion prior consistent with the coefficient prior strength. Weak values carry no computational penalty for dIndependent_Normal_Gamma sampling: the $\tau^2$ truncation window comes from limiting-posterior quantiles independent of the prior strength (see pop.ing_prior below).
pop.dispersion.nprior	Optional absolute effective prior sample size(s) (in group units) for the Block~2 $\tau^2_k$ dispersion prior. A positive scalar, or a list / numeric vector with one value per random-effect component (named or positional). See pop.dispersion.pwt.
pop.max_disp_perc	Scalar in $(0.5, 1)$, default 0.99, or a named/positional length-$p_{\mathrm{re}}$ vector (one value per random-effect component). Same tail-probability role as group.max_disp_perc, but for the Block~2 $\tau^2_k$ truncation window stored in pop.ing_prior. Only operationally consumed when the sampler's dispersion pfamily is dIndependent_Normal_Gamma (ignored for dNormal).
pop.intercept_source	Character string controlling the prior mean for the global intercept hyperparameter (Intercept)::(Intercept) only. One of "null_model" (default) or "full_model". When "null_model", the prior mean is taken from a random-intercept-only reference fit y ~ 1 + (1 | group) that omits all fixed-effect predictors (analogous to Prior_Setup with intercept_source = "null_model"). When "full_model", the full-model MLE intercept is used.
pop.effects_source	Character string controlling the prior mean for all other Block~2 hyperparameters (including population-mean slopes stored as (Intercept) columns in non-intercept RE components, and any non-intercept columns in X_hyper). One of "null_effects" (default) or "full_model". When "null_effects", prior means are set to zero. When "full_model", full-model MLE values are used.
pop.mu	Optional prior mean override(s) for the population-level (Block~2, fixef) coefficients – a finer-grained alternative to pop.intercept_source/pop.effects_source; prefer those first, and reach for pop.mu only when a specific mean needs to differ from either derived default. Same shape as pop.pwt: a scalar, or a list with one element per random-effect component (named or unnamed positional), each element NULL (keep the pop.intercept_source/pop.effects_source-derived default mean for that component – a partial override across components), a scalar, or a vector of length $p_k$ (optionally named with the predictor column names of X_hyper[[k]], reordered to match). Replaces mu wholesale for the components supplied; components left NULL keep their usual pop.intercept_source/pop.effects_source default. The returned pop.mu field (see Value below) is always populated with the values actually used (derived defaults when pop.mu is not supplied), in the same shape this argument expects – call once without pop.mu to inspect the defaults, then feed a modified copy of ps$pop.mu back in to override.
group.dispersion.pwt	Optional relative prior weight(s) in $(0, 1)$ for the Block~1 observation $\sigma^2$ Gamma prior (Gaussian models only), decoupled from pop.pwt (Block~2 fixef) and from pop.dispersion.pwt (Block~2 $\tau^2_k$). Either a scalar converted to $n_{\mathrm{prior}} = w/(1-w)\times n$ on the total observation count $n$ (pooled group.ing_prior), or a named or positional vector of length $J$ ($J =$ nlevels(groups)) with one weight per group level for per-group dGamma_list() calibration ($n_{\mathrm{prior},j} = w_j/(1-w_j)\times n_j$). When group.dispersion.pwt is a vector, the pooled group.ing_prior continues to use the default w = 0.01 on total $n$ (unless group.dispersion.nprior is supplied explicitly). At most one of group.dispersion.pwt and group.dispersion.nprior may be supplied; when neither is, scalar group.dispersion.pwt = 0.01 is used for both pooled and per-group paths.
group.dispersion.nprior	Optional positive scalar: absolute effective prior sample size for the Block~1 $\sigma^2$ prior (observation units, not groups). See group.dispersion.pwt.
group.alpha_target	Scalar in $(0, 1)$, default 0.01, or NULL to disable. Gaussian models with per-group dispformula only: target percentage of posterior beta_j draws allowed to fall outside the per-group Mahalanobis/log-concavity ellipsoid of inst/BLOCK_GIBBS_ERGODICITY_ING.md Section 16.6. This is a different quantity from group.dispersion.pwt: group.dispersion.pwt is the Gamma prior weight (an input to the prior), while group.alpha_target is a target violation rate (an output being aimed at). When non-NULL and dispformula requests per-group dispersion, Prior_Setup_GLMM() searches, for each group, for the smallest group.dispersion.pwt driving the predicted violation rate down to group.alpha_target, floored at the group.dispersion.pwt resolved from the group.dispersion.pwt/group.dispersion.nprior arguments above (so calibration only ever sharpens, never loosens, the per-group prior; supplying group.dispersion.pwt explicitly still guarantees a minimum). Active by default as a safeguard against poorly-behaved per-group priors; set to NULL to opt out and use group.dispersion.pwt/ group.dispersion.nprior directly, unmodified. Silently ignored (not an error) when left at its default on a model that does not support per-group calibration (e.g.\ dispformula = ~1); it is only an error to explicitly supply a non-NULL value in that case. See group.pwt_calibration below.
group.max_disp_perc	Scalar in $(0.5, 1)$, default 0.99, or a named/positional length-$J$ vector (one value per group level, Gaussian models with per-group dispformula only). Tail probability used to compute the Block~1 $\sigma^2$ truncation window stored in group.ing_prior – the pooled scalar default when a vector is supplied here (dispformula = ~1), or the resolved per-group value (per group, per-group dispformula). The window is the central $2 \times \mathrm{group.max\_disp\_perc} - 1$ mass interval: disp_lower = $(1-\mathrm{group.max\_disp\_perc})$ quantile and disp_upper = $\mathrm{group.max\_disp\_perc}$ quantile of the relevant Gamma (precision) distribution, inverted to the dispersion scale. Tighter values (e.g.\ 0.95) shrink the truncation window and typically improve Block~1 acceptance rates at the cost of slightly less envelope coverage; looser values (e.g.\ 0.999) widen the window. Passed to dGamma() as max_disp_perc. A per-group vector lets outlier groups be given a tighter (or looser) window than the rest; see dGamma_list().

Details

Block 1 (per-group, independent):

$p(\mathbf{b}_j \mid \mathbf{y}, \mathrm{fixef}, \sigma^2, \Sigma_b) = \mathcal{N}(\boldsymbol{\mu}_{b,j}^*, \boldsymbol{\Sigma}_{b,j}^*)$

$\boldsymbol{\Sigma}_{b,j}^{*-1} = \mathbf{Z}_j'\mathbf{Z}_j / \sigma^2 + \mathrm{diag}(1/\tau^2_k)$

when family = gaussian(). For non-Gaussian families there is no observation-level dispersion; Block~1 uses dNormal with ddef = TRUE (see dNormal).

Block 2 (per-RE coefficient $k$, independent):

$p(\mathrm{fixef}_k \mid \mathbf{b}_k, \tau^2_k) = \mathcal{N}(\boldsymbol{\mu}_{\mathrm{fixef},k}^*, \boldsymbol{\Sigma}_{\mathrm{fixef},k}^*)$

$\boldsymbol{\Sigma}_{\mathrm{fixef},k}^{*-1} = \mathbf{X}_k'\mathbf{X}_k / \tau^2_k + \boldsymbol{\Sigma}_{\mathrm{fixef},k}^{-1}$

Naming convention. Top-level fields use pop.* for Block~2 (population / $\gamma$, $\tau^2_k$) and group.* for Block~1 (within-group / observation $\sigma^2$, RE covariance $\Psi$). Nested entries in pop.prior_list[[k]] use glmbayes-style mu, Sigma, and dispersion (where dispersion is the plug-in RE variance $\tau^2_k$).

Why default calibration depends on classical estimates. Prior_Setup_GLMM scales Block~2 covariances from vcov(fit_ref) by $(1-\mathrm{pwt})/\mathrm{pwt}$ and plugs in RE variances from the full reference fit, where fit_ref is the pooled lmer/glmer fit when dispformula = ~1, or an equivalent glmmTMB fit with the same dispformula when per-group dispersion is requested (see calibration_source and mer_fit above). By default the global intercept prior mean comes from a random-intercept-only null fit; all other prior means are zero (pop.effects_source = "null_effects"). This requires:

1. Converged reference lmer/glmer fit from model_setup on the full formula (and a random-intercept-only null fit when pop.intercept_source = "null_model"). Fits with lme4 checkConv failures (e.g.\ large max|grad|) are rejected.

2. Every X_hyper[[k]] column maps to a fixef(fit_ref) term.

3. Each RE variance $\tau^2_k$ from the reference fit is strictly positive.

Value

Object of class "Prior_Setup_GLMM" with fields:

formula

Model formula.

family

Family object.

pop.pwt

Prior weight(s) used: the scalar as supplied, the canonical named list of per-predictor weight vectors when a list was supplied, or the dimension-adaptive default list (one weight per Block~2 predictor, per RE component) when pop.pwt was NULL.

pop.nprior

Named list with one named numeric vector per random-effect component – the effective prior sample size (group units) actually implied, $n_i = J\,w_i/(1-w_i)$, whether derived from pop.pwt/pop.sd (the usual case) or overridden via pop.nprior. Always present, and in the same shape the pop.nprior argument expects, so it can be inspected from a first call and fed back in (after editing) as pop.nprior to a second call.

pop.sd

Named list with one named numeric vector per random-effect component – the effective coefficient-scale prior standard deviations actually used ($\sqrt{\mathrm{diag}(\Sigma_{\mathrm{fixef},k})}$), whether derived from pop.pwt/pop.nprior (the usual case) or overridden via pop.sd. Always present, and in the same shape the pop.sd argument expects, so it can be inspected from a first call and fed back in (after editing) as pop.sd to a second call.

pop.mu

Named list with one named numeric vector per random-effect component – the prior means actually used (mu), whether derived from pop.intercept_source/pop.effects_source (the usual case) or overridden via pop.mu. Always present, and in the same shape the pop.mu argument expects, so it can be inspected from a first call and fed back in (after editing) as pop.mu to a second call.

pop.dispersion.pwt

Named per-component vector of relative dispersion prior weights (always present; consistent with pop.dispersion.nprior via $w_k = n_k/(n_k + J)$).

pop.dispersion.nprior

Named per-component vector of effective prior sample sizes for the Block~2 dispersion prior (always present; used by pfamily_list() to calibrate dIndependent_Normal_Gamma components).

group.dispersion.pwt

Gaussian models only: scalar or length-$J$ vector of relative prior weights for Block~1 $\sigma^2$ – a scalar unless a per-group vector was supplied or group.alpha_target calibration was active, in which case it is the resolved length-$J$ vector (see dGamma_list()).

group.dispersion.nprior

Gaussian models only: effective prior sample size for Block~1 $\sigma^2$, on the same pooled-vs.-vector basis as group.dispersion.pwt (a scalar on the observation scale, or a named length-$J$ vector of per-group $n_{\mathrm{prior},j}$ reflecting the group.alpha_target-calibrated group.dispersion.pwt when calibration was active; see group.pwt_calibration).

group.alpha_target

The value used, or NULL if calibration was not active.

group.pwt_calibration

NULL unless calibration was active; otherwise a per-group diagnostics data.frame (group, n_j, pwt_floor, pct_outside_before, w_star, clipped_at_ceiling, w_final, floor_binds, pct_outside_after) from the search described under group.alpha_target.

group.max_disp_perc

Gaussian models only: scalar or length-$J$ vector of the tail probability used for the Block~1 $\sigma^2$ truncation window (group.ing_prior) – the scalar as supplied, or the resolved length-$J$ vector when a vector was supplied.

pop.max_disp_perc

Named length-$p_{\mathrm{re}}$ vector of resolved per-component values used for the Block~2 $\tau^2_k$ truncation windows in pop.ing_prior (always present).

pop.intercept_source

How the global intercept prior mean was chosen ("null_model" or "full_model").

pop.effects_source

How other Block~2 prior means were chosen ("null_effects" or "full_model").

block_formula

Within-group Block~1 formula: response regressed on the random-coefficient structure only (columns of design$groupef.names); level-2 hyper covariates are excluded. Used by dGamma_list().

group.tau_sd

Named vector sqrt(Psi) from the reference fit; shared population RE standard deviations for per-group calibration.

data

Data frame passed to Prior_Setup_GLMM() (reference for dGamma_list() diagnostics).

design

Full model_setup object (all groups), including length-n weights and offset vectors aligned with design$y.

mer_fit

Reference lmer/glmer fit on all groups (the full-formula fit from model_setup), always present regardless of dispformula (backs group.dispersion and x$lmer/x$glmer in lmerb()/glmerb()).

fit_ref

The calibration reference for Block~2 (fixef/$\tau^2_k$) and the per-group Block~1 inputs (group.tau_sd, BLUP coefficients): identical to mer_fit when dispformula = ~1; otherwise an equivalent glmmTMB fit with the same dispformula (Gaussian models only). See calibration_source.

group.dispersion.fit

NULL when dispformula = ~1; otherwise the same glmmTMB object as fit_ref.

group.dispersion.ref

NULL when dispformula = ~1; otherwise a named length-$J$ vector of per-group observation-level dispersion read from fit_ref's dispersion linear predictor (glmmTMB::predict(fit_ref, type = "disp"), aggregated by group). Diagnostic only – not the value fed to the sampler; compare against group.ing_prior's sigma2_hat.

calibration_source

"lme4" or "glmmTMB"; which reference fit produced fit_ref (and therefore pop.prior_list, pop.ing_prior, group.tau_sd, and group.ing_prior).

group.dispersion

For Gaussian models: resolved Block~1 observation $\sigma^2$ used to mean-match the Gamma prior. When dispformula = ~1, the pooled A12 §3.3.4 aggregate $\hat\sigma^2 = (\sum_j S_{\mathrm{marg},j})/(n - J p_{\mathrm{re}})$ (Part VI Omega_j included); classical lmer residual in attr(., "classical"). When dispformula is per-group, a named length-$J$ vector of per-group group.ing_prior[[j]]$sigma2_hat. attr(., "source") records provenance. NULL for non-Gaussian families. Pass to rlmerb()/rglmerb() as group.dispersion for fixed-$\sigma^2$ routes.

group.Sigma

Diagonal RE covariance matrix (Block~1), classical reference-fit $\Psi$ diagonals.

pop.dispersion.ref

Named vector of classical reference-fit $\hat\tau^2_k$ (VarCorr / Psi) before the A12 §3.3.4 hyper-regression remargin; compare to pop.prior_list[[k]]$dispersion.

pop.prior_list

Named Block~2 prior list per RE coefficient. For Gaussian models dispersion is the A12 §3.3.4 marginal hyper-regression center for $\tau^2_k$ (not the raw VarCorr plug-in; see pop.dispersion.ref).

pop.ing_prior

Named per-component list of the prospective dIndependent_Normal_Gamma calibration: Gamma precision-prior shape $= (n_0 + 1 + p_k)/2$ and rate $= \hat\tau^2_k (n_0 + p_k - 1)/2$ (the lmebayesCore default calibration with $n_0 =$ pop.dispersion.nprior; since rate $= \hat\tau^2_k (\code{shape} - 1)$, the implied inverse-Gamma prior on $\tau^2_k$ has mean exactly $\hat\tau^2_k$), where $\hat\tau^2_k$ is the A12 §3.3.4 marginal hyper-regression center in pop.prior_list[[k]]$dispersion, and the default $\tau^2_k$ truncation window disp_lower / disp_upper: the $(1-\mathrm{pop.max\_disp\_perc})$ / $\mathrm{pop.max\_disp\_perc}$ quantiles of the limiting posterior $\Gamma((J+1)/2,\; \hat\tau^2_k (J-1)/2)$ – the weak-prior ($n_0 \to 0$) limit of the Block~2 posterior Gamma for the precision (lmebayesCore Chapter A12, Theorem 2; inverted to a $\tau^2$ interval). This window is identical for all $n_0$, covers $\ge$ $2 \times \mathrm{pop.max\_disp\_perc} - 1$ of the exact posterior for every prior strength, and keeps the envelope sampler's cost stable as priors weaken; see inst/ING_TRUNCATION_WINDOW.md. Used by pfamily_list() when ptypes = "dIndependent_Normal_Gamma"; ignored for dNormal priors.

group.ing_prior

Gaussian models only; NULL for non-Gaussian models. Shape depends on dispformula:

dispformula = ~1 (pooled)

A single list: prospective dGamma() group.dispersion calibration for Block~1 ING (observation $\sigma^2$ shared across all group levels) – mean-matched shape / rate with $n_{\mathrm{prior}} = \mathrm{group.dispersion.pwt}/(1-\mathrm{group.dispersion.pwt})\times n$, $p = p_{\mathrm{re}}$, and $\hat\sigma^2$ = group.dispersion (pooled A12 §3.3.4 aggregate of per-group $S_{\mathrm{marg},j}$ by default; same ING algebra as pop.ing_prior for $\tau^2_k$; requires group.dispersion.pwt $\le 0.5$), plus disp_lower / disp_upper as the central $2 \times \mathrm{group.max\_disp\_perc} - 1$ prior-mass interval from the same calibrated shape/rate (always the pooled/scalar group.max_disp_perc value, even when a per-group vector was supplied). Pass these fields to dGamma().

dispformula requesting per-group dispersion

A named list (one entry per group level) of per-group dGamma() density calibration (sigma2_hat, shape_ING, rate, rate_gamma, n_prior, n_j, n_combined, disp_lower, disp_upper, max_disp_perc, omega_j, ...). shape_ING/rate already fold in the Part VI model-derived Omega_j (fixed-effect/gamma uncertainty about b_j's prior mean, propagated through group j's hyper-design row and pop.prior_list[[k]]$Sigma), so this marginal integrates out both random effects (b_j) and fixed effects (gamma); see inst/DGAMMA_LIST_MARGINAL_AND_BOUNDS.md Part VI. disp_lower/disp_upper are $(1-w_j)$/$w_j$ quantiles of Gamma(shape_ING + n_j/2, rate_post), where rate_post is sigma2_hat mean-matched at that same inflated shape ($n_j/2$ added only for the window, so it reflects the posterior spread the sampler's own envelope machinery actually draws $\sigma^2_j$ from – see EnvelopeDispersionBuild.cpp's shape2 = Shape + n_w/2 fallback – while shape_ING/rate themselves, the values actually fed to the sampler as the Gamma prior, are unchanged); $w_j$ is this group's own resolved group.max_disp_perc value (stored per-group as max_disp_perc). Used directly by dGamma_list().

Which shape is returned is fully determined by dispformula; see dGamma_list() for how consumers tell the two apart.

dispformula

The dispformula supplied.

x invisibly.

See Also

model_setup, Prior_Setup, build_mu_all, pfamily_list, dGamma_list

---