| Title: | Bayesian Linear Mixed-Effects Models via Two-Block Gibbs Sampling |
|---|---|
| Description: | Provides Bayesian linear and generalized linear mixed-effects model fitting with near-independent posterior samples. The main functions mirror the 'lme4' lmer() and glmer() interfaces (X/beta for fixed effects, Z/b for random effects, D for the random-effects covariance) while adding prior family specifications for Gaussian, Poisson, binomial, and Gamma models with log-concave likelihoods. Sampling uses efficient two-block Gibbs sampling procedures with parallel chains and theoretically driven stopping rules. Gibbs sampling engines are provided by the 'lmebayesCore' package. |
| Authors: | Kjell Nygren [aut, cre], The R Core Team [ctb, cph] (R Mathlib sources, R stats modeling code, and derived/adapted routines), The R Foundation [cph] (Portions of R Mathlib and R source code), Ross Ihaka [ctb, cph] (R Mathlib and original R modeling infrastructure), Robert Gentleman [ctb, cph] (Portions of R Mathlib), Simon Davies [ctb] (Original R glm implementation), Morten Welinder [ctb, cph] (Portions of R Mathlib), Martin Maechler [ctb] (Portions of R Mathlib) |
| Maintainer: | Kjell Nygren <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.1.0 |
| Built: | 2026-10-02 21:45:27 UTC |
| Source: | https://github.com/knygren/lmebayes |
Two-block Gibbs samplers for Bayesian linear and generalized linear mixed-effects models, following lme4 notation. Builds on glmbayes for iid GLM sampling within blocks.
Row-block interfaces include lmbBlock and glmbBlock;
mixed-model setup from lme4 formulas via model_setup.
Lower-level simulation uses simfunction and envelope
utilities from lmebayesCore.
See the package README at https://github.com/knygren/lmebayes for examples.
The vignettes are a curriculum, meant to be read roughly in order. Each assumes the previous ones.
vignette("Chapter-00")Roadmap: what the package does and which chapter answers which question.
vignette("Chapter-01")Getting started: one model end to
end on lme4::sleepstudy.
vignette("Chapter-02")From lm/glm to mixed
models: pooling, the centered parameterization, and the model-class
restrictions.
vignette("Chapter-03")The prior workflow:
model_setup, Prior_Setup_GLMM, pfamily_list.
vignette("Chapter-04")Linear mixed models with group-level covariates and cross-level moderation.
vignette("Chapter-05")Reading the output: every accessor and container.
vignette("Chapter-06")Convergence and tv_tol.
vignette("Chapter-07")Estimated variance components and the pilot stage.
vignette("Chapter-08")Observation dispersion, including
per-group .
vignette("Chapter-09")Poisson GLMMs.
vignette("Chapter-10")Binomial GLMMs.
vignette("Chapter-11")Point estimates without simulation.
vignette("Chapter-12")Row-block models.
vignette("Chapter-13")Comparison with lme4 and glmmTMB.
Engine internals, derivations, and convergence theory are documented in
the Chapter-B* vignettes of lmebayesCore; the iid GLM stack
underneath is documented in glmbayes.
In interactive sessions, attaching the package with library(lmebayes)
may emit a short packageStartupMessage
when has_opencl() is FALSE (typical for CRAN binaries) but a
GPU or OpenCL stack appears available on the host. OpenCL modelling paths
require a source install with OpenCL at compile time;
has_opencl() then reports whether that build succeeded.
Set options(glmbayes.quiet_opencl_startup = TRUE) to suppress attach
notes (recommended for CI and R CMD check).
Kjell Nygren
There are no references for Rd macro \insertAllCites on this help page.
lmerb, model_setup, lmbBlock, glmbBlock;
simfunction, EnvelopeBuild;
lmb and glmb for fixed-effects-only Bayesian
linear and generalized linear models;
rlmb and rglmb for iid
posterior draws.
Useful links:
set.seed(333) ## Dobson (1990) Page 93: Randomized Controlled Trial : counts <- c(18, 17, 15, 20, 10, 20, 25, 13, 12) outcome <- gl(3, 1, 9) treatment <- gl(3, 3) print(d.AD <- data.frame(treatment, outcome, counts)) glm.D93 <- glm(counts ~ outcome + treatment, family = poisson()) ## Plain iid GLM fitting engine (rglmb) is not re-exported from lmebayes; ## it lives in glmbayesCore, so we source it from there directly. ps <- Prior_Setup(counts ~ outcome + treatment, family = poisson()) rglmb.D93 <- glmbayesCore::rglmb( n = 200, y = ps$y, x = as.matrix(ps$x), pfamily = dNormal(mu = ps$mu, Sigma = ps$Sigma), family = poisson(), weights = rep(1, nrow(ps$x)) ) print(rglmb.D93) summary(rglmb.D93)set.seed(333) ## Dobson (1990) Page 93: Randomized Controlled Trial : counts <- c(18, 17, 15, 20, 10, 20, 25, 13, 12) outcome <- gl(3, 1, 9) treatment <- gl(3, 3) print(d.AD <- data.frame(treatment, outcome, counts)) glm.D93 <- glm(counts ~ outcome + treatment, family = poisson()) ## Plain iid GLM fitting engine (rglmb) is not re-exported from lmebayes; ## it lives in glmbayesCore, so we source it from there directly. ps <- Prior_Setup(counts ~ outcome + treatment, family = poisson()) rglmb.D93 <- glmbayesCore::rglmb( n = 200, y = ps$y, x = as.matrix(ps$x), pfamily = dNormal(mu = ps$mu, Sigma = ps$Sigma), family = poisson(), weights = rep(1, nrow(ps$x)) ) print(rglmb.D93) summary(rglmb.D93)
Data with information on 17 overdoses of the drug amitriptyline
data(AMI)data(AMI)
This data frame contains the following columns:
TOTtotal TCAD plasma level
AMIamount of amitriptyline present in the TCAD plasma level
GENgender (male = 0, female = 1)
AMTamount of drug taken at time of overdose
PRPR wave measurement
DIAPdiastolic blood pressure
QRSQRS wave measurement
Each row is one overdose episode. Variables include total tricyclic antidepressant level, amitriptyline component, gender, reported amount ingested, and ECG-related measures (PR interval, QRS duration, diastolic blood pressure). The dataset is used in package examples for binomial and related regression; see (Dobson 1990) for analogous generalized linear modelling of clinical outcomes.
Dobson A~J (1990). An Introduction to Generalized Linear Models. Chapman and Hall, London.
############################### Start of AMI dataset example #################### data(AMI) summary(AMI) ############################################################################### ## End of AMI dataset example ############################################################################################################## Start of AMI dataset example #################### data(AMI) summary(AMI) ############################################################################### ## End of AMI dataset example ###############################################################################
ranef.lmerb objects to long data.frame formatCoerce ranef.lmerb objects to long data.frame format
## S3 method for class 'ranef.lmerb' as.data.frame(x, ...)## S3 method for class 'ranef.lmerb' as.data.frame(x, ...)
x |
An object of class |
... |
Ignored. |
A data.frame with columns grpvar, term,
grp, condval, and condsd (when postVar is
available on x).
A processed version of the UCI Bike Sharing Dataset (hourly data). The data include derived variables for part of day, quarter, and Fourier terms for cyclic effects of hour and month.
BikeSharingBikeSharing
A data frame with 17,379 observations and 25 variables (original plus derived):
Record index.
Date.
Season (1: spring, 2: summer, 3: fall, 4: winter).
Year (0: 2011, 1: 2012).
Month (1–12).
Hour (0–23).
Whether the day is a holiday (0/1).
Day of week (0–6).
Working day (0/1).
Weather situation (1–4).
Normalized temperature.
Normalized feeling temperature.
Normalized humidity.
Normalized wind speed.
Count of casual users.
Count of registered users.
Total count (casual + registered).
Hour as numeric (same as hr).
Month as integer.
Factor: Night (0–5h), Morning (6–11h), Afternoon (12–17h), Evening (18–23h).
Factor: Q1–Q4.
Sine term for 24-hour cycle.
Cosine term for 24-hour cycle.
Sine term for 12-month cycle.
Cosine term for 12-month cycle.
UCI Machine Learning Repository: Bike Sharing Dataset. https://archive.ics.uci.edu/dataset/275/bike+sharing+dataset
############################### Start of BikeSharing dataset example #################### data("BikeSharing") head(BikeSharing) dim(BikeSharing) ############################################################################### ## End of BikeSharing dataset example ############################################################################################################## Start of BikeSharing dataset example #################### data("BikeSharing") head(BikeSharing) dim(BikeSharing) ############################################################################### ## End of BikeSharing dataset example ###############################################################################
A copy of Boston where all predictors (every column except
medv) have been mean-centered (subtract column means, no scaling).
data("Boston_centered")data("Boston_centered")
A data frame with 506 observations and 14 variables (same names as
Boston). See ?MASS::Boston for variable
descriptions.
Derived from MASS::Boston. Original data described in Harrison
and Rubinfeld (1978); see ?Boston in MASS.
############################### Boston_centered dataset example #################### data("Boston_centered") head(Boston_centered) summary(Boston_centered) ## Predictors are mean-centered (column means ~0) predictors <- setdiff(names(Boston_centered), "medv") colMeans(Boston_centered[predictors]) form <- medv ~ crim + zn + indus + chas + nox + age + dis + rad + tax + ptratio + black + lstat + rm lm.boston <- lm(form, data = Boston_centered, x = TRUE, y = TRUE) summary(lm.boston) ############################################################################### ## End of Boston_centered dataset example ############################################################################################################## Boston_centered dataset example #################### data("Boston_centered") head(Boston_centered) summary(Boston_centered) ## Predictors are mean-centered (column means ~0) predictors <- setdiff(names(Boston_centered), "medv") colMeans(Boston_centered[predictors]) form <- medv ~ crim + zn + indus + chas + nox + age + dis + rad + tax + ptratio + black + lstat + rm lm.boston <- lm(form, data = Boston_centered, x = TRUE, y = TRUE) summary(lm.boston) ############################################################################### ## End of Boston_centered dataset example ###############################################################################
The data give the Canadian automobile insurance experience for policy years 1956 and 1957 as of June 30, 1959. The data includes virtually every insurance company operating in Canada and was collated by the Statistical Agency (Canadian Underwriters' Association - Statistical Department) acting under instructions from the Superintendent of Insurance. The data given here is for private passenger automobile liability for non-farmers for all of Canada excluding Saskatchewan.
data(carinsca)data(carinsca)
A data frame with 20 observations on the following 6 variables:
MeritMerit Rating:
3 - licensed and accident free 3 or more years
2 - licensed and accident free 2 years
1 - licensed and accident free 1 year
0 - all others
Class1 - pleasure, no male operator under 25
2 - pleasure, non-principal male operator under 25
3 - business use
4 - unmarried owner or principal operator under 25
5 - married owner or principal operator under 25
InsuredEarned car years
PremiumEarned premium in 1000's
(adjusted to what the premium would have been had all cars been written at 01 rates)
ClaimsNumber of claims
CostTotal cost of the claim in 1000's of dollars
One could apply Poisson regression to the number of claims and gamma regression to the cost per claim.
Bailey, R. A., and Simon, LeRoy J. (1960). Two studies in automobile insurance ratemaking. ASTIN Bulletin, 192-217.
Data downloaded from http://www.statsci.org/data/general/carinsca.html. That site also contains classical Poisson and Gamma regression examples.
############################### Start of carinsca dataset example #################### data(carinsca) str(carinsca) head(carinsca) ############################################################################### ## End of carinsca dataset example ############################################################################################################## Start of carinsca dataset example #################### data(carinsca) str(carinsca) head(carinsca) ############################################################################### ## End of carinsca dataset example ###############################################################################
A cleaned version of the Cleveland heart disease dataset from the UCI Machine
Learning Repository. This version contains only complete cases and includes a
derived binary outcome variable hd indicating the presence ("Yes")
or absence ("No") of heart disease.
data("Cleveland")data("Cleveland")
A data frame with 297 observations and 15 variables:
Age in years (numeric).
Sex (0 = female, 1 = male).
Chest pain type (numeric code 1–4).
Resting blood pressure (mm Hg).
Serum cholesterol (mg/dl).
Fasting blood sugar > 120 mg/dl (1 = true, 0 = false).
Resting electrocardiographic results (numeric code).
Maximum heart rate achieved.
Exercise-induced angina (1 = yes, 0 = no).
ST depression induced by exercise relative to rest.
Slope of the peak exercise ST segment.
Number of major vessels colored by fluoroscopy (0–3).
Thalassemia status (numeric code).
Original UCI disease score (0–4).
Binary heart disease indicator: "No" (num = 0) or "Yes" (num > 0).
UCI Machine Learning Repository: Heart Disease Data Set. https://archive.ics.uci.edu/dataset/45/heart+disease
############################### Start of Cleveland dataset example #################### data("Cleveland") head(Cleveland) summary(Cleveland) ############################################################################### ## End of Cleveland dataset example ############################################################################################################## Start of Cleveland dataset example #################### data("Cleveland") head(Cleveland) summary(Cleveland) ############################################################################### ## End of Cleveland dataset example ###############################################################################
Mirrors coef.merMod: Block~1 stored coefficients
(ranef.lmerb) are already on the per-group coefficient scale used
internally by the sampler (not lme4-style deviations).
## S3 method for class 'lmerb' coef(object, type = c("mean", "mode"), ...) ## S3 method for class 'glmerb' coef(object, type = c("mean", "mode"), ...)## S3 method for class 'lmerb' coef(object, type = c("mean", "mode"), ...) ## S3 method for class 'glmerb' coef(object, type = c("mean", "mode"), ...)
object |
An |
type |
|
... |
Ignored. |
A named list (one element per grouping factor) of data.frames
with the same layout as coef(merMod).
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
priors for lmerb(..., group.dispersion = ...).
dGamma_list(object, ...)dGamma_list(object, ...)
object |
A prior-specification object. |
... |
Additional arguments passed to methods (e.g.
|
Where pfamily_list configures the between-group
variance components , dGamma_list() configures the
within-group observation variance. It returns one
dGamma prior per group level, so each group
gets its own rather than sharing a pooled value.
Each element places a Gamma prior on the precision
(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.
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.
dGamma, Prior_Setup_GLMM,
pfamily_list, rLMM_reg
## dGamma_list() — per-group dGamma from Prior_Setup_GLMM (Ex_13b data) ## ## Requires dispformula = ~school_id so group.ing_prior is per-group. ## Requires suggested package bayesrules. if (requireNamespace("bayesrules", quietly = TRUE)) { data(big_word_club, package = "bayesrules") dat <- big_word_club dat$school_id <- factor(dat$school_id) dat <- subset( dat, !is.na(score_ppvt) & !is.na(invalid_ppvt) & invalid_ppvt == 0L & complete.cases(dat[, c( "score_ppvt", "distracted_a1", "distracted_ppvt", "private_school", "title1", "free_reduced_lunch", "school_id" )]) ) form_lmer <- score_ppvt ~ private_school + title1 + free_reduced_lunch + distracted_ppvt + distracted_a1 + free_reduced_lunch:distracted_a1 + (1 + distracted_ppvt + distracted_a1 || school_id) design0 <- model_setup(form_lmer, data = dat) dat <- subset(dat, school_id %in% names(design0$groupef.rank)[design0$groupef.rank]) dat$school_id <- droplevels(dat$school_id) ps <- Prior_Setup_GLMM( form_lmer, data = dat, pop.pwt = 0.01, dispformula = ~school_id, group.max_disp_perc = 0.8, group.dispersion.pwt = 0.1, group.alpha_target = NULL ) disp_pf <- dGamma_list(ps) ## Group names (school_id levels) g <- names(disp_pf) utils::head(g, 10L) ## Print a subset (recommended when J is large) print(disp_pf, groups = c("2", "18", "47")) ## Variation: rebuild truncation windows; print one group disp_pf_tight <- dGamma_list(ps, max_disp_perc_measurement = 0.9) print(disp_pf_tight, groups = "2") ## Downstream attributes (not shown by print.dGamma_list) names(attributes(disp_pf)) utils::head(attr(disp_pf, "window_diagnostics"), 3L) }## dGamma_list() — per-group dGamma from Prior_Setup_GLMM (Ex_13b data) ## ## Requires dispformula = ~school_id so group.ing_prior is per-group. ## Requires suggested package bayesrules. if (requireNamespace("bayesrules", quietly = TRUE)) { data(big_word_club, package = "bayesrules") dat <- big_word_club dat$school_id <- factor(dat$school_id) dat <- subset( dat, !is.na(score_ppvt) & !is.na(invalid_ppvt) & invalid_ppvt == 0L & complete.cases(dat[, c( "score_ppvt", "distracted_a1", "distracted_ppvt", "private_school", "title1", "free_reduced_lunch", "school_id" )]) ) form_lmer <- score_ppvt ~ private_school + title1 + free_reduced_lunch + distracted_ppvt + distracted_a1 + free_reduced_lunch:distracted_a1 + (1 + distracted_ppvt + distracted_a1 || school_id) design0 <- model_setup(form_lmer, data = dat) dat <- subset(dat, school_id %in% names(design0$groupef.rank)[design0$groupef.rank]) dat$school_id <- droplevels(dat$school_id) ps <- Prior_Setup_GLMM( form_lmer, data = dat, pop.pwt = 0.01, dispformula = ~school_id, group.max_disp_perc = 0.8, group.dispersion.pwt = 0.1, group.alpha_target = NULL ) disp_pf <- dGamma_list(ps) ## Group names (school_id levels) g <- names(disp_pf) utils::head(g, 10L) ## Print a subset (recommended when J is large) print(disp_pf, groups = c("2", "18", "47")) ## Variation: rebuild truncation windows; print one group disp_pf_tight <- dGamma_list(ps, max_disp_perc_measurement = 0.9) print(disp_pf_tight, groups = "2") ## Downstream attributes (not shown by print.dGamma_list) names(attributes(disp_pf)) utils::head(attr(disp_pf, "window_diagnostics"), 3L) }
Computes the directional tail probability based on posterior draws and prior mean, using whitening transformation and projection onto the direction of disagreement. This diagnostic identifies directional disagreement between posterior and prior, and is especially useful for visualizing rejection regions in whitened space. The whitening uses Mahalanobis distance (Mahalanobis 1936) in posterior-precision-scaled coordinates.
directional_tail(fit, mu0 = NULL)directional_tail(fit, mu0 = NULL)
fit |
A fitted model object of class 'glmb' or 'lmb' |
mu0 |
An optional argument containing a reference vector relative to which the directional tail is computed. Defaults to the prior mean. |
An object of class 'directional_tail' containing:
mahalanobis_shift |
Measures the standardized Mahalanobis distance between the posterior and prior means, using posterior precision for scaling. In the Gaussian case, this directly determines the directional tail probability via Phi(-||w||). |
p_directional |
Directional tail probability (proportion of draws in the direction of disagreement) |
delta |
Mean shift in whitened space |
draws |
List containing whitened draws, raw draws, and tail flags |
Mahalanobis PC (1936).
“On the generalized distance in statistics.”
Proceedings of the National Institute of Sciences of India, 2(1), 49–55.
Nygren K (2025).
“Chapter A04: Directional Tail Diagnostics for Prior-Posterior Disagreement.”
Vignette in the glmbayes R package.
R vignette name: Chapter-A04.
lmerb / glmerb fitsBy default returns fixef from the reference
lmer/glmer fit (population-level regression coefficients).
Use type = "hyper" or "hyper.mean" for Block~2 hyperparameter
summaries stored on the Bayesian fit.
## S3 method for class 'lmerb' fixef(object, type = c("mer", "hyper", "hyper.mean"), ...) ## S3 method for class 'glmerb' fixef(object, type = c("mer", "hyper", "hyper.mean"), ...)## S3 method for class 'lmerb' fixef(object, type = c("mer", "hyper", "hyper.mean"), ...) ## S3 method for class 'glmerb' fixef(object, type = c("mer", "hyper", "hyper.mean"), ...)
object |
An |
type |
|
... |
Ignored. |
A named numeric vector.
glmb is used to fit Bayesian generalized linear models, specified by giving a symbolic descriptions of
the linear predictor, the error distribution, and the prior distribution.
glmb( formula, family = binomial, pfamily = dNormal(mu, Sigma, dispersion = 1), n = 1000, data, weights, use_parallel = TRUE, use_opencl = FALSE, verbose = FALSE, subset, offset, na.action, Gridtype = 2, n_envopt = NULL, start = NULL, etastart, mustart, control = list(...), model = TRUE, method = "glm.fit", x = FALSE, y = TRUE, contrasts = NULL, ... )glmb( formula, family = binomial, pfamily = dNormal(mu, Sigma, dispersion = 1), n = 1000, data, weights, use_parallel = TRUE, use_opencl = FALSE, verbose = FALSE, subset, offset, na.action, Gridtype = 2, n_envopt = NULL, start = NULL, etastart, mustart, control = list(...), model = TRUE, method = "glm.fit", x = FALSE, y = TRUE, contrasts = NULL, ... )
formula |
an object of class |
family |
a description of the error distribution and link
function to be used in the model. For |
pfamily |
a description of the prior distribution and associated constants to be used in the model. This
should be a pfamily function (see |
n |
number of draws to generate. If |
data |
an optional data frame, list or environment (or object
coercible by |
weights |
an optional vector of ‘prior weights’ to be used
in the fitting process. Should be |
use_parallel |
Logical. Whether to use parallel processing during simulation. |
use_opencl |
Logical. Whether to use OpenCL acceleration during Envelope construction. |
verbose |
Logical. Whether to print progress messages. |
subset |
an optional vector specifying a subset of observations to be used in the fitting process. |
offset |
this can be used to specify an a priori known component to be included in the linear
predictor during fitting. This should be |
na.action |
a function which indicates what should happen when the data contain |
Gridtype |
an optional argument specifying the method used to determine the number of tangent points used to construct the enveloping function. |
n_envopt |
Effective sample size passed to EnvelopeOpt for grid
construction. Defaults to match |
start |
starting values for the parameters in the linear predictor. |
etastart |
starting values for the linear predictor. |
mustart |
starting values for the vector of means. |
control |
a list of parameters for controlling the fitting
process. For |
model |
a logical value indicating whether model frame should be included as a component of the returned value. |
method |
the method to be used in fitting the model. The default
method User-supplied fitting functions can be supplied either as a function
or a character string naming a function, with a function which takes
the same arguments as |
x, y
|
For For |
contrasts |
an optional list. See the |
... |
For For |
The function glmb is a Bayesian version of the classical
glm function. The original R implementation of
glm was written by Simon Davies (under Ross Ihaka at the
University of Auckland) and has since been extensively rewritten by
members of the R Core Team; its design was inspired by the S
function described in (Hastie and Pregibon 1992), which in turn relies on the
formula framework described in (Wilkinson and Rogers 1973).
Setup (including the use of formulas and families) mirrors that of glm but adds a
required pfamily argument to specify the prior distribution. The design of the pfamily family of
functions was created by Kjell Nygren and is modeled on how glm
uses family to specify the likelihood.
For any implemented combination of family, link, and pfamily,
glmb generates independent draws from the posterior density-
no MCMC chains are required. Results can be printed or summarized
with methods that mirror those for glm (e.g.\ print.glmb,
summary.glmb), as well as all the usual glm/lm
generics (predict, residuals, etc.).
A helper, Prior_Setup, assists users in choosing prior
parameters. It ships with sensible defaults but also allows full
customization. In particular, the default for dNormal is a
reparameterization of Zellner's g-prior (Zellner 1986).
Currently supported response families are
gaussian (identity link), poisson and quasipoisson
(log link), gamma (log link), and binomial and
quasibinomial (logit, probit, cloglog). All families support a
dNormal prior; the Gaussian family also offers
dNormalGamma and dIndependent_Normal_Gamma.
Two conjugate pfamilies add closed-form IID sampling for intercept-only
models with an identity link: dBeta for
binomial(link = "identity") (Beta–Binomial conjugacy) and
dGamma(Inv_Dispersion = FALSE) for
poisson(link = "identity") and Gamma(link = "identity")
(Gamma–Poisson and Gamma–Gamma rate conjugacy).
For dispersion estimation with fixed coefficients, dGamma
(default Inv_Dispersion = TRUE) places a Gamma prior on the
inverse dispersion for Gaussian and Gamma(log) models.
For the Gaussian family, draws under dNormal and
dNormalGamma come from posterior distributions resulting from conjugate
prior distributions (Raiffa and Schlaifer 1961). For all other priors or response families,
we use an accept-reject sampler built on the likelihood-subgradient envelope
method (Nygren and Nygren 2006). The
Gridtype argument controls how many tangent points are used
in the envelope-trading off envelope tightness against construction
cost-and iters reports candidate counts before acceptance.
By default, glmb draws n = 1000 samples, uses parallel
CPU simulation, and-if use_opencl = TRUE-GPU-accelerated
envelope building. "glmb" comes with many of the same kinds of method functions
that come with "glm" and "lm", so you can still call extractAIC,
fitted.values, or any other standard method.
The lmb function is a Bayesian version of the lm function that can
be used to estimate models from the Gaussian family without the need for a family argument.
rglmb and rlmb are functions with more minimalistic interfaces for estimating the same
models without most of the internal overhead (these functions are called internally by glmb and lmb).
The reduced overhead may be beneficial for Gibbs sampling implementations.
glmb returns an object of class "glmb". The function summary (i.e.,
summary.glmb) can be used to obtain or print a summary of the results. The generic accessor functions
coefficients, fitted.values, residuals, and extractAIC can be used
to extract various useful features of the value returned by glmb.
An object of class "glmb" is a list containing at least the following components:
glm |
an object of class |
coefficients |
a matrix of dimension |
coef.means |
a vector of |
coef.mode |
a vector of |
dispersion |
Either a constant provided as part of the call, or a vector of length |
Prior |
A list with the priors specified for the model in question. Items in list may vary based on the type of prior |
fitted.values |
a matrix of dimension |
family |
the |
linear.predictors |
an |
deviance |
an |
pD |
An Estimate for the effective number of parameters |
Dbar |
Expected value for minus twice the log-likelihood function |
Dthetabar |
Value of minus twice the log-likelihood function evaluated at the mean value for the coefficients |
DIC |
Estimated Deviance Information criterion |
prior.weights |
a vector of weights specified or implied by the model |
y |
a vector with the dependent variable |
x |
a matrix with the implied design matrix for the model |
model |
if requested (the default),the model frame |
call |
the matched call |
formula |
the formula supplie |
terms |
the |
data |
the |
famfunc |
Family functions used during estimation process |
iters |
an |
contrasts |
(where relevant) the contrasts used. |
xlevels |
(where relevant) a record of the levels of the factors used in fitting |
digits |
the number of significant digits to use when printing. |
In addition, non-empty fits will have (yet to be implemented) components qr, R
and effects relating to the final weighted linear fit for the posterior mode.
Objects of class "glmb" are normall of class c("glmb","glm","lm"),
that is inherit from classes glm and lm and well-designed
methods from those classed will be applied when appropriate.
If a binomial glmb model was specified by giving a two-column
response, the weights returned by prior.weights are the total number of
cases (factored by the supplied case weights) and the component of y
of the result is the proportion of successes.
Hastie T~J, Pregibon D (1992).
“Generalized Linear Models.”
In Chambers J~M, Hastie T~J (eds.), Statistical Models in S, chapter 6.
Wadsworth & Brooks/Cole, Belmont, CA.
Nygren K~N, Nygren L~M (2006).
“Likelihood Subgradient Densities.”
Journal of the American Statistical Association, 101(475), 1144–1156.
doi:10.1198/016214506000000357.
Raiffa H, Schlaifer R (1961).
Applied Statistical Decision Theory.
Clinton Press, Inc., Boston.
Wilkinson GN, Rogers CE (1973).
“Symbolic Descriptions of Factorial Models for Analysis of Variance.”
Applied Statistics, 22(3), 392–399.
doi:10.2307/2346786.
Zellner A (1986).
“On Assessing Prior Distributions and Bayesian Regression Analysis with g‐Prior Distributions.”
In Goel P~K, Zellner A (eds.), Bayesian Inference and Decision Techniques: Essays in Honor of Bruno de Finetti, volume 6 of Studies in Bayesian Econometrics and Statistics, 233–243.
Elsevier.
Dobson A~J (1990).
An Introduction to Generalized Linear Models.
Chapman and Hall, London.
lm, glm, family, formula
for classical modeling functions, family objects, and formula syntax
pfamily for documentation of pfamily functions used to specify priors.
Prior_Setup, Prior_Check for functions used to initialize and to check priors,
EnvelopeBuild for envelope construction methods.
Further reading: (Nygren and Nygren 2006); (Nygren 2025, 2025); OpenCL/GPU: (Nygren 2025, 2025).
summary.glmb, predict.glmb, residuals.glmb, simulate.glmb,
extractAIC.glmb, dummy.coef.glmb and methods(class="glmb") for glmb
and the methods and generic functions for classes glm and lm from which class glmb inherits.
glmbayes Modeling Functions
lmb()
set.seed(333) ## Dobson (1990) Page 93: Randomized Controlled Trial : counts <- c(18,17,15,20,10,20,25,13,12) outcome <- gl(3,1,9) treatment <- gl(3,3) print(d.AD <- data.frame(treatment, outcome, counts)) ## Classical Model glm.D93 <- glm(counts ~ outcome + treatment, family = poisson(link=log)) summary(glm.D93) ## Poisson Prior and Model ps=Prior_Setup(counts ~ outcome + treatment,family = poisson()) mu=ps$mu V=ps$Sigma # Step 2: Call the glmb function glmb.D93<-glmb(counts ~ outcome + treatment, family=poisson(), pfamily=dNormal(mu=mu,Sigma=V)) summary(glmb.D93) # Menarche Binomial Data Example data(menarche,package="MASS") Age2=menarche$Age-13 ## Classical Model glm.out<-glm(cbind(Menarche, Total-Menarche) ~ Age2, family=binomial(logit), data=menarche) summary(glm.out) ## Logit Prior and Model ps1=Prior_Setup(cbind(Menarche, Total-Menarche) ~ Age2,family=binomial(logit), data=menarche) glmb.out1<-glmb(cbind(Menarche, Total-Menarche) ~ Age2, family=binomial(logit),pfamily=dNormal(mu=ps1$mu,Sigma=ps1$Sigma), data=menarche) summary(glmb.out1) ## Posterior mean fitted probabilities on the response scale (see also ## \code{vignette("Chapter-05", package = "glmbayes")}) require(graphics) pred1 <- predict(glmb.out1, type = "response") pred1_m <- colMeans(pred1) plot( Menarche / Total ~ Age, data = menarche, main = "Proportion with menarche (data and posterior mean fit)" ) lines(menarche$Age, pred1_m, col = "blue", lwd = 2) ## Probit Prior and Model ps2=Prior_Setup(cbind(Menarche, Total-Menarche) ~ Age2,family=binomial(probit), data=menarche) glmb.out2<-glmb(cbind(Menarche, Total-Menarche) ~ Age2, family=binomial(probit),pfamily=dNormal(mu=ps2$mu,Sigma=ps2$Sigma), data=menarche) summary(glmb.out2) ## clog-log Prior and Model ps3=Prior_Setup(cbind(Menarche, Total-Menarche) ~ Age2,family=binomial(cloglog), data=menarche) glmb.out3<-glmb(cbind(Menarche, Total-Menarche) ~ Age2, family=binomial(cloglog),pfamily=dNormal(mu=ps3$mu,Sigma=ps3$Sigma), data=menarche) summary(glmb.out3) ## Comparison of DIC Statistics DIC_Out=rbind(extractAIC(glmb.out1),extractAIC(glmb.out2),extractAIC(glmb.out3)) rownames(DIC_Out)=c("logit","probit","clog-log") colnames(DIC_Out)=c("pD","DIC") DIC_Out ### Gamma regression data(carinsca) carinsca$Merit <- ordered(carinsca$Merit) carinsca$Class <- factor(carinsca$Class) oldopt <- options(contrasts = c("contr.treatment", "contr.treatment")) Claims=carinsca$Claims Insured=carinsca$Insured Merit=carinsca$Merit Class=carinsca$Class Cost=carinsca$Cost out <- glm(Cost/Claims~Merit+Class,family=Gamma(link="log"),weights=Claims,x=TRUE) summary(out) disp=gamma.dispersion(out) ps=Prior_Setup(Cost/Claims~Merit+Class,family=Gamma(link="log"),weights=Claims) mu=ps$mu V=ps$Sigma out3 <- glmb(Cost/Claims~Merit+Class,family=Gamma(link="log"), pfamily=dNormal(mu=mu,Sigma=V,dispersion=disp),weights=Claims) summary(out) summary(out3) options(oldopt) ## glmb with dGamma prior (dispersion-only; coefficients fixed) ctl <- c(4.17, 5.58, 5.18, 6.11, 4.50, 4.61, 5.17, 4.53, 5.33, 5.14) trt <- c(4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69) group <- gl(2, 10, 20, labels = c("Ctl", "Trt")) weight <- c(ctl, trt) ps_dg <- Prior_Setup(weight ~ group, family = gaussian()) rate_dg <- if (!is.null(ps_dg$rate_gamma)) ps_dg$rate_gamma else ps_dg$rate out_glmb_dGamma <- glmb(weight ~ group, family = gaussian(), pfamily = dGamma(shape = ps_dg$shape, rate = rate_dg, beta = ps_dg$coefficients)) summary(out_glmb_dGamma)set.seed(333) ## Dobson (1990) Page 93: Randomized Controlled Trial : counts <- c(18,17,15,20,10,20,25,13,12) outcome <- gl(3,1,9) treatment <- gl(3,3) print(d.AD <- data.frame(treatment, outcome, counts)) ## Classical Model glm.D93 <- glm(counts ~ outcome + treatment, family = poisson(link=log)) summary(glm.D93) ## Poisson Prior and Model ps=Prior_Setup(counts ~ outcome + treatment,family = poisson()) mu=ps$mu V=ps$Sigma # Step 2: Call the glmb function glmb.D93<-glmb(counts ~ outcome + treatment, family=poisson(), pfamily=dNormal(mu=mu,Sigma=V)) summary(glmb.D93) # Menarche Binomial Data Example data(menarche,package="MASS") Age2=menarche$Age-13 ## Classical Model glm.out<-glm(cbind(Menarche, Total-Menarche) ~ Age2, family=binomial(logit), data=menarche) summary(glm.out) ## Logit Prior and Model ps1=Prior_Setup(cbind(Menarche, Total-Menarche) ~ Age2,family=binomial(logit), data=menarche) glmb.out1<-glmb(cbind(Menarche, Total-Menarche) ~ Age2, family=binomial(logit),pfamily=dNormal(mu=ps1$mu,Sigma=ps1$Sigma), data=menarche) summary(glmb.out1) ## Posterior mean fitted probabilities on the response scale (see also ## \code{vignette("Chapter-05", package = "glmbayes")}) require(graphics) pred1 <- predict(glmb.out1, type = "response") pred1_m <- colMeans(pred1) plot( Menarche / Total ~ Age, data = menarche, main = "Proportion with menarche (data and posterior mean fit)" ) lines(menarche$Age, pred1_m, col = "blue", lwd = 2) ## Probit Prior and Model ps2=Prior_Setup(cbind(Menarche, Total-Menarche) ~ Age2,family=binomial(probit), data=menarche) glmb.out2<-glmb(cbind(Menarche, Total-Menarche) ~ Age2, family=binomial(probit),pfamily=dNormal(mu=ps2$mu,Sigma=ps2$Sigma), data=menarche) summary(glmb.out2) ## clog-log Prior and Model ps3=Prior_Setup(cbind(Menarche, Total-Menarche) ~ Age2,family=binomial(cloglog), data=menarche) glmb.out3<-glmb(cbind(Menarche, Total-Menarche) ~ Age2, family=binomial(cloglog),pfamily=dNormal(mu=ps3$mu,Sigma=ps3$Sigma), data=menarche) summary(glmb.out3) ## Comparison of DIC Statistics DIC_Out=rbind(extractAIC(glmb.out1),extractAIC(glmb.out2),extractAIC(glmb.out3)) rownames(DIC_Out)=c("logit","probit","clog-log") colnames(DIC_Out)=c("pD","DIC") DIC_Out ### Gamma regression data(carinsca) carinsca$Merit <- ordered(carinsca$Merit) carinsca$Class <- factor(carinsca$Class) oldopt <- options(contrasts = c("contr.treatment", "contr.treatment")) Claims=carinsca$Claims Insured=carinsca$Insured Merit=carinsca$Merit Class=carinsca$Class Cost=carinsca$Cost out <- glm(Cost/Claims~Merit+Class,family=Gamma(link="log"),weights=Claims,x=TRUE) summary(out) disp=gamma.dispersion(out) ps=Prior_Setup(Cost/Claims~Merit+Class,family=Gamma(link="log"),weights=Claims) mu=ps$mu V=ps$Sigma out3 <- glmb(Cost/Claims~Merit+Class,family=Gamma(link="log"), pfamily=dNormal(mu=mu,Sigma=V,dispersion=disp),weights=Claims) summary(out) summary(out3) options(oldopt) ## glmb with dGamma prior (dispersion-only; coefficients fixed) ctl <- c(4.17, 5.58, 5.18, 6.11, 4.50, 4.61, 5.17, 4.53, 5.33, 5.14) trt <- c(4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69) group <- gl(2, 10, 20, labels = c("Ctl", "Trt")) weight <- c(ctl, trt) ps_dg <- Prior_Setup(weight ~ group, family = gaussian()) rate_dg <- if (!is.null(ps_dg$rate_gamma)) ps_dg$rate_gamma else ps_dg$rate out_glmb_dGamma <- glmb(weight ~ group, family = gaussian(), pfamily = dGamma(shape = ps_dg$shape, rate = rate_dg, beta = ps_dg$coefficients)) summary(out_glmb_dGamma)
Fits one glmb per observation block. Same row partition as
lmbBlock, but supports GLM family objects.
Counterpart to lmbBlock; see summary.bglmb for
the summary method and rNormalGLM_reg_group for Gibbs sampling.
glmbBlock( formula, block, family = gaussian(), pfamily = NULL, pfamily_list = NULL, n = 1000, data, subset, weights, na.action, offset, contrasts = NULL, model = TRUE, x = FALSE, y = TRUE, Gridtype = 2, n_envopt = NULL, use_parallel = TRUE, use_opencl = FALSE, verbose = FALSE, ... ) ## S3 method for class 'bglmb' print(x, digits = max(3, getOption("digits") - 3), ...)glmbBlock( formula, block, family = gaussian(), pfamily = NULL, pfamily_list = NULL, n = 1000, data, subset, weights, na.action, offset, contrasts = NULL, model = TRUE, x = FALSE, y = TRUE, Gridtype = 2, n_envopt = NULL, use_parallel = TRUE, use_opencl = FALSE, verbose = FALSE, ... ) ## S3 method for class 'bglmb' print(x, digits = max(3, getOption("digits") - 3), ...)
formula |
an object of class |
block |
Block partition: |
family |
a description of the error distribution and link
function to be used in the model. For |
pfamily |
Recycled to all blocks, or use |
pfamily_list |
Optional list of |
n |
number of draws to generate. If |
data |
an optional data frame, list or environment (or object
coercible by |
subset |
an optional vector specifying a subset of observations to be used in the fitting process. |
weights |
an optional vector of weights to be used in the fitting
process. Should be |
na.action |
a function which indicates what should happen when the data contain |
offset |
this can be used to specify an a priori known component to be
included in the linear predictor during fitting. This should be |
contrasts |
an optional list. See the |
model, x, y
|
For |
Gridtype |
an optional argument specifying the method used to determine the number of tangent points used to construct the enveloping function. |
n_envopt |
Effective sample size passed to EnvelopeOpt for grid
construction. Defaults to match |
use_parallel |
Logical. Whether to use parallel processing during simulation. |
use_opencl |
Logical. Whether to use OpenCL acceleration during Envelope construction. |
verbose |
Logical. Whether to print progress messages. |
... |
For |
digits |
Number of significant digits to use when printing. |
A named list of class "bglmb" (list of "glmb" fits from glmbayes).
glmbBlock(): glmb fit per row block.
glmbayes Modeling Functions
lmbBlock()
Bayesian generalized linear mixed-effects model fit
glmerb( formula, data = NULL, family = gaussian(), pfamily_list, dispersion_ranef = NULL, dispformula = ~1, n = 1000L, gap_tol = 0.0196, mode_gap_max = 1, tv_tol = 0.01, simulate = TRUE, REML = TRUE, control = NULL, start = NULL, verbose = 0L, subset, weights, na.action, offset, contrasts = NULL, devFunOnly = FALSE, fixef = NULL, progbar = FALSE, sim_method = "DEFAULT", ... ) ## S3 method for class 'glmerb' print( x, digits = max(3L, getOption("digits") - 3L), sweep_history = FALSE, sweep_history_stage = c("main", "pilot", "both"), max_sweeps = Inf, ... )glmerb( formula, data = NULL, family = gaussian(), pfamily_list, dispersion_ranef = NULL, dispformula = ~1, n = 1000L, gap_tol = 0.0196, mode_gap_max = 1, tv_tol = 0.01, simulate = TRUE, REML = TRUE, control = NULL, start = NULL, verbose = 0L, subset, weights, na.action, offset, contrasts = NULL, devFunOnly = FALSE, fixef = NULL, progbar = FALSE, sim_method = "DEFAULT", ... ) ## S3 method for class 'glmerb' print( x, digits = max(3L, getOption("digits") - 3L), sweep_history = FALSE, sweep_history_stage = c("main", "pilot", "both"), max_sweeps = Inf, ... )
formula |
Mixed-model formula (single grouping factor; same constraints
as |
data |
Data frame containing all variables in |
family |
A |
pfamily_list |
Required named list of
|
dispersion_ranef |
Observation-level measurement dispersion, treated
as known during sampling. For families with a dispersion parameter
(e.g. |
dispformula |
One-sided formula selecting the measurement-dispersion
structure: |
n |
Number of iid draws per group (default |
gap_tol |
Legacy mode–mean gap tolerance. When |
mode_gap_max |
Maximum per-coordinate mode–mean gap (in posterior
standard deviation units) used to calibrate pilot inner sweeps
(default |
tv_tol |
Total variation tolerance per stored draw, in (0, 1)
(default |
simulate |
Logical (default |
REML |
Logical; passed to |
control |
Optional |
start |
Optional starting values; passed to |
verbose |
Verbosity flag; passed to |
subset |
Optional subset; passed to |
weights |
Optional weights; passed to |
na.action |
Missing-data handler; passed to |
offset |
Optional offset; passed to |
contrasts |
Optional contrasts; passed to |
devFunOnly |
If |
fixef |
Optional named list of hyper-parameter vectors (Block 2 state).
When |
progbar |
Logical. When |
sim_method |
Sampling engine: |
... |
Reserved for future use. |
x |
Object of class |
digits |
Number of significant digits to use when printing. |
sweep_history |
If |
sweep_history_stage |
Which stage to print when |
max_sweeps |
Passed to |
Draft entry point for lmebayes GLMM models with a glmer-like
interface, analogous to lmerb for Gaussian responses and to
glmb for fixed-effects GLMs.
It takes the same arguments as lmerb plus family.
With family = gaussian() it fits the same model as
lmerb, differing only in that the embedded reference fit
comes from glmer rather than lmer.
The model being estimated.
For group with link ,
where collects the level-2 fixed effects, is
diagonal with entries (one per random-effect coefficient),
and is built by model_setup from the
level-2 predictors. As in lmerb, only a single grouping
factor and uncorrelated (||) random-effect terms are supported.
Marginalizing recovers the usual glmer form with
, , and mean-zero random
effects .
Supported families are gaussian(), poisson(),
binomial(), and Gamma() with the standard links. What they
have in common is a log-concave likelihood, which is the condition
the envelope accept-reject sampler requires (Nygren and Nygren 2006).
Because the link is applied within groups, and
live on the link scale — log rates for poisson(),
log-odds for binomial() — so any summary on the response scale
must apply the inverse link to the draws rather than to their averages.
What the priors decide.
pfamily_list carries one component per random-effect coefficient
and declares whether the between-group variance is known
(dNormal) or estimated (dIndependent_Normal_Gamma, under a
Gamma prior on truncated to
[disp_lower, disp_upper]). dispersion_ranef is the
residual variance and applies to gaussian() only; it must be
NULL for poisson() and binomial(), which have no
free dispersion parameter. Both are normally built from one
Prior_Setup_GLMM object.
What the returned draws are.
As in lmerb, glmerb runs n independent
replicate chains and stores one draw from each, taken at that chain's
final sweep. The draws are therefore mutually independent — an iid
sample, not an autocorrelated trace — so ordinary Monte Carlo standard
errors apply and there is nothing to thin or diagnose for
autocorrelation.
The accompanying guarantee is weaker here than for a Gaussian response,
and the difference matters. With family = gaussian() and fixed
variances the posterior is exactly normal and the sweep count is
certified. A Poisson or binomial posterior is not normal, so the same
calibration is applied at local curvature near the posterior mode. The
resulting m_convergence is a well-motivated sweep budget
rather than a proof, and tv_tol should be read that way.
Tightening it is cheap: cost grows only logarithmically in
tv_tol.
The same fact — that the mode of a skewed posterior is not its mean —
is why a short pilot stage runs first. Its cross-chain average
gives the main chains a better starting point than the mode does, which
shortens the main stage. Pilot draws are diagnostic only and are never
part of the returned sample. gap_tol and mode_gap_max are
legacy controls used when tv_tol is NULL.
Point estimates (simulate = FALSE).
simulate = FALSE returns the joint posterior mode by
iterated conditional modes (the exact Gaussian mean when
family = gaussian()), at fixed variance-component plug-ins. No
sweeps run and no RNG is consumed, so the result is deterministic —
useful while iterating on model specification.
Object of class "glmerb": same fixef.* structure as
"lmerb", including sigma2 and sigma2.mean for
family = gaussian() (see lmerb), with additional
family, glmer (reference
glmer fit, always the plain pooled-dispersion fit
regardless of dispformula), dispformula (as supplied),
dispersion_fit (NULL unless dispformula requests
per-group dispersion, in which case the diagnostic-only
glmmTMB fit), fixef.init (main-chain start from
pilot colMeans when a pilot runs; NULL when no pilot runs),
pilot_chisq (Hotelling chi-squared test of
pilot mean vs ICM mode), gap_tol, and mode_gap_max.
Nygren, K. (2020). On the total variation distance between multivariate normal densities with applications to two-block Gibbs samplers. Unpublished manuscript.
Nygren, K. N. and Nygren, L. M. (2006). Likelihood Subgradient Densities. Journal of the American Statistical Association 101(475), 1144–1156.
Jones, G. L. and Hobert, J. P. (2001). Honest exploration of intractable probability distributions via Markov chain Monte Carlo. Statistical Science 16, 312–334.
lmerb, glmerb_posterior_mode,
glmb; demo for the full sampling workflow
(demo("Ex_14_glmerb_airbnb_small", package = "lmebayes")).
## glmerb: neighborhood random effects on bayesrules::airbnb_small (Poisson) ## ## Fast man-page example: prior setup, glmer reference fit, and ICM posterior ## mode only (simulate = FALSE). No Gibbs draws — suitable for R CMD check. ## ## The full sampling workflow (pilot stage, main draws, centering diagnostics) ## is preserved as a demo: ## demo("Ex_14_glmerb_airbnb_small", package = "lmebayes") if (requireNamespace("bayesrules", quietly = TRUE)) { data(airbnb_small, package = "bayesrules") dat <- airbnb_small dat$rating_c <- dat$rating - mean(dat$rating, na.rm = TRUE) dat$walk_c <- dat$walk_score - mean(dat$walk_score, na.rm = TRUE) dat <- dat[complete.cases(dat[, c("reviews", "rating_c", "walk_c", "neighborhood")]), ] dat$neighborhood <- droplevels(factor(dat$neighborhood)) form <- reviews ~ walk_c + rating_c + (1 + rating_c || neighborhood) ps <- Prior_Setup_GLMM(form, data = dat, family = poisson(), pop.pwt = 0.01) fit <- glmerb( form, data = dat, family = poisson(), pfamily_list = pfamily_list(ps), simulate = FALSE ) lmebayes:::print_coef_means(fit) print(fit) lme4::fixef(fit$glmer) }## glmerb: neighborhood random effects on bayesrules::airbnb_small (Poisson) ## ## Fast man-page example: prior setup, glmer reference fit, and ICM posterior ## mode only (simulate = FALSE). No Gibbs draws — suitable for R CMD check. ## ## The full sampling workflow (pilot stage, main draws, centering diagnostics) ## is preserved as a demo: ## demo("Ex_14_glmerb_airbnb_small", package = "lmebayes") if (requireNamespace("bayesrules", quietly = TRUE)) { data(airbnb_small, package = "bayesrules") dat <- airbnb_small dat$rating_c <- dat$rating - mean(dat$rating, na.rm = TRUE) dat$walk_c <- dat$walk_score - mean(dat$walk_score, na.rm = TRUE) dat <- dat[complete.cases(dat[, c("reviews", "rating_c", "walk_c", "neighborhood")]), ] dat$neighborhood <- droplevels(factor(dat$neighborhood)) form <- reviews ~ walk_c + rating_c + (1 + rating_c || neighborhood) ps <- Prior_Setup_GLMM(form, data = dat, family = poisson(), pop.pwt = 0.01) fit <- glmerb( form, data = dat, family = poisson(), pfamily_list = pfamily_list(ps), simulate = FALSE ) lmebayes:::print_coef_means(fit) print(fit) lme4::fixef(fit$glmer) }
Returns whether this lmebayes installation was built with OpenCL
support (link-time / USE_OPENCL). This is independent of
has_opencl() in opencltools,
which probes the host runtime (ICD, drivers, headers).
For workstation diagnostics (drivers, PATH, ICD), use
diagnose_glmbayes() from
opencltools. For kernel loading from inst/cl/, use
load_kernel_source() with the
appropriate package argument once opencltools is installed.
has_opencl()has_opencl()
Logical scalar.
opencltools, glmbayes
lmb is used to fit Bayesian linear models, specified by giving a symbolic descriptions of the linear
predictor and the prior distribution.
lmb( formula, pfamily, n = 1000, data, subset, weights, na.action, method = "qr", model = TRUE, x = TRUE, y = TRUE, qr = TRUE, singular.ok = TRUE, contrasts = NULL, offset, Gridtype = 2, n_envopt = NULL, use_parallel = TRUE, use_opencl = FALSE, verbose = FALSE, ... )lmb( formula, pfamily, n = 1000, data, subset, weights, na.action, method = "qr", model = TRUE, x = TRUE, y = TRUE, qr = TRUE, singular.ok = TRUE, contrasts = NULL, offset, Gridtype = 2, n_envopt = NULL, use_parallel = TRUE, use_opencl = FALSE, verbose = FALSE, ... )
formula |
an object of class |
pfamily |
a description of the prior distribution and associated constants to be used in the model.
For a single-response formula this should be a single |
n |
number of draws to generate. If |
data |
an optional data frame, list or environment (or object
coercible by |
subset |
an optional vector specifying a subset of observations to be used in the fitting process. |
weights |
an optional vector of weights to be used in the fitting
process. Should be |
na.action |
a function which indicates what should happen when the data contain |
method |
the method to be used in fitting the classical model during a call to |
model, x, y, qr
|
logicals. If |
singular.ok |
logical. If |
contrasts |
an optional list. See the |
offset |
this can be used to specify an a priori known component to be
included in the linear predictor during fitting. This should be |
Gridtype |
an optional argument specifying the method used to determine the number of tangent points used to construct the enveloping function. |
n_envopt |
Effective sample size passed to EnvelopeOpt for grid
construction. Defaults to match |
use_parallel |
Logical. Whether to use parallel processing during simulation. |
use_opencl |
Logical. Whether to use OpenCL acceleration during Envelope construction. |
verbose |
Logical. Whether to print progress messages. |
... |
For |
The function lmb is a Bayesian extension of the classical lm function.
It retains the familiar formula interface and model setup used in lm, while introducing
posterior simulation and prior specification via the pfamily argument. Internally,
lmb calls lm to obtain the classical least squares fit, then generates independent
draws from the posterior distribution using either multivariate normal simulation (for Gaussian priors)
or accept-reject sampling via likelihood-subgradient envelopes (Nygren and Nygren 2006).
The symbolic formula interface follows (Wilkinson and Rogers 1973), and the overall design of lm was inspired by the S system
(Chambers 1992). lmb comes with many of the same types of generic methods
that are available to lm and glm, including predict, residuals, extractAIC,
and summary. Many of these are inherited from glmb.
Prior specification is handled via the pfamily argument, which defines the prior mean,
covariance, and dispersion. The design of the pfamily family of functions was created by Kjell Nygren
and is modeled on how glm uses family to specify the likelihood. A helper function,
Prior_Setup, assists users in choosing prior parameters. It ships with sensible defaults but
also allows full customization. All models support the dNormal prior; the Gaussian family also supports
dNormalGamma and dIndependent_Normal_Gamma, which allow for more flexible prior structures
including independent priors on variance components.
Posterior draws are generated using the prior specification provided via pfamily. For Gaussian models
with conjugate priors, draws are obtained directly from the posterior distribution using standard simulation
procedures for multivariate normal densities (Raiffa and Schlaifer 1961). For non-conjugate setups, the
function uses envelope-based accept-reject sampling, where the Gridtype parameter controls the granularity
of the envelope construction. The number of candidates generated before acceptance is returned in the
iters component.
The output includes both the classical lm fit and Bayesian diagnostics such as the Deviance
Information Criterion (DIC), effective number of parameters (pD), and posterior summaries.
The DIC, introduced by (Spiegelhalter et al. 2002), provides a Bayesian analog to AIC
by balancing model fit and complexity using posterior expectations. This dual structure allows users
to compare classical and Bayesian fits side-by-side, and to leverage familiar modeling workflows
while gaining access to richer inferential tools.
The lmb function is a specialized version of glmb for Gaussian models,
and does not require a family argument. For conjugate models, it uses standard simulation methods for
posterior draws, avoiding the need for envelope construction or subgradient sampling.
Like glmb, it returns objects compatible with many standard methods from lm and glm,
including extractAIC, fitted.values, and residuals.
For more minimalistic workflows, rlmb and rglmb offer stripped-down interfaces
for posterior sampling without the overhead of full model objects. rlmb is called from within
lmb. The functions rlmb might be useful in Gibbs sampling or simulation-heavy contexts.
lmb returns an object of class "lmb". The function summary (i.e.,
summary.glmb) can be used to obtain or print a summary of the results. The generic accessor functions
coefficients, fitted.values, residuals, and extractAIC can be used
to extract various useful features of the value returned by lmb.
An object of class "lmb" is a list containing at least the following components:
lm |
an object of class |
coefficients |
a matrix of dimension |
coef.means |
a vector of |
coef.mode |
a vector of |
dispersion |
Either a constant provided as part of the call, or a vector of length |
Prior |
A list with the priors specified for the model in question. Items in list may vary based on the type of prior |
residuals |
a matrix of dimension |
fitted.values |
a matrix of dimension |
linear.predictors |
an |
deviance |
an |
pD |
An Estimate for the effective number of parameters |
Dbar |
Expected value for minus twice the log-likelihood function |
Dthetabar |
Value of minus twice the log-likelihood function evaluated at the mean value for the coefficients |
DIC |
Estimated Deviance Information criterion |
weights |
a vector of weights specified or implied by the model |
prior.weights |
a vector of weights specified or implied by the model |
y |
a vector of observations of length |
x |
a design matrix of dimension |
model |
if requested (the default),the model frame |
call |
the matched call |
formula |
the formula supplied |
terms |
the |
data |
the |
famfunc |
family functions used during estimation and post processing |
iters |
an |
contrasts |
(where relevant) the contrasts used. |
xlevels |
(where relevant) a record of the levels of the factors used in fitting |
pfamily |
the prior family specified |
digits |
the number of significant digits to use when printing. |
In addition, non-empty fits will have (yet to be implemented) components qr, R
and effects relating to the final weighted linear fit for the posterior mode.
Objects of class "lmb" are normally of class c("lmb","glmb","glm","lm"),
that is inherit from classes glmb. glm and lm and well-designed
methods from those classed will be applied when appropriate.
Chambers JM (1992).
“Linear Models.”
In Chambers JM, Hastie TJ (eds.), Statistical Models in S, chapter 4, 85–124.
Wadsworth & Brooks/Cole, Pacific Grove, CA.
Nygren K~N, Nygren L~M (2006).
“Likelihood Subgradient Densities.”
Journal of the American Statistical Association, 101(475), 1144–1156.
doi:10.1198/016214506000000357.
Raiffa H, Schlaifer R (1961).
Applied Statistical Decision Theory.
Clinton Press, Inc., Boston.
Spiegelhalter DJ, Best NG, Carlin BP, van der Linde A (2002).
“Bayesian Measures of Model Complexity and Fit.”
Journal of the Royal Statistical Society: Series B (Statistical Methodology), 64(4), 583–639.
doi:10.1111/1467-9868.00353.
Wilkinson GN, Rogers CE (1973).
“Symbolic Descriptions of Factorial Models for Analysis of Variance.”
Applied Statistics, 22(3), 392–399.
doi:10.2307/2346786.
The classical modeling functions lm and glm.
glmb, rglmb, rlmb
for related Bayesian GLM/linear interfaces;
EnvelopeBuild for envelope construction when accept–reject sampling is used.
pfamily for documentation of pfamily functions used to specify priors.
Prior_Setup, Prior_Check for functions used to initialize and to check priors,
Further reading: (Nygren and Nygren 2006); (Nygren 2025, 2025); independent Normal–Gamma sampler: (Nygren 2025).
summary.glmb, predict.glmb, simulate.glmb,
extractAIC.glmb, dummy.coef.glmb and methods(class="glmb") for methods
inherited from class glmb and the methods and generic functions for classes glm and
lm from which class lmb also inherits.
glmbayes Modeling Functions
glmb()
## Annette Dobson (1990) "An Introduction to Generalized Linear Models". ## Page 9: Plant Weight Data. ctl <- c(4.17, 5.58, 5.18, 6.11, 4.50, 4.61, 5.17, 4.53, 5.33, 5.14) trt <- c(4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69) group <- gl(2, 10, 20, labels = c("Ctl", "Trt")) weight <- c(ctl, trt) ps <- Prior_Setup(weight ~ group, gaussian()) # Prior_Setup supplies the prior mean and covariance components used below. ## Classical model lm.D9 <- lm(weight ~ group, x = TRUE, y = TRUE) summary(lm.D9) vcov(lm.D9) s <- summary(lm.D9) disp_classical <- s$sigma^2 cat("Classical lm dispersion (sigma^2 = RSS/(n-p)):", disp_classical, "\n") ## Conjugate Normal Prior (fixed dispersion) lmb.D9 <- lmb( weight ~ group, pfamily = dNormal(mu = ps$mu, ps$Sigma, dispersion = ps$dispersion) ) summary(lmb.D9) vcov(lmb.D9) ## Conjugate Normal_Gamma Prior (second argument is Sigma_0 from Prior_Setup) lmb.D9_v2 <- lmb( weight ~ group, pfamily = dNormal_Gamma( ps$mu, Sigma_0 = ps$Sigma_0, shape = ps$shape, rate = ps$rate ) ) summary(lmb.D9_v2) vcov(lmb.D9_v2) ## Independent_Normal_Gamma_Prior (same mu, Sigma, rate as dNormal_Gamma; shape = ps$shape_ING) lmb.D9_v3 <- lmb( weight ~ group, dIndependent_Normal_Gamma( ps$mu, ps$Sigma, shape = ps$shape_ING, rate = ps$rate ) ) summary(lmb.D9_v3) ## anova anova(lmb.D9) ## lmb with dGamma prior (dispersion-only; coefficients fixed) rate_dg <- if (!is.null(ps$rate_gamma)) ps$rate_gamma else ps$rate out_lmb_dGamma <- lmb( weight ~ group, pfamily = dGamma(shape = ps$shape, rate = rate_dg, beta = ps$coefficients)) summary(out_lmb_dGamma)## Annette Dobson (1990) "An Introduction to Generalized Linear Models". ## Page 9: Plant Weight Data. ctl <- c(4.17, 5.58, 5.18, 6.11, 4.50, 4.61, 5.17, 4.53, 5.33, 5.14) trt <- c(4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69) group <- gl(2, 10, 20, labels = c("Ctl", "Trt")) weight <- c(ctl, trt) ps <- Prior_Setup(weight ~ group, gaussian()) # Prior_Setup supplies the prior mean and covariance components used below. ## Classical model lm.D9 <- lm(weight ~ group, x = TRUE, y = TRUE) summary(lm.D9) vcov(lm.D9) s <- summary(lm.D9) disp_classical <- s$sigma^2 cat("Classical lm dispersion (sigma^2 = RSS/(n-p)):", disp_classical, "\n") ## Conjugate Normal Prior (fixed dispersion) lmb.D9 <- lmb( weight ~ group, pfamily = dNormal(mu = ps$mu, ps$Sigma, dispersion = ps$dispersion) ) summary(lmb.D9) vcov(lmb.D9) ## Conjugate Normal_Gamma Prior (second argument is Sigma_0 from Prior_Setup) lmb.D9_v2 <- lmb( weight ~ group, pfamily = dNormal_Gamma( ps$mu, Sigma_0 = ps$Sigma_0, shape = ps$shape, rate = ps$rate ) ) summary(lmb.D9_v2) vcov(lmb.D9_v2) ## Independent_Normal_Gamma_Prior (same mu, Sigma, rate as dNormal_Gamma; shape = ps$shape_ING) lmb.D9_v3 <- lmb( weight ~ group, dIndependent_Normal_Gamma( ps$mu, ps$Sigma, shape = ps$shape_ING, rate = ps$rate ) ) summary(lmb.D9_v3) ## anova anova(lmb.D9) ## lmb with dGamma prior (dispersion-only; coefficients fixed) rate_dg <- if (!is.null(ps$rate_gamma)) ps$rate_gamma else ps$rate out_lmb_dGamma <- lmb( weight ~ group, pfamily = dGamma(shape = ps$shape, rate = rate_dg, beta = ps$coefficients)) summary(out_lmb_dGamma)
Fits one lmb per observation block (SAS BY-style split on
rows), sharing the same formula on each subset. Contrast with
lmb on a cbind(...) response (several response columns) and
rNormalGLM_reg_group (Gibbs conditional draws, matrix API).
lmbBlock( formula, block, pfamily = NULL, pfamily_list = NULL, n = 1000, data, subset, weights, na.action, method = "qr", model = TRUE, x = TRUE, y = TRUE, qr = TRUE, singular.ok = TRUE, contrasts = NULL, offset, Gridtype = 2, n_envopt = NULL, use_parallel = TRUE, use_opencl = FALSE, verbose = FALSE, ... ) ## S3 method for class 'blmb' print(x, digits = max(3, getOption("digits") - 3), ...)lmbBlock( formula, block, pfamily = NULL, pfamily_list = NULL, n = 1000, data, subset, weights, na.action, method = "qr", model = TRUE, x = TRUE, y = TRUE, qr = TRUE, singular.ok = TRUE, contrasts = NULL, offset, Gridtype = 2, n_envopt = NULL, use_parallel = TRUE, use_opencl = FALSE, verbose = FALSE, ... ) ## S3 method for class 'blmb' print(x, digits = max(3, getOption("digits") - 3), ...)
formula |
an object of class |
block |
Block partition: |
pfamily |
A single |
pfamily_list |
Optional list of |
n |
number of draws to generate. If |
data |
an optional data frame, list or environment (or object
coercible by |
subset |
an optional vector specifying a subset of observations to be used in the fitting process. |
weights |
an optional vector of weights to be used in the fitting
process. Should be |
na.action |
a function which indicates what should happen when the data contain |
method |
the method to be used in fitting the classical model during a call to |
model, x, y, qr
|
For |
singular.ok |
logical. If |
contrasts |
an optional list. See the |
offset |
this can be used to specify an a priori known component to be
included in the linear predictor during fitting. This should be |
Gridtype |
an optional argument specifying the method used to determine the number of tangent points used to construct the enveloping function. |
n_envopt |
Effective sample size passed to EnvelopeOpt for grid
construction. Defaults to match |
use_parallel |
Logical. Whether to use parallel processing during simulation. |
use_opencl |
Logical. Whether to use OpenCL acceleration during Envelope construction. |
verbose |
Logical. Whether to print progress messages. |
... |
For |
digits |
Number of significant digits to use when printing. |
A named list of class "blmb" (list of "lmb" fits).
lmbBlock(): Gaussian lmb fit per row block.
glmbayes Modeling Functions
glmbBlock()
## Row-block lmb (iris): BY Species — SAS-style separate regressions set.seed(42) data("iris", package = "datasets") ps_block <- Prior_SetupGroup( Sepal.Length ~ Sepal.Width + Petal.Length, group = "Species", data = iris, family = gaussian() ) out_blmb <- lmbBlock( Sepal.Length ~ Sepal.Width + Petal.Length, block = "Species", pfamily_list = pfamily_list(ps_block), data = iris, n = 150L, use_parallel = FALSE ) print(out_blmb) summary(out_blmb)## Row-block lmb (iris): BY Species — SAS-style separate regressions set.seed(42) data("iris", package = "datasets") ps_block <- Prior_SetupGroup( Sepal.Length ~ Sepal.Width + Petal.Length, group = "Species", data = iris, family = gaussian() ) out_blmb <- lmbBlock( Sepal.Length ~ Sepal.Width + Petal.Length, block = "Species", pfamily_list = pfamily_list(ps_block), data = iris, n = 150L, use_parallel = FALSE ) print(out_blmb) summary(out_blmb)
Bayesian linear mixed-effects model fit
lmerb( formula, data = NULL, pfamily_list, dispersion_ranef, dispformula = ~1, n = 1000L, tv_tol = 0.01, gap_tol = 0.0196, mode_gap_max = 1, diag_sweeps = FALSE, progbar = NULL, sim_method = "DEFAULT", simulate = TRUE, REML = TRUE, control = lme4::lmerControl(), start = NULL, verbose = 0L, subset, weights, na.action, offset, contrasts = NULL, devFunOnly = FALSE, fixef = NULL, ... ) ## S3 method for class 'lmerb' print( x, digits = max(3L, getOption("digits") - 3L), sweep_history = FALSE, sweep_history_stage = c("main", "pilot", "both"), max_sweeps = Inf, components = NULL, covariate = NULL, ... )lmerb( formula, data = NULL, pfamily_list, dispersion_ranef, dispformula = ~1, n = 1000L, tv_tol = 0.01, gap_tol = 0.0196, mode_gap_max = 1, diag_sweeps = FALSE, progbar = NULL, sim_method = "DEFAULT", simulate = TRUE, REML = TRUE, control = lme4::lmerControl(), start = NULL, verbose = 0L, subset, weights, na.action, offset, contrasts = NULL, devFunOnly = FALSE, fixef = NULL, ... ) ## S3 method for class 'lmerb' print( x, digits = max(3L, getOption("digits") - 3L), sweep_history = FALSE, sweep_history_stage = c("main", "pilot", "both"), max_sweeps = Inf, components = NULL, covariate = NULL, ... )
formula |
Mixed-model formula (single grouping factor; same constraints
as |
data |
Data frame containing all variables in |
pfamily_list |
Required named list of
|
dispersion_ranef |
Observation-level measurement dispersion
|
dispformula |
One-sided formula selecting the measurement-dispersion
structure: |
n |
Number of iid draws per group (default |
tv_tol |
Total variation tolerance per stored draw, in (0, 1)
(default |
gap_tol |
Legacy mode–mean gap tolerance for the pilot stage when
any Block~2 component uses |
mode_gap_max |
Maximum per-coordinate mode–mean gap (in posterior
standard deviation units) used to calibrate pilot inner sweeps when ING
Block~2 components run a pilot stage (default |
diag_sweeps |
Temporary diagnostic flag (ING models with pilot).
Non- |
progbar |
Logical. Show text progress bars during sampling (passed to
|
sim_method |
Sampling engine: |
simulate |
Logical (default |
REML |
Logical; passed to |
control |
|
start |
Optional starting values; passed to |
verbose |
Verbosity flag; passed to |
subset |
Optional subset; passed to |
weights |
Optional weights; passed to |
na.action |
Missing-data handler; passed to |
offset |
Optional offset; passed to |
contrasts |
Optional contrasts; passed to |
devFunOnly |
If |
fixef |
Optional named list of hyper-parameter vectors (Block 2 state).
When |
... |
Reserved for future use. |
x |
Object of class |
digits |
Number of significant digits to use when printing. |
sweep_history |
If |
sweep_history_stage |
Which stage to print when |
max_sweeps |
Passed to |
components |
Optional RE components passed to sweep-history |
covariate |
Optional covariate names passed to sweep-history |
Fits a Bayesian linear mixed model from an lmer-style formula,
returning posterior draws. It is the mixed-model counterpart of
lmb for fixed-effects models.
lmerb builds the design with model_setup, embeds the
companion lmer fit for reference, and passes both to
the lmebayesCore sampler. Priors are supplied as
pfamily_list (population level) and dispersion_ranef
(residual variance), normally built from one
Prior_Setup_GLMM object.
For group , lmerb estimates
The second line is the model for how groups differ: each group's own
coefficient vector varies around a population expectation
, and measures how widely. This is a
centered parameterization — the full coefficient
enters the likelihood, not a mean-zero deviation — so coef() on
the fit reports directly and ranef() reports the
deviation that lme4 calls
. Marginalizing recovers the usual lme4 model, with
playing the role of the fixed effects.
is diagonal because formula may only contain
uncorrelated (||) random-effect terms and a single grouping
factor; there is no random-effect correlation to estimate. See
vignette("Chapter-02", package = "lmebayes") for the formula rules.
The two prior arguments are where you declare which variances are known and which are estimated, and that choice determines what comes back.
pfamily_listOne component per random-effect coefficient.
A dNormal component fixes the between-group variance
at the supplied value; a
dIndependent_Normal_Gamma component estimates it, under a
Gamma prior on truncated to
[disp_lower, disp_upper].
dispersion_ranefThe residual variance :
a fixed scalar, a pooled dGamma() prior to estimate one
common value, a dGamma_list() to estimate a separate
per group, or a named vector of fixed per-group
values. dispformula must agree with the shape supplied.
When every variance is fixed, the posterior is exactly multivariate normal and the draws are exact. Estimating any variance makes the posterior non-Gaussian and the draws approximate in the sense described next.
lmerb does not return one long chain. It runs n
independent replicate chains and keeps one draw from
each, taken at that chain's final sweep. Two things follow:
The draws are mutually independent. They form an iid
sample, so ordinary Monte Carlo standard errors apply. Nothing needs
thinning, there is no effective sample size to compute, and trace
plots or on the returned draws are not
meaningful.
Each draw comes from a distribution deliberately made
close to the posterior. The sweep count is not a burn-in guess: it
is calibrated before sampling so that the distribution reached at the
final sweep is within tv_tol of the exact posterior in total
variation. The resulting m_convergence is reported on the fit.
n and tv_tol therefore control different errors — Monte
Carlo noise and per-draw bias — and can be set independently. To
certify all n draws jointly at level , pass
tv_tol = alpha / n; cost grows only logarithmically in
tv_tol, so this usually adds a few sweeps.
The calibration itself, and what to do if a warning reports slow
convergence, are covered in
vignette("Chapter-06", package = "lmebayes").
Object of class "lmerb": a list with the following
components (parallel to glmb and lmb).
Every draw slot holds one row per replicate chain, taken at that
chain's final sweep, so the rows are mutually independent and can be
averaged directly. Accessors fixef(), ranef(),
coef() and summary() are usually more convenient than
reaching into these fields; see
vignette("Chapter-05", package = "lmebayes").
callThe matched call.
formulaThe formula supplied.
lmerlmer fit from
model_setup (full formula), embedded as a sub-object —
analogous to glmb$glm and lmb$lm. Use
coef(fit$lmer) for per-group classical coefficients. Always
the plain pooled-dispersion lmer fit, regardless of
dispformula.
dispformulaThe dispformula supplied.
dispersion_fitNULL when dispformula = ~1;
otherwise the diagnostic-only glmmTMB fit with
per-group residual dispersion (dispformula). Identical to
glmmTMB below; kept for backward compatibility.
glmmTMBSame value as dispersion_fit (an alias):
NULL when dispformula = ~1, otherwise the
glmmTMB per-group-dispersion reference fit.
Prefer this field name in new code; downstream consumers that need
a per-group-dispersion-aware reference fit (e.g.\
summary.lmerb()) use glmmTMB when non-NULL and
fall back to lmer otherwise.
priorNormalized prior container: pfamily_list
(as supplied, reordered to the RE coefficient names),
dispersion_ranef, the reconstructed Sigma_ranef, and
the per-component prior_list (mu_fixef,
Sigma_fixef, dispersion_fixef) — analogous to
glmb$Prior.
model_setupThe model_setup object built
inside lmerb from formula and data.
popef.modeNamed list of exact posterior mode (= mean,
since the joint posterior is Gaussian) vectors for the level-2 fixed
effects , computed by
lmerb_posterior_mean (ICM).
groupef.mode numeric matrix
of exact posterior mode random effects from ICM. Rows are group
levels (levels(design$group)); columns are
design$re_coef_names.
popefNamed list of matrices of
population draws , one per random-effect coefficient.
NULL when simulate = FALSE.
popef.meansNamed list of posterior mean vectors, the
column means of popef. NULL when
simulate = FALSE.
groupefdata.frame with n * J rows:
draw, the grouping-factor column, and one column per RE
variable. These are full group coefficients ,
comparable to lme4::coef(), not mean-zero deviations.
Average over draw within each group for posterior
means (see Examples). NULL when simulate = FALSE.
popef.dispersion matrix of
the between-group variances at each stored draw:
sampled values for dIndependent_Normal_Gamma components,
constant columns for dNormal components, where the value was
treated as known. NULL when simulate = FALSE.
popef.dispersion.meanNamed vector of posterior means of
(colMeans(fixef.dispersion)). NULL when
simulate = FALSE.
group.dispersionObservation-level residual variance :
a scalar when dispersion_ranef is fixed, a length-n vector
of final-sweep draws when dispersion_ranef is dGamma();
NULL when not applicable. When simulate = FALSE, the
prior plug-in scalar (prior$dispersion_ranef) is returned for
Gaussian models.
group.dispersion.meanPosterior mean of ; equals
group.dispersion when fixed. NULL when group.dispersion is
NULL.
draw_engineName of the Block~1 sampling engine used
(e.g. "rGLMM_sweep_ing_block1_ind"). NULL when
simulate = FALSE.
fixef.iters matrix of the
total number of Block~2 candidates generated per stored draw,
summed over the m_convergence inner sweeps.
dIndependent_Normal_Gamma components count envelope
accept-reject candidates until acceptance (as iters in
glmbayes samplers); dNormal components count exactly
one conjugate draw per sweep. NULL when
simulate = FALSE.
fixef.iters.meanNamed vector of average candidates per
accepted Block~2 draw
(colMeans(fixef.iters)/m_convergence; equals 1 for
dNormal components, and approximately the reciprocal
acceptance rate of the ING envelope sampler otherwise).
NULL when simulate = FALSE.
fixef.muNumeric matrix p_re x J of Block 1 prior
means at the final Gibbs state (from
build_mu_all).
convergenceList describing the sweep-count calibration:
method ("exact", "exact_iid" when
sim_method_used = "DEFAULT", or "local_gaussian_mode"
for non-Gaussian glmerb), tv_tol,
lambda_star, eigenvalues, m_min (derived
minimum sweeps), and m_convergence (sweeps actually used).
NULL when simulate = FALSE.
sim_method_used"DEFAULT" (exact iid draws) or
"TWO_BLOCK_GIBBS" (two-block Gibbs sweeps), whichever engine
actually ran. NULL when simulate = FALSE.
Nygren, K. (2020). On the total variation distance between multivariate normal densities with applications to two-block Gibbs samplers. Unpublished manuscript.
Jones, G. L. and Hobert, J. P. (2001). Honest exploration of intractable probability distributions via Markov chain Monte Carlo. Statistical Science 16, 312–334.
Prior_Setup_GLMM, model_setup,
build_mu_all,
two_block_rNormal_reg,
lmerb_posterior_mean,
rNormal_reg_group,
lmb, glmb
## lmerb: school random effects on bayesrules::big_word_club ## ## Random intercept and distracted_ppvt random slope by school; the school's ## share of free/reduced-lunch students (free_reduced_lunch, constant within ## school) is a level-2 predictor of the school intercepts. ## ## Fast man-page example: prior setup, lmer reference fit, and ICM posterior ## mean only (simulate = FALSE). No Gibbs draws — suitable for R CMD check. ## ## The full sampling workflows are preserved as demos: ## demo("Ex_14_lmerb_Sleepstudy", package = "lmebayes") # lme4-style small model ## demo("Ex_12_lmerb_BigWordClub", package = "lmebayes") # full formula if (requireNamespace("bayesrules", quietly = TRUE)) { data(big_word_club, package = "bayesrules") dat <- big_word_club dat$school_id <- factor(dat$school_id) dat <- subset(dat, !is.na(score_ppvt) & !is.na(invalid_ppvt) & invalid_ppvt == 0L) dat <- dat[complete.cases(dat[, c("score_ppvt", "distracted_ppvt", "free_reduced_lunch", "school_id")]), ] form <- score_ppvt ~ free_reduced_lunch + distracted_ppvt + (1 + distracted_ppvt || school_id) ## Default hyperpriors calibrated from a reference lmer fit (weak prior). ps <- Prior_Setup_GLMM(form, data = dat, pop.pwt = 0.01) fit <- lmerb( form, data = dat, pfamily_list = pfamily_list(ps), dispersion_ranef = ps$group.dispersion, simulate = FALSE ) ## lmer MLE vs ICM posterior mean (no MCMC means when simulate = FALSE). lmebayes:::print_coef_means(fit) print(fit) summary(fit) }## lmerb: school random effects on bayesrules::big_word_club ## ## Random intercept and distracted_ppvt random slope by school; the school's ## share of free/reduced-lunch students (free_reduced_lunch, constant within ## school) is a level-2 predictor of the school intercepts. ## ## Fast man-page example: prior setup, lmer reference fit, and ICM posterior ## mean only (simulate = FALSE). No Gibbs draws — suitable for R CMD check. ## ## The full sampling workflows are preserved as demos: ## demo("Ex_14_lmerb_Sleepstudy", package = "lmebayes") # lme4-style small model ## demo("Ex_12_lmerb_BigWordClub", package = "lmebayes") # full formula if (requireNamespace("bayesrules", quietly = TRUE)) { data(big_word_club, package = "bayesrules") dat <- big_word_club dat$school_id <- factor(dat$school_id) dat <- subset(dat, !is.na(score_ppvt) & !is.na(invalid_ppvt) & invalid_ppvt == 0L) dat <- dat[complete.cases(dat[, c("score_ppvt", "distracted_ppvt", "free_reduced_lunch", "school_id")]), ] form <- score_ppvt ~ free_reduced_lunch + distracted_ppvt + (1 + distracted_ppvt || school_id) ## Default hyperpriors calibrated from a reference lmer fit (weak prior). ps <- Prior_Setup_GLMM(form, data = dat, pop.pwt = 0.01) fit <- lmerb( form, data = dat, pfamily_list = pfamily_list(ps), dispersion_ranef = ps$group.dispersion, simulate = FALSE ) ## lmer MLE vs ICM posterior mean (no MCMC means when simulate = FALSE). lmebayes:::print_coef_means(fit) print(fit) summary(fit) }
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.
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, ... )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, ... )
formula |
A two-sided linear formula, as in |
data |
An optional data frame containing the variables named in
|
family |
As in |
REML |
Logical scalar, as in |
control |
|
start |
Optional starting values, passed to the reference
|
verbose |
Integer scalar passed to the reference |
subset |
An optional expression indicating the subset of rows of
|
weights |
An optional vector of prior weights, as in |
na.action |
A function indicating what should happen when |
offset |
An optional a priori known component to include in the
linear predictor, as in |
contrasts |
An optional list of contrasts (see the
|
devFunOnly |
Logical, as in |
fit_mer |
If |
dispformula |
As in |
... |
Ignored. |
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 , 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")).
Object of class "model_setup" with:
yResponse vector, length nrow(D).
weightsNumeric prior-weight vector of length
length(y) (ones when unspecified), aligned with y/
D/group after subset/na.action.
offsetNumeric offset vector of length length(y)
(zeros when unspecified and no formula offset()), aligned
with y/D/group.
DGroup-effects model matrix (n_obs x
p_re): per-observation loadings on the within-group
group-effect coefficient vector (columns groupef.names).
groupFactor 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_nameGrouping factor name.
groupef.namesGroup-effect coefficient names from
lme4; always identical to names(W) (and to
colnames(D)).
WNamed list of group-level hyper-design matrices (one
per groupef.names entry, one row per group level).
popef.moderationCross-level moderation metadata; see
extract_re_hyper_matrices.
formulaThe formula argument, as supplied.
familyThe normalized family
object.
callThe matched call.
termsThe terms object for
formula, as in lm/glm (and
glmbayes::glmb/lmb).
dataThe data argument, as supplied (NULL
when omitted), as in lm/glm (and
glmbayes::glmb/lmb).
dispformulaThe dispformula argument, as
supplied.
lmer / glmer
The reference
lmer (Gaussian) or glmer
(otherwise) fit on formula; NULL when
fit_mer = FALSE.
varcorrThe raw VarCorr object from
the reference fit.
PsiNamed numeric vector (names groupef.names) of
group-effect variances.
dispersionResidual (measurement) variance from the reference fit.
glmmTMB_fitNULL 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.rankNamed logical (names levels(group)):
TRUE if that group's is full column rank; see
check_identifiability.
groupef.estimableNamed logical, same names as
groupef.rank: TRUE if the group is additionally
estimable (family-specific MLE/dispersion check); see
check_identifiability.
groupef.glm_checkPer-group diagnostic data frame from
check_identifiability (family-dependent; may be
NULL).
popef.rankNamed logical (names groupef.names):
TRUE if W[[k]] is full column rank when restricted to
groupef.estimable groups.
popef.deficientNegation of popef.rank.
popef.rank_okScalar TRUE only when every
popef.rank entry is TRUE; required by
Prior_Setup_GLMM to derive default hyperpriors
automatically.
check_identifiability,
extract_re_hyper_matrices
## model_setup() — design + reference lmer (Ex_13b / big_word_club) ## ## Same formula and data cleaning as ## demo("Ex_13b_rLMM_estimated_dispersion_known_vcov_BigWordClub_PartVI"). ## Requires suggested package bayesrules. if (requireNamespace("bayesrules", quietly = TRUE)) { data(big_word_club, package = "bayesrules") dat <- big_word_club dat$school_id <- factor(dat$school_id) dat <- subset( dat, !is.na(score_ppvt) & !is.na(invalid_ppvt) & invalid_ppvt == 0L & complete.cases(dat[, c( "score_ppvt", "distracted_a1", "distracted_ppvt", "private_school", "title1", "free_reduced_lunch", "school_id" )]) ) form_lmer <- score_ppvt ~ private_school + title1 + free_reduced_lunch + distracted_ppvt + distracted_a1 + free_reduced_lunch:distracted_a1 + (1 + distracted_ppvt + distracted_a1 || school_id) ## Default: design matrices + reference lmer + identifiability design <- model_setup(form_lmer, data = dat) print(design) ## Full-rank filter (same as Ex_13b before prior setup / sampling) full_rank <- names(design$groupef.rank)[design$groupef.rank] dat_fr <- subset(dat, school_id %in% full_rank) dat_fr$school_id <- droplevels(dat_fr$school_id) design_fr <- model_setup(form_lmer, data = dat_fr) stopifnot(all(design_fr$groupef.rank)) print(design_fr) ## Variation: design only (no reference merMod fit) design_only <- model_setup(form_lmer, data = dat_fr, fit_mer = FALSE) stopifnot(is.null(design_only$lmer)) names(design_only$W) }## model_setup() — design + reference lmer (Ex_13b / big_word_club) ## ## Same formula and data cleaning as ## demo("Ex_13b_rLMM_estimated_dispersion_known_vcov_BigWordClub_PartVI"). ## Requires suggested package bayesrules. if (requireNamespace("bayesrules", quietly = TRUE)) { data(big_word_club, package = "bayesrules") dat <- big_word_club dat$school_id <- factor(dat$school_id) dat <- subset( dat, !is.na(score_ppvt) & !is.na(invalid_ppvt) & invalid_ppvt == 0L & complete.cases(dat[, c( "score_ppvt", "distracted_a1", "distracted_ppvt", "private_school", "title1", "free_reduced_lunch", "school_id" )]) ) form_lmer <- score_ppvt ~ private_school + title1 + free_reduced_lunch + distracted_ppvt + distracted_a1 + free_reduced_lunch:distracted_a1 + (1 + distracted_ppvt + distracted_a1 || school_id) ## Default: design matrices + reference lmer + identifiability design <- model_setup(form_lmer, data = dat) print(design) ## Full-rank filter (same as Ex_13b before prior setup / sampling) full_rank <- names(design$groupef.rank)[design$groupef.rank] dat_fr <- subset(dat, school_id %in% full_rank) dat_fr$school_id <- droplevels(dat_fr$school_id) design_fr <- model_setup(form_lmer, data = dat_fr) stopifnot(all(design_fr$groupef.rank)) print(design_fr) ## Variation: design only (no reference merMod fit) design_only <- model_setup(form_lmer, data = dat_fr, fit_mer = FALSE) stopifnot(is.null(design_only$lmer)) names(design_only$W) }
Prior family objects provide a convenient way to specify the details of the priors
used by matrix-input samplers such as rglmb and rlmb. See the documentation for
rglmb and rlmb for the details of how such model fitting
takes place.
Under a Beta(shape1, shape2) prior on the binomial probability
and a Binomial(, ) likelihood with identity link (
directly), the posterior is:
dNormal(mu, Sigma, dispersion = NULL) dNormal_Gamma(mu, Sigma_0, shape, rate) dIndependent_Normal_Gamma( mu, Sigma, shape, rate, max_disp_perc = 0.99, disp_lower = NULL, disp_upper = NULL ) dGamma( shape, rate, beta, Inv_Dispersion = TRUE, lik_shape = 1, max_disp_perc = 0.99, disp_lower = NULL, disp_upper = NULL )dNormal(mu, Sigma, dispersion = NULL) dNormal_Gamma(mu, Sigma_0, shape, rate) dIndependent_Normal_Gamma( mu, Sigma, shape, rate, max_disp_perc = 0.99, disp_lower = NULL, disp_upper = NULL ) dGamma( shape, rate, beta, Inv_Dispersion = TRUE, lik_shape = 1, max_disp_perc = 0.99, disp_lower = NULL, disp_upper = NULL )
mu |
a prior mean vector for the the modeling coefficients used in several pfamilies |
Sigma |
a prior variance-covariance matrix for |
dispersion |
the dispersion to be assumed when it is not given a prior. Should be provided
when the Normal prior is for the |
Sigma_0 |
prior variance-covariance on the precision-weighted coefficient scale for
|
shape |
The prior shape parameter for the gamma piece (inverse dispersion / precision).
When taking defaults from |
rate |
The prior rate parameter paired with |
max_disp_perc |
Specifies the percentile used to truncate the posterior dispersion
distribution when constructing the envelope for accept-reject sampling. This determines
the lower and upper bounds for the dispersion ( |
disp_lower |
lower bound truncation for dispersion |
disp_upper |
upper bound truncation for dispersion |
beta |
Initial coefficient matrix (1 |
Inv_Dispersion |
Logical (default
|
lik_shape |
Known shape parameter |
pfamily is a generic with methods for fitted objects such as rglmb and
rlmb. The dNormal() prior is supported for all response families.
The gaussian() family additionally supports dNormal_Gamma(),
dIndependent_Normal_Gamma(), and dGamma() (precision prior).
Intercept-only models with an identity link support two closed-form conjugate priors:
dBeta() for binomial(link = "identity") and
dGamma(Inv_Dispersion = FALSE) for poisson(link = "identity") and
Gamma(link = "identity").
A pfamily object represents a structured prior specification for use in Bayesian generalized linear modeling.
Each constructor function (e.g., dNormal(), dGamma(), dNormal_Gamma(), dBeta()) returns an object of
class "pfamily" containing the prior parameters, supported likelihood families, compatible link functions,
and a simulation function for posterior sampling.
These priors are designed to integrate seamlessly with modeling functions such as rglmb() and rlmb() in the
glmbayes package, which consume the pfamily object to define the prior distribution over model parameters.
The pfamily() generic retrieves the embedded prior from a fitted model object, while print.pfamily() displays its structure.
prior_list, simfun, and pfun. The named list prior_list holds the hyperparameters for the chosen
prior family. When a model function draws from the posterior, it passes prior_list into the element
simfun (e.g., rNormal_reg, rGamma_reg) so the low-level sampler receives
one consistent list structure regardless of which constructor built the pfamily.
Each constructor also stores a matching pfun (see prior_simfuncs) for i.i.d.\ draws from
the prior via pfun(n, prior_list, params), mirroring how simfun is used for the posterior.
Prior_Setup and default hyperparameters. Prior_Setup() fits an auxiliary GLM and returns
default mu, Sigma / Sigma_0, dispersion, Gamma shape and rate, and
related fields aligned with the data and prior-weight (pwt) choices. Those values can be supplied as
arguments to the pfamily constructors when you want package-default priors on the same scale as the model
matrix. Recommended use of shape and rate is not identical across constructors: for
dIndependent_Normal_Gamma(), pass shape = ps$shape_ING from Prior_Setup (not the
scalar ps$shape used by dNormal_Gamma()). For dGamma() with fixed coefficients
(beta), pass rate = ps$rate_gamma when that field is present (otherwise ps$rate); see
Prior_Setup and compute_gaussian_prior.
dNormal(): Specifies a multivariate normal prior over regression coefficients. It is conjugate for
Gaussian likelihoods with an identity link function, and serves as the primary implemented prior for all
other supported likelihood families in the current framework. This structure facilitates efficient posterior
sampling and analytical tractability. The returned prior_list includes ddef: TRUE when
dispersion was omitted or NULL (so the default 1 was used), FALSE when
dispersion was supplied explicitly (including 1).
For models with log-concave likelihood functions-such as Poisson, Binomial, and Gamma families-
posterior sampling under a dNormal prior is performed using a (Nygren and Nygren 2006)
likelihood subgradient approach. This method constructs tight enveloping functions around the posterior
using subgradients of the log-likelihood, enabling efficient accept-reject sampling even in high dimensions.
When the posterior distribution is approximately normal (typically the case for large sample sizes), the
area under the enveloping function is bounded above by a constant factor-approximately
in the univariate case, and in -dimensional models. These bounds ensure that
the rejection rate remains manageable and that the sampler remains computationally efficient.
The concept of conjugate priors was first formalized by (Raiffa and Schlaifer 1961), and further developed for regression models using g-prior structures by (Zellner 1986).
dGamma(): A Gamma prior with two distinct roles controlled by Inv_Dispersion:
Inv_Dispersion = TRUE (default): prior on the inverse dispersion (precision
or shape ). Used for dispersion estimation in Gaussian and
Gamma(log) models, typically in a Gibbs step with beta held fixed
(Gelman et al. 2013; Dobson 1990; McCullagh and Nelder 1989).
With Gaussian Prior_Setup output, prefer rate_gamma for rate
(see Details above).
Inv_Dispersion = FALSE: conjugate Gamma prior on the rate parameter
directly. Supports intercept-only models with an identity link:
Poisson (Gamma–Poisson conjugacy) and Gamma (Gamma–Gamma conjugacy).
Posterior draws are closed-form IID samples via rGamma_Conjugate_reg.
The lik_shape argument specifies the known Gamma likelihood shape (default 1,
i.e.\ exponential). Prior_Setup returns calibrated conj_poisson
hyperparameters for this path.
dBeta(): A Beta prior on the binomial probability for intercept-only
binomial(link = "identity") models. The posterior is a closed-form Beta draw
(Beta–Binomial conjugacy) produced by rBeta_reg. Arguments shape1
and shape2 are the prior pseudo-success and pseudo-failure counts.
Prior_Setup returns calibrated conj_beta hyperparameters for this path.
dNormal_Gamma(): Combines a multivariate normal prior on coefficients with a gamma prior on precision,
forming a conjugate structure for Gaussian models with unknown variance. The second argument is Sigma_0
(precision-weighted scale); it is aliased internally to Sigma in prior_list.
This formulation parallels classical Normal-Gamma models and is compatible with hierarchical extensions
(Gelman et al. 2013; Raiffa and Schlaifer 1961).
dIndependent_Normal_Gamma(): Similar to dNormal_Gamma(), but assumes independence between the
coefficient and precision priors. This structure is useful for models where prior independence is desired
or analytically convenient. With Prior_Setup on a Gaussian model, pass shape_ING as
the shape argument (see Details above).
Each pfamily object includes:
pfamily, prior_list, okfamilies, plinks, simfun, and pfun (see Value).
mu / Sigma: the surrogate Normal mean is shape1/(shape1+shape2) and
the surrogate variance is the Beta variance
shape1*shape2/((shape1+shape2)^2*(shape1+shape2+1)).
An object of class "pfamily" (with a concise print method). A list with elements:
pfamily |
Character string: the constructor name ( |
prior_list |
Named list of prior hyperparameters. It is passed into
|
okfamilies |
Character vector of implemented |
plinks |
Function of one |
simfun |
Function used to generate posterior draws (e.g., |
pfun |
Prior-simulation function paired with this constructor (e.g.,
|
Dobson A~J (1990).
An Introduction to Generalized Linear Models.
Chapman and Hall, London.
Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, Rubin DB (2013).
Bayesian Data Analysis, 3rd edition.
CRC Press.
McCullagh P, Nelder J~A (1989).
Generalized Linear Models.
Chapman and Hall, London.
Nygren K~N, Nygren L~M (2006).
“Likelihood Subgradient Densities.”
Journal of the American Statistical Association, 101(475), 1144–1156.
doi:10.1198/016214506000000357.
Raiffa H, Schlaifer R (1961).
Applied Statistical Decision Theory.
Clinton Press, Inc., Boston.
Zellner A (1986).
“On Assessing Prior Distributions and Bayesian Regression Analysis with g‐Prior Distributions.”
In Goel P~K, Zellner A (eds.), Bayesian Inference and Decision Techniques: Essays in Honor of Bruno de Finetti, volume 6 of Studies in Bayesian Econometrics and Statistics, 233–243.
Elsevier.
rglmb, rlmb for modeling functions that consume pfamily objects.
rNormal_reg, rNormalGamma_reg, rGamma_reg, rGamma_Conjugate_reg, rindepNormalGamma_reg for lower-level posterior sampling functions used by pfamily constructors.
prior_simfuncs for the paired prior-simulation functions (pfun).
Prior_Setup, Prior_Check for initializing and validating prior specifications.
EnvelopeBuild for envelope construction methods used in likelihood subgradient sampling (Nygren and Nygren 2006).
See also (Hastie and Pregibon 1992) for the original S modeling framework that inspired the design of pfamily.
Generic for converting a prior-setup object into a named list of
pfamily objects (class "pfamily_list") for block or
mixed-model samplers.
pfamily_list(object, ptypes = NULL, ...)pfamily_list(object, ptypes = NULL, ...)
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 |
... |
Passed to methods (unused for |
For a Prior_Setup_GLMM object, each group-effect
coefficient uses parameters from
object$pop.prior_list[[k]]:
"dNormal": dNormal(mu, Sigma, dispersion) with
group-effect variance 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
from object$pop.ing_prior[[k]] (calibrated by
Prior_Setup_GLMM via pop.dispersion.pwt /
pop.dispersion.nprior, mean-matched at that same marginal
). With effective prior sample size
and population predictors for component
,
Requires (equivalently
pop.dispersion.pwt ). 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 priors: those
use dGamma_list (per-group dispformula) or a
single dGamma built from
object$group.ing_prior when dispersion is pooled.
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.
Prior_SetupGroup, Prior_Setup_GLMM,
pfamily,
dNormal, dNormal_Gamma,
dIndependent_Normal_Gamma
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
decays like (linear in l,
vs. Claim~3's for variance). Each coefficient's per-sweep
deviation from its reference mean is standardized by that same sweep's
cross-chain SD, , so coefficients on very different scales can be overlaid
on one chart and Z(l) is expected to shrink toward 0 as l
grows.
plot_mean_convergence(hist, ...)plot_mean_convergence(hist, ...)
hist |
Object of class |
... |
Passed to methods; see |
The x-axis is the inner Gibbs sweep index , and each point is a
cross-chain statistic computed over the n_chains
independent replicate chains, all of which sit at sweep
simultaneously. Nothing here is a running statistic along a single
chain.
The plotted quantity is the standardized deviation
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
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.
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.
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.
hist invisibly.
plot_var_convergence,
lmerb_posterior_mean
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 satisfies
for a convergent (symmetric, eigenvalues in ) matrix .
Var(l) / Var_ref ratios are therefore bounded above by 1 and
increase toward it as l grows.
plot_var_convergence(hist, ...)plot_var_convergence(hist, ...)
hist |
Object of class |
... |
Passed to methods (e.g. |
The x-axis is the inner Gibbs sweep index , 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 — the variance across the n_chains independent
replicate chains, all of which are at sweep 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
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
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.
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.
sweep_history is recorded only when all chains advance together,
one sweep at a time, so that the cross-chain spread at sweep is
meaningful. Two situations produce no history:
The exact route (fixed , 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.
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.
hist invisibly.
plot_mean_convergence,
lmerb_posterior_covariance
Print random effects in lme4 layout
## S3 method for class 'ranef.lmerb' print(x, ...)## S3 method for class 'ranef.lmerb' print(x, ...)
x |
An object of class |
... |
Ignored. |
x, invisibly.
Helper function to facilitate the Setup of Prior Distributions for glm models.
Prior_Setup( formula, family = gaussian(), data = NULL, weights = NULL, subset = NULL, na.action = na.fail, offset = NULL, contrasts = NULL, pwt = NULL, pwt_default_low = 0.01, pwt_default_high = 0.05, n_prior = NULL, sd = NULL, dispersion = NULL, intercept_source = c("null_model", "full_model"), effects_source = c("null_effects", "full_model"), mu = NULL, k = 1, ... )Prior_Setup( formula, family = gaussian(), data = NULL, weights = NULL, subset = NULL, na.action = na.fail, offset = NULL, contrasts = NULL, pwt = NULL, pwt_default_low = 0.01, pwt_default_high = 0.05, n_prior = NULL, sd = NULL, dispersion = NULL, intercept_source = c("null_model", "full_model"), effects_source = c("null_effects", "full_model"), mu = NULL, k = 1, ... )
formula |
an object of class |
family |
a description of the error distribution and link function to be used in the model. |
data |
an optional data frame, list or environment (or object
coercible by |
weights |
an optional vector of ‘prior weights’ to be used
in the fitting process. Should be |
subset |
an optional vector specifying a subset of observations to be used in the fitting process. (See additional details about how this argument interacts with data-dependent bases in the ‘Details’ below.) |
na.action |
how |
offset |
this can be used to specify an a priori known
component to be included in the linear predictor during fitting.
This should be |
contrasts |
an optional list. See the |
pwt |
Weight on the prior relative to the likelihood function at the maximum likelihood
estimate. If supplied, this value is used directly (scalar or one value per coefficient).
If |
pwt_default_low |
Default prior weight used when |
pwt_default_high |
Default prior weight used when |
n_prior |
Optional scalar effective prior sample size (on the |
sd |
Optional vector argument with the prior standard deviations for the coefficients |
dispersion |
Optional scalar dispersion override (default |
intercept_source |
Specifies the method through which the prior mean for the intercept term is set. Options are based on the null intercept only model (null_model) or full_models. The default is the null model which is safer if variables are not centered. |
effects_source |
Specifies the method through which the prior means for the effects terms are set. Options are null_effects (prior means set to zero) or full_model (effect means set to match maximum likelihood estimates). |
mu |
Optional vector argument with the prior means for the coefficients |
k |
Scalar (default |
... |
For For |
Inputs to the function
The inputs to Prior_Setup() fall into three conceptual categories:
1. Model specification
formula: structure of the GLM (response and predictors).
family: error distribution and link.
data, weights, subset, na.action,
offset, contrasts, control, ...: as in
glm.
2. Prior variance–covariance specification
pwt: prior weight relative to the likelihood. If scalar, used to
construct a Zellner-type g-prior. If vector, applied elementwise.
n_prior: optional scalar effective prior sample size. Replaces
scalar pwt only when pwt is scalar and sd is not
used; otherwise supplies precision-prior / calibration only.
sd: optional vector of prior standard deviations. If provided,
used to compute pwt from the diagonal of vcov(glm_full).
pwt_default_low, pwt_default_high: defaults for pwt
when not supplied.
3. Prior mean specification
intercept_source: method for setting the prior mean of the intercept
("null_model" or "full_model").
effects_source: method for setting the prior mean of the effects
("null_effects" or "full_model").
mu: optional user-specified prior mean vector; overrides other
centering logic if provided.
Prior covariance and Zellner scaling
Let be the covariance matrix of the
full-model GLM coefficients. For non-Gaussian families, the prior covariance
is:
where denotes elementwise multiplication.
Intercept-only Poisson(link = "identity") conjugate prior on the rate
When the design is a single column (intercept only), family = poisson(),
link = "identity", scalar pwt, and offsets are zero, the effective
conjugate prior observation count already satisfies
(otherwise conj_poisson remains NULL with a warning). Writing the
weighted mean , the output list
component conj_poisson stores
and , so the prior mean for the rate matches .
Omitting optional mu resets the surrogate Normal summaries mu and the
sole diagonal element of Sigma to these Gamma moments (Sigma_{11} =
).
For Gaussian families, Prior_Setup() also constructs the dispersion-free
covariance
which under scalar pwt and the default calibration reduces to
Gaussian Normal–Gamma calibration and
For family = gaussian(), the function performs the Normal–Gamma
calibration described in (Nygren 2025). Let:
,
,
the weighted least-squares estimator,
the dispersion-free prior covariance.
The marginal quadratic term is
where is the weighted residual sum of squares at
. Under the default scalar-pwt Zellner mapping
,
this simplifies to
which makes the limiting behavior as transparent.
The calibrated dispersion is
and the Normal–Gamma hyperparameters are
The independent Normal–Gamma shape is
Posterior summaries for the conjugate Normal–Gamma prior
Under the conjugate Normal–Gamma prior (used by dNormal_Gamma()),
the posterior has:
Posterior mean
Posterior expectation of
Posterior covariance
Weak-prior limits (Theorems 2 and 3)
As (equivalently ),
, and the conjugate Normal–Gamma
posterior converges to the classical weighted least-squares limit:
For the independent Normal–Gamma prior used by
dIndependent_Normal_Gamma(), neither the posterior mean nor the
posterior covariance is available in closed form; the posterior must be
obtained by numerical integration or sampling (e.g.,
rindepNormalGamma_reg()). Theorem 3 in
(Nygren 2025) shows that the ING posterior has
the same weak-prior limit as the conjugate Normal–Gamma posterior:
A list of class "PriorSetup" with components:
mu |
Prior mean vector (length equal to the number of coefficients). |
Sigma |
Coefficient-scale prior variance–covariance matrix. |
Sigma_0 |
For |
dispersion |
Calibrated dispersion (Gaussian models only), equal to
|
shape |
Derived prior Gamma shape parameter for the Normal–Gamma prior
on precision (Gaussian only), |
shape_ING |
For |
rate |
Derived prior Gamma rate parameter (Gaussian only), using the
calibrated |
rate_gamma |
For |
coefficients |
Named numeric vector of returned coefficient values.
For |
model |
The model frame used to construct the design matrix (if
|
x |
The model matrix used (if |
y |
The response vector used (if |
call |
The matched call to |
PriorSettings |
A list containing prior configuration details, including
|
conj_poisson |
|
Nygren K (2025). “Core 10: Technical Derivations for Priors Returned by Prior_ Setup.” Vignette in the glmbayesCore R package. R vignette name: Core-10; formerly glmbayes appendix Chapter-A12.
pfamily for prior-family objects and the constructors
dNormal, dNormal_Gamma, dGamma,
and dIndependent_Normal_Gamma.
rglmb, rlmb for matrix-input fits with a
pfamily built from Prior_Setup() output; rglmb,
rlmb for matrix-based sampling that consumes the same prior
structure; simfuncs for functions that take a prior_list
assembled from those components (including
rindepNormalGamma_reg for
dIndependent_Normal_Gamma()).
multi_prior_setup for a matrix/cbind response with Gaussian;
use with rlmb
Prior_Setup per column.
(Zellner 1986); (Raiffa and Schlaifer 1961); (Gelman et al. 2013); (McCullagh and Nelder 1989); (Nygren 2025); (Nygren 2025).
Other prior:
Prior_Check(),
multi_prior_setup()
## Not run: ## Full runnable examples are maintained in \pkg{glmbayesCore}: example(Prior_Setup, package = "glmbayesCore", ask = FALSE, echo = TRUE) ## End(Not run)## Not run: ## Full runnable examples are maintained in \pkg{glmbayesCore}: example(Prior_Setup, package = "glmbayesCore", ask = FALSE, echo = TRUE) ## End(Not run)
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:
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 )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 )
formula |
Mixed-model formula passed to |
data |
Data frame containing all variables in |
family |
Model |
REML |
Logical scalar passed to |
control |
Optional |
start |
Optional starting values passed to |
verbose |
Integer scalar passed to |
subset, na.action
|
Optional arguments passed through to
|
weights, offset
|
Optional arguments passed through to
Not fully active for mixed-model sampling yet. They affect the
classical reference fit used for prior calibration. Low-level
|
contrasts |
Optional contrast list passed to
|
dispformula |
One-sided formula selecting the Block~1 measurement-
dispersion structure: |
pop.pwt |
Prior weight(s) in |
pop.pwt_default_low |
Scalar in |
pop.pwt_default_high |
Scalar in |
pop.nprior |
Optional absolute effective prior sample
size(s), in group units ( |
pop.sd |
Optional prior standard deviation(s), on the coefficient
scale of |
pop.dispersion.pwt |
Optional relative prior weight(s) in
|
pop.dispersion.nprior |
Optional absolute effective prior
sample size(s) (in group units) for the Block~2 |
pop.max_disp_perc |
Scalar in |
pop.intercept_source |
Character string controlling the prior mean for
the global intercept hyperparameter |
pop.effects_source |
Character string controlling the prior mean for
all other Block~2 hyperparameters (including population-mean slopes
stored as |
pop.mu |
Optional prior mean override(s) for the population-level
(Block~2, |
group.dispersion.pwt |
Optional relative prior weight(s) in |
group.dispersion.nprior |
Optional positive scalar: absolute effective
prior sample size for the Block~1 |
group.alpha_target |
Scalar in |
group.max_disp_perc |
Scalar in |
Block 1 (per-group, independent):
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 , independent):
Naming convention. Top-level fields use pop.* for Block~2
(population / , ) and group.* for Block~1
(within-group / observation , RE covariance ).
Nested entries in pop.prior_list[[k]] use glmbayes-style
mu, Sigma, and dispersion (where dispersion
is the plug-in RE variance ).
Why default calibration depends on classical estimates.
Prior_Setup_GLMM scales Block~2 covariances from
vcov(fit_ref) by 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:
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.
Every X_hyper[[k]] column maps to a fixef(fit_ref) term.
Each RE variance from the reference fit is strictly positive.
Object of class "Prior_Setup_GLMM" with fields:
formulaModel formula.
familyFamily object.
pop.pwtPrior 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.npriorNamed list with one named numeric vector per
random-effect component – the effective prior sample size (group
units) actually implied, , 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.sdNamed list with one named numeric vector per
random-effect component – the effective coefficient-scale prior
standard deviations actually used
(), 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.muNamed 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.pwtNamed per-component vector of relative
dispersion prior weights (always present; consistent with
pop.dispersion.nprior via ).
pop.dispersion.npriorNamed 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.pwtGaussian models only: scalar or length-
vector of relative prior weights for Block~1 – a
scalar unless a per-group vector was supplied or
group.alpha_target calibration was active, in which case it
is the resolved length- vector (see dGamma_list()).
group.dispersion.npriorGaussian models only: effective prior
sample size for Block~1 , on the same pooled-vs.-vector
basis as group.dispersion.pwt (a scalar on the observation scale, or a
named length- vector of per-group
reflecting the group.alpha_target-calibrated group.dispersion.pwt
when calibration was active; see group.pwt_calibration).
group.alpha_targetThe value used, or NULL
if calibration was not active.
group.pwt_calibrationNULL 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_percGaussian models only: scalar or
length- vector of the tail probability used for the Block~1
truncation window (group.ing_prior) – the
scalar as supplied, or the resolved length- vector when a
vector was supplied.
pop.max_disp_percNamed length-
vector of resolved per-component values used for the Block~2
truncation windows in pop.ing_prior (always
present).
pop.intercept_sourceHow the global intercept prior mean
was chosen ("null_model" or "full_model").
pop.effects_sourceHow other Block~2 prior means were
chosen ("null_effects" or "full_model").
block_formulaWithin-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_sdNamed vector sqrt(Psi) from the reference
fit; shared population RE standard deviations for per-group calibration.
dataData frame passed to Prior_Setup_GLMM()
(reference for dGamma_list() diagnostics).
designFull model_setup object (all groups),
including length-n weights and offset vectors
aligned with design$y.
mer_fitReference 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_refThe calibration reference for Block~2
(fixef/) 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.fitNULL when dispformula = ~1;
otherwise the same glmmTMB object as fit_ref.
group.dispersion.refNULL when dispformula = ~1;
otherwise a named length- 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.dispersionFor Gaussian models: resolved Block~1
observation used to mean-match the Gamma prior.
When dispformula = ~1, the pooled A12 §3.3.4 aggregate
(Part VI
Omega_j included); classical lmer residual in
attr(., "classical"). When dispformula is per-group,
a named length- 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- routes.
group.SigmaDiagonal RE covariance matrix (Block~1),
classical reference-fit diagonals.
pop.dispersion.refNamed vector of classical reference-fit
(VarCorr / Psi) before the A12 §3.3.4
hyper-regression remargin; compare to
pop.prior_list[[k]]$dispersion.
pop.prior_listNamed Block~2 prior list per RE coefficient.
For Gaussian models dispersion is the A12 §3.3.4 marginal
hyper-regression center for (not the raw VarCorr
plug-in; see pop.dispersion.ref).
pop.ing_priorNamed per-component list of the prospective
dIndependent_Normal_Gamma calibration: Gamma precision-prior
shape and rate
(the lmebayesCore default
calibration with pop.dispersion.nprior; since
rate , the implied
inverse-Gamma prior on has mean exactly
), where is the A12 §3.3.4
marginal hyper-regression center in
pop.prior_list[[k]]$dispersion, and the default
truncation window disp_lower / disp_upper: the
/
quantiles of the limiting posterior
– the weak-prior
() limit of the Block~2 posterior Gamma for the
precision (lmebayesCore Chapter A12, Theorem 2; inverted to a
interval). This window is identical for all
, covers
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_priorGaussian 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 shared across all group
levels) – mean-matched shape / rate with
,
, and =
group.dispersion (pooled A12 §3.3.4 aggregate of
per-group by default; same ING algebra
as pop.ing_prior for ; requires
group.dispersion.pwt ),
plus disp_lower / disp_upper as the central
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 dispersionA 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 /
quantiles of Gamma(shape_ING + n_j/2, rate_post), where
rate_post is sigma2_hat mean-matched at that same
inflated shape ( added only for the window, so it
reflects the posterior spread the sampler's own envelope
machinery actually draws 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); 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.
dispformulaThe dispformula supplied.
x invisibly.
model_setup, Prior_Setup,
build_mu_all,
pfamily_list, dGamma_list
## Prior_Setup_GLMM() — Block~1 / Block~2 prior calibration (Ex_13b data) ## ## Same formula and cleaning as ## demo("Ex_13b_rLMM_estimated_dispersion_known_vcov_BigWordClub_PartVI"). ## Requires suggested package bayesrules. if (requireNamespace("bayesrules", quietly = TRUE)) { data(big_word_club, package = "bayesrules") dat <- big_word_club dat$school_id <- factor(dat$school_id) dat <- subset( dat, !is.na(score_ppvt) & !is.na(invalid_ppvt) & invalid_ppvt == 0L & complete.cases(dat[, c( "score_ppvt", "distracted_a1", "distracted_ppvt", "private_school", "title1", "free_reduced_lunch", "school_id" )]) ) form_lmer <- score_ppvt ~ private_school + title1 + free_reduced_lunch + distracted_ppvt + distracted_a1 + free_reduced_lunch:distracted_a1 + (1 + distracted_ppvt + distracted_a1 || school_id) ## Keep full-rank schools (sampler / per-group ING path) design0 <- model_setup(form_lmer, data = dat) dat <- subset(dat, school_id %in% names(design0$groupef.rank)[design0$groupef.rank]) dat$school_id <- droplevels(dat$school_id) ## Pooled measurement dispersion (dispformula = ~1, default) ps_pooled <- Prior_Setup_GLMM( form_lmer, data = dat, pop.pwt = 0.01 ) print(ps_pooled) names(ps_pooled$pop.prior_list) ## Per-group Block~1 calibration (Ex_13b-style dispformula). ## group.alpha_target = NULL skips the ellipsoid pwt search (faster for examples); ## Ex_13b uses group.alpha_target = 0.01 with group.dispersion.pwt = 0.1 as floor. ps_group <- Prior_Setup_GLMM( form_lmer, data = dat, pop.pwt = 0.01, dispformula = ~school_id, group.max_disp_perc = 0.8, group.dispersion.pwt = 0.1, group.alpha_target = NULL ) print(ps_group) stopifnot(is.list(ps_group$group.ing_prior[[1L]])) length(ps_group$group.ing_prior) }## Prior_Setup_GLMM() — Block~1 / Block~2 prior calibration (Ex_13b data) ## ## Same formula and cleaning as ## demo("Ex_13b_rLMM_estimated_dispersion_known_vcov_BigWordClub_PartVI"). ## Requires suggested package bayesrules. if (requireNamespace("bayesrules", quietly = TRUE)) { data(big_word_club, package = "bayesrules") dat <- big_word_club dat$school_id <- factor(dat$school_id) dat <- subset( dat, !is.na(score_ppvt) & !is.na(invalid_ppvt) & invalid_ppvt == 0L & complete.cases(dat[, c( "score_ppvt", "distracted_a1", "distracted_ppvt", "private_school", "title1", "free_reduced_lunch", "school_id" )]) ) form_lmer <- score_ppvt ~ private_school + title1 + free_reduced_lunch + distracted_ppvt + distracted_a1 + free_reduced_lunch:distracted_a1 + (1 + distracted_ppvt + distracted_a1 || school_id) ## Keep full-rank schools (sampler / per-group ING path) design0 <- model_setup(form_lmer, data = dat) dat <- subset(dat, school_id %in% names(design0$groupef.rank)[design0$groupef.rank]) dat$school_id <- droplevels(dat$school_id) ## Pooled measurement dispersion (dispformula = ~1, default) ps_pooled <- Prior_Setup_GLMM( form_lmer, data = dat, pop.pwt = 0.01 ) print(ps_pooled) names(ps_pooled$pop.prior_list) ## Per-group Block~1 calibration (Ex_13b-style dispformula). ## group.alpha_target = NULL skips the ellipsoid pwt search (faster for examples); ## Ex_13b uses group.alpha_target = 0.01 with group.dispersion.pwt = 0.1 as floor. ps_group <- Prior_Setup_GLMM( form_lmer, data = dat, pop.pwt = 0.01, dispformula = ~school_id, group.max_disp_perc = 0.8, group.dispersion.pwt = 0.1, group.alpha_target = NULL ) print(ps_group) stopifnot(is.list(ps_group$group.ing_prior[[1L]])) length(ps_group$group.ing_prior) }
Runs Prior_Setup on each block subset of the data.
Typical callers are rNormal_reg_group /
rNormalGLM_reg_group or the lmbBlock / glmbBlock
wrappers in lmebayes.
Prior_SetupGroup( formula, group, family = gaussian(), data = NULL, weights = NULL, subset = NULL, na.action = na.fail, offset = NULL, contrasts = NULL, pwt = NULL, pwt_default_low = 0.01, pwt_default_high = 0.05, n_prior = NULL, sd = NULL, dispersion = NULL, intercept_source = c("null_model", "full_model"), effects_source = c("null_effects", "full_model"), mu = NULL, k = 1, ... )Prior_SetupGroup( formula, group, family = gaussian(), data = NULL, weights = NULL, subset = NULL, na.action = na.fail, offset = NULL, contrasts = NULL, pwt = NULL, pwt_default_low = 0.01, pwt_default_high = 0.05, n_prior = NULL, sd = NULL, dispersion = NULL, intercept_source = c("null_model", "full_model"), effects_source = c("null_effects", "full_model"), mu = NULL, k = 1, ... )
formula |
A |
group |
Group partition: |
family |
a description of the error distribution and link function to be used in the model. |
data |
an optional data frame, list or environment (or object
coercible by |
weights |
an optional vector of ‘prior weights’ to be used
in the fitting process. Should be |
subset |
an optional vector specifying a subset of observations to be used in the fitting process. (See additional details about how this argument interacts with data-dependent bases in the ‘Details’ below.) |
na.action |
how |
offset |
this can be used to specify an a priori known
component to be included in the linear predictor during fitting.
This should be |
contrasts |
an optional list. See the |
pwt, n_prior, sd, mu, dispersion, k
|
Each of these six calibration
arguments (passed to
The two forms may be mixed freely across the six arguments in one call. |
pwt_default_low |
Default prior weight used when |
pwt_default_high |
Default prior weight used when |
intercept_source |
Specifies the method through which the prior mean for the intercept term is set. Options are based on the null intercept only model (null_model) or full_models. The default is the null model which is safer if variables are not centered. |
effects_source |
Specifies the method through which the prior means for the effects terms are set. Options are null_effects (prior means set to zero) or full_model (effect means set to match maximum likelihood estimates). |
... |
Not used. |
Prior_SetupGroup() calls Prior_Setup once per block,
on that block's rows only, and returns the results as a named list keyed by
block ID.
Prior_Setup_GLMM calibrates a
hierarchical prior (one reference lmer/glmer fit,
Block~2 hyperpriors, cross-group shrinkage). Prior_SetupGroup()
calibrates independent priors per row block with no shared
or ; row-block engines draw blocks iid.
Names for per-block calibration lists must match normalize_group
block IDs exactly.
After calibration, convert to sampler-ready pfamily objects with
pfamily_list(x, ptypes = ...) (see Prior families via
pfamily_list() below for ptypes and allowed families).
A named list of class "Prior_SetupGroup". Each element is a
Prior_Setup result for one block.
x invisibly.
An object of class "pfamily_list": a named list of
pfamily objects with attribute "ptypes".
Prior_Setup, multi_prior_setup,
normalize_group, pfamily_list,
rNormal_reg_group, rNormalGLM_reg_group,
Prior_Setup_GLMM
Methods mirroring ranef for lmerb and
glmerb fits. Values summarize Block~1 random-effects draws
stored in fit$groupef (posterior chain mean by default,
or the exact ICM fit$groupef.mode).
## S3 method for class 'lmerb' ranef(object, condVar = FALSE, drop = FALSE, type = c("mean", "mode"), ...) ## S3 method for class 'glmerb' ranef(object, condVar = FALSE, drop = FALSE, type = c("mean", "mode"), ...)## S3 method for class 'lmerb' ranef(object, condVar = FALSE, drop = FALSE, type = c("mean", "mode"), ...) ## S3 method for class 'glmerb' ranef(object, condVar = FALSE, drop = FALSE, type = c("mean", "mode"), ...)
object |
An |
condVar |
Logical; if |
drop |
Logical; if |
type |
|
... |
Ignored. |
An object of class "ranef.lmerb": a named list with one
data.frame per grouping factor (rows = levels, columns = RE terms).
Attributes include type, n_draws, and optionally
postVar.
ranef, coef.lmerb,
fixef.lmerb, lmerb.
rglmerb draws from the posterior of a Bayesian generalized linear
mixed model — Poisson, binomial, Gamma, or Gaussian — with a single
grouping factor and uncorrelated random effects, given a design from
model_setup and priors for the population level
(pfamily_list) and, for Gaussian responses only, the residual
variance (dispprior_list).
rglmerb( n, design, pfamily_list, family = poisson(), dispprior_list = NULL, offset = NULL, weights = 1, gap_tol = 0.0196, tv_tol = 0.01, mode_gap_max = 1, collect_block1 = TRUE, verbose = TRUE, progbar = FALSE, sim_method = "DEFAULT", simulate = TRUE )rglmerb( n, design, pfamily_list, family = poisson(), dispprior_list = NULL, offset = NULL, weights = 1, gap_tol = 0.0196, tv_tol = 0.01, mode_gap_max = 1, collect_block1 = TRUE, verbose = TRUE, progbar = FALSE, sim_method = "DEFAULT", simulate = TRUE )
n |
Integer. Number of independent chains in the main stage.
Required when |
design |
A |
pfamily_list |
Named list of Block~2 |
family |
A |
dispprior_list |
Block~1 observation-dispersion prior, in the same
shapes as the |
offset, weights
|
Observation offset and prior weights, as in
|
gap_tol |
Legacy mode–mean gap for deriving the pilot chain count when
|
tv_tol |
Total variation tolerance for convergence calibration.
Inner Gibbs sweeps and pilot chain counts are derived internally.
Ignored when |
mode_gap_max |
Pilot inner-sweep calibration (non-Gaussian and
Gaussian+ING only). Ignored when |
collect_block1 |
Collect |
verbose |
Print stage headers and diagnostics. |
progbar |
Progress bars when |
sim_method |
Simulation method for |
simulate |
Logical (default |
The returned sample is n mutually independent draws, one per
replicate chain, taken at that chain's final sweep. It is analyzed like an
iid sample rather than like an MCMC trace.
This is the matrix-level interface; for a formula interface use
glmerb() in lmebayes. rlmerb is the
Gaussian-only counterpart.
rglmerb estimates the same two-stage structure as
rlmerb, with a link function between the linear predictor
and the response mean. For group with link :
The across-group stage is unchanged and remains Gaussian: group
coefficients scatter around a population expectation
with diagonal covariance
. What changes is the within-group
stage and, with it, the interpretation of , which now lives
on the link scale — log rates for poisson(), log-odds for
binomial(). Summaries that are meaningful on the response scale
need the inverse link applied to the draws.
Supported families are gaussian(), poisson(),
binomial(), and Gamma() with their standard links. The
restriction is log-concavity of the likelihood, which the sampler
requires (Nygren and Nygren 2006).
designA model_setup object built with the
matching family. For non-Gaussian responses its estimability
report matters more than usual: a group with no variation in the
response (all-zero counts, a single binomial outcome level) has no
finite maximum likelihood estimate, and is reported separately from
the rank check. See
vignette("Chapter-B08", package = "lmebayesCore").
pfamily_listThe population-level prior, one component per
random-effect coefficient, declaring whether each between-group
variance is known (dNormal) or estimated
(dIndependent_Normal_Gamma).
dispprior_listOnly meaningful for gaussian(),
where it gives the residual variance. poisson() and
binomial() have no free dispersion parameter, so it must be
NULL.
As in rlmerb, rglmerb runs n independent
replicate chains and stores one draw from each, taken at the final sweep,
so the n draws are mutually independent and ordinary Monte Carlo
standard errors apply. There is nothing to thin and no autocorrelation to
diagnose.
The strength of the accompanying guarantee is weaker here than in the
Gaussian case, and the difference is worth understanding rather than
glossing over. The calibration behind m_convergence is derived for
Gaussian targets. A Poisson or binomial posterior is not Gaussian, so the
engine applies the same calibration at curvature evaluated near the
posterior mode. The sweep count is therefore a well-motivated
budget rather than a certificate, and tv_tol should be
read in that spirit. Only family = gaussian() with fixed
and all-dNormal components yields exact draws.
A second difference follows from the same fact. Because the posterior mode
is not the posterior mean when the target is skewed, chains started at the
mode begin at an unknown distance from the target. rglmerb
therefore runs a short pilot stage first and uses its cross-chain
mean as the starting point for the main chains, which shortens the main
stage. Pilot draws are diagnostic only and are discarded — they are
never part of the returned sample.
vignette("Chapter-B03", package = "lmebayesCore") covers the
trade-off between pilot size and main-stage length.
simulate = FALSE returns the joint posterior mode by iterated
conditional modes (or the exact mean for gaussian()), at fixed
variance-component plug-ins. No sweeps run and no RNG is consumed, so the
result is deterministic — useful while iterating on model
specification. All draw, pilot, and convergence slots are NULL.
An object of class c("rglmerb", "list"), with the same
layout as rlmerb. Each draw slot holds one row per
replicate chain, taken at that chain's final sweep.
popefNamed list of population draws , one
n matrix per random-effect coefficient, on
the link scale. popef.means and popef.mode are the
draw mean and the posterior mode.
groupefLong-format data frame of group coefficients
(n rows: draw, the
grouping factor, then one column per coefficient), also on the link
scale. groupef.mode is the mode matrix.
Collected only when collect_block1 = TRUE.
popef.dispersionDraws of the between-group variances
; constant for dNormal components, varying for
estimated ones. popef.dispersion.mean is the column mean.
group.dispersionResidual variance, Gaussian responses
only; absent for poisson() and binomial().
m_convergence, convergence_info
Sweeps run per
chain and the calibration behind that count (method,
lambda_star, draw_engine,
sim_method_used, icm_info).
pilotPresent when a pilot stage ran, holding its size
n, the chisq calibration constant, and its
draws. Diagnostic only — these are not posterior draws.
Inputs are echoed back in design, family, n,
call, pfamily_list, dispprior_list,
prior.weights and offset, alongside the packed
prior and thin Prior summaries.
When simulate = FALSE, every draw, pilot and convergence slot is
NULL and only the mode fields are populated.
rlmerb, rglmerb_point,
rLMMNormal_reg_known_vcov,
rLMMNormal_reg_estimated_vcov,
rLMMindepNormalGamma_reg_known_vcov,
rLMMindepNormalGamma_reg_estimated_vcov, rGLMM_reg,
Prior_Setup_GLMM
rlmerb draws from the posterior of a Bayesian linear mixed model
with a single grouping factor and uncorrelated random effects, given a
design from model_setup and priors for the population level
(pfamily_list) and the residual variance (dispprior_list).
rlmerb( n, design, pfamily_list, dispprior_list, offset = NULL, weights = 1, tv_tol = 0.01, progbar = TRUE, verbose = TRUE, print_icm_table = TRUE, gap_tol = 0.0196, mode_gap_max = 1, diag_sweeps = FALSE, sim_method = "DEFAULT", simulate = TRUE )rlmerb( n, design, pfamily_list, dispprior_list, offset = NULL, weights = 1, tv_tol = 0.01, progbar = TRUE, verbose = TRUE, print_icm_table = TRUE, gap_tol = 0.0196, mode_gap_max = 1, diag_sweeps = FALSE, sim_method = "DEFAULT", simulate = TRUE )
n |
Integer. Number of stored draws (each draw is one full pass through
|
design |
A |
pfamily_list |
Named list of Block~2 |
dispprior_list |
Block~1 observation-dispersion prior, in the same
shapes accepted by the |
offset, weights
|
Observation offset and prior weights, as in
|
tv_tol |
Single numeric in |
progbar |
Logical. Show a text progress bar during sampling.
Default |
verbose |
Logical. Print the reference-vs-ICM table and the convergence
calibration line. Default |
print_icm_table |
Logical. When |
gap_tol |
Legacy mode–mean gap tolerance for the pilot stage when
any population component uses |
mode_gap_max |
Pilot inner-sweep calibration for ING population models
(default |
diag_sweeps |
Diagnostic flag for ING population models with a pilot stage.
When |
sim_method |
Simulation method: |
simulate |
Logical (default |
The returned sample is n mutually independent draws — one per
replicate chain, taken at that chain's final sweep, where the number of
sweeps is calibrated in advance so the final sweep is within
tv_tol of the target posterior. It is therefore analyzed like an
iid sample rather than like an MCMC trace.
This is the matrix-level interface. For a formula interface, use
lmerb() in lmebayes. The Poisson/binomial/Gamma counterpart
is rglmerb, and the non-hierarchical analogue is
rlmb.
rlmerb estimates the two-stage Gaussian mixed model in a
centered (Lindley–Smith) parameterization
(Lindley and Smith 1972). For group
:
Read the second line as the model for how groups differ: each group's own
coefficient vector is drawn around a population expectation
, with measuring how far groups
scatter from it. The estimand is therefore the pair
at the population level together with all
group vectors , and for the residual scale.
What distinguishes this from the lme4 presentation is that the
full coefficient appears in the likelihood, not a
mean-zero deviation. Marginalizing recovers the familiar form with
, , and
, so corresponds to
lme4's fixed effects and to its random effects. The
practical consequence is that groupef holds directly
(comparable to lme4::coef()), not .
is diagonal, , because model_setup accepts only
uncorrelated (||) random-effect terms and a single grouping factor.
There is no random-effect correlation parameter to estimate.
Three arguments carry the statistical content; the rest are controls.
designA model_setup object, which fixes
, , the grouping factor, and the level-2 design
. It also reports whether the data can identify
the model at all; see
vignette("Chapter-B08", package = "lmebayesCore").
pfamily_listThe prior on the population level, one
component per random-effect coefficient. This is where you declare
whether the between-group variance is
known (dNormal, treating the supplied dispersion as
fixed) or estimated (dIndependent_Normal_Gamma,
placing a Gamma prior on the precision ).
dispprior_listThe prior on the residual variance
: a fixed value (scalar, or one per group), a single
dGamma() to estimate one pooled
, or a named list of dGamma() objects to estimate
a separate per group.
Both priors are normally produced from one
Prior_Setup_GLMM object, via pfamily_list()
and either ps$group.dispersion or dGamma_list(ps).
Those two choices are also what determine which internal engine runs.
You never name an engine; it follows from the priors. The statistically
meaningful consequence is that when is fixed and every
component is dNormal, the posterior is exactly multivariate
normal and can be sampled directly; every other combination requires
iterative sampling. See
vignette("Chapter-B01", package = "lmebayesCore") for the full
mapping and vignette("Chapter-B04", package = "lmebayesCore")
for the exact route.
This is the part most likely to differ from other samplers, and it governs how the output should be analyzed.
rlmerb runs n independent replicate chains and
keeps one draw from each — the state after that chain's final
sweep. It does not run one long chain and thin it. Two consequences:
The n draws are mutually independent. They are a
genuine iid sample, not an autocorrelated series. Ordinary Monte Carlo
standard errors apply directly; there is no effective sample size to
compute, nothing to thin, and no trace plot or to
inspect on the returned sample.
Each draw comes from a distribution deliberately made
close to the target. The sweep count m_convergence is not a
burn-in guess. It is calibrated before sampling so that the
distribution reached at the final sweep is within tv_tol of the
exact posterior in total variation. tv_tol is a per-draw
guarantee, and it is reported on the fit.
So n and tv_tol control different errors and can be set
independently: n controls Monte Carlo noise, tv_tol controls
the residual bias of each draw. To certify the whole n-draw sample
jointly at level , pass tv_tol = alpha / n; because
cost grows only logarithmically in tv_tol, this typically
adds a few sweeps.
In the exactly-Gaussian case (fixed , all dNormal)
the draws are exact rather than approximate: m_convergence is
1 and no approximation is involved.
How m_convergence is derived, and what to do when a warning
reports strong cross-block coupling, is covered in
vignette("Chapter-B02", package = "lmebayesCore") and
vignette("Chapter-B03", package = "lmebayesCore").
simulate = FALSE returns only the posterior mean of
at fixed variance-component plug-ins, computed by
iterated conditional modes. No sweeps run and no RNG is consumed, so the
result is deterministic and fast — useful when iterating on model
specification. All draw and convergence slots are NULL. This same
quantity is where the replicate chains start when simulate = TRUE,
which is why it is always computed.
An object of class c("rlmerb", "list"). Each draw slot holds
one row (or one value) per replicate chain, taken at that chain's final
sweep; see the "What the returned draws are" section above.
Population level (, the analogue of lme4
fixed effects):
popefNamed list, one entry per random-effect
coefficient. Entry is an n matrix of
draws, columns named by colnames(design$W[[k]]).
popef.means, popef.mode
Named lists in the same
shape: the mean of the draws, and the posterior mode used to start
every chain. popef.init records the starting values.
Group level (, one vector per group):
groupefData frame in long format with
n rows: a draw index, a column named
for the grouping factor, then one column per random-effect
coefficient. These are full coefficients, comparable to
lme4::coef(), not mean-zero deviations.
groupef.mode matrix of posterior
modes, rows named by group level.
Variance components:
popef.dispersionn matrix of
between-group variance draws . Constant down a column
for a dNormal component (the value was treated as known);
varying for an dIndependent_Normal_Gamma component.
popef.dispersion.mean is the column mean.
group.dispersionResidual variance : a
fixed scalar or per-group vector, or draws when a dGamma
measurement prior was used. group.dispersion.mean is the
corresponding mean.
Convergence reporting:
m_convergenceInteger sweeps run per chain; 1 on
the exact route.
convergence_infoList with method
("exact_iid" when the draws are exact), lambda_star
and eigenvalues (the contraction rate behind the sweep
count), draw_engine, sim_method_used, and
icm_info for the mode solve. convergence is the same
list.
Inputs echoed back: design, family, n,
call, pfamily_list, dispprior_list,
prior.weights, offset, plus the packed prior and
thin Prior summaries.
When simulate = FALSE, every draw and convergence slot is
NULL and only the mode fields are populated.
rglmerb, rlmerb_point,
rLMMNormal_reg_known_vcov,
rLMMNormal_reg_estimated_vcov,
rLMMindepNormalGamma_reg_known_vcov,
rLMMindepNormalGamma_reg_estimated_vcov,
Prior_Setup_GLMM,
rlmb
Per-group observation-level summaries when
dispersion_ranef is per-group: either a named list of
dGamma() priors (dispersion_mode == "gamma_list") or
a named numeric vector of fixed, known values
(dispersion_mode == "fixed_vector").
For "gamma_list", the layout mirrors the Block~2
tables in summary.lmerb: a prior reference table (with a
common lmer/glmer residual column), an optional posterior
overview from stored draws (including per-group Cand/draw,
envelope candidates per inner sweep), and per-group distribution
percentiles of the draws (1\
99\
the prior table shows whether the simulation tends to pile up against the
truncation bounds. The prior table also reports each group's implied
effective prior sample size (n_prior, back-computed from
shape_ING under the Prior_Setup() default k = 1
calibration) and observed data sample size (n_data); the overview
table reports their sum (n_total).
For "fixed_vector", there is no prior to calibrate and nothing is
sampled (each group's is known and constant throughout),
so only a single table of the fixed per-group values next to the
lmer/glmer residual variance is returned; overview
and percentiles are NULL.
summary_sigma2(object, ...) ## S3 method for class 'lmerb' summary_sigma2( object, type = c("both", "prior", "overview"), digits = max(3L, getOption("digits") - 3L), ... ) ## S3 method for class 'glmerb' summary_sigma2( object, type = c("both", "prior", "overview"), digits = max(3L, getOption("digits") - 3L), ... ) ## S3 method for class 'summary.sigma2.lmerb' print(x, digits = NULL, ...) ## S3 method for class 'summary.sigma2.lmerb' as.data.frame(x, row.names = NULL, ...)summary_sigma2(object, ...) ## S3 method for class 'lmerb' summary_sigma2( object, type = c("both", "prior", "overview"), digits = max(3L, getOption("digits") - 3L), ... ) ## S3 method for class 'glmerb' summary_sigma2( object, type = c("both", "prior", "overview"), digits = max(3L, getOption("digits") - 3L), ... ) ## S3 method for class 'summary.sigma2.lmerb' print(x, digits = NULL, ...) ## S3 method for class 'summary.sigma2.lmerb' as.data.frame(x, row.names = NULL, ...)
object |
An |
type |
|
digits |
Number of significant digits for printing. |
x |
An object of class |
row.names |
Ignored. |
... |
Ignored. |
An object of class "summary.sigma2.lmerb" with components
prior, overview, percentiles (overview and
percentiles are always NULL for "fixed_vector"),
group_name, n_groups, simulated, mode, and
mer_label. When window_diagnostics were attached by
dGamma_list(), the prior table includes
blup_infl, R_lo, R_hi, and asymmetric_window.
Summary method for row-block glmb fits (class "bglmb").
summary.bglmb applies summary.glmb to each block;
print.summary.bglmb follows summary.mlmb /
summary.blmb with per-block sections.
## S3 method for class 'bglmb' summary(object, ...) ## S3 method for class 'summary.bglmb' print(x, digits = max(3, getOption("digits") - 3), ...)## S3 method for class 'bglmb' summary(object, ...) ## S3 method for class 'summary.bglmb' print(x, digits = max(3, getOption("digits") - 3), ...)
object |
An object of class |
x |
An object of class |
digits |
Number of significant digits for printing. |
... |
Passed to |
summary.bglmb returns a named list of "summary.glmb"
objects with class "summary.bglmb".
glmbBlock, lmbBlock,
summary.mlmb, summary.blmb,
summary.glmb
Summary method for row-block lmb fits (class "blmb").
## S3 method for class 'blmb' summary(object, ...) ## S3 method for class 'summary.blmb' print(x, digits = max(3, getOption("digits") - 3), ...)## S3 method for class 'blmb' summary(object, ...) ## S3 method for class 'summary.blmb' print(x, digits = max(3, getOption("digits") - 3), ...)
object |
An object of class |
x |
An object of class |
digits |
Number of significant digits for printing. |
... |
Passed to |
Methods for lmerb fits. summary.lmerb builds
Block~2 (level-2 fixed effect) tables per random-effects component,
following the layout of summary.glmb and the
multi-response structure of summary.mlmb.
## S3 method for class 'lmerb' summary(object, groups = NULL, digits = max(3L, getOption("digits") - 3L), ...) ## S3 method for class 'glmerb' summary(object, groups = NULL, digits = max(3L, getOption("digits") - 3L), ...) ## S3 method for class 'summary.lmerb' print(x, digits = max(3L, getOption("digits") - 3L), ...)## S3 method for class 'lmerb' summary(object, groups = NULL, digits = max(3L, getOption("digits") - 3L), ...) ## S3 method for class 'glmerb' summary(object, groups = NULL, digits = max(3L, getOption("digits") - 3L), ...) ## S3 method for class 'summary.lmerb' print(x, digits = max(3L, getOption("digits") - 3L), ...)
object |
An object of class |
groups |
Optional character vector of grouping levels for which to
include a per-group Block~1 (random effects) detail table. When
|
digits |
Number of significant digits for printing. |
x |
An object of class |
... |
Ignored. |
summary.lmerb returns an object of class
"summary.lmerb", a list with components call,
formula, n, simulated, varcor,
fixef_prior_overview (stacked prior and glmer/lmer
reference across RE components),
fixef_overview (Block~2 hyperparameters with posterior summaries and
Pr(Prior_tail)),
fixef_percentiles_overview (stacked distribution percentiles across
RE components, when simulated),
popef (per-RE-component tables; not printed, available on the
returned object),
ranef_overview, ranef.iters.mean (Block~1 envelope candidates
per inner sweep, averaged over groups; printed separately from the overview
table), any_non_normal, tau2_prior_overview (per-component
prior reference: Block~2 pfamily name (Prior),
, , truncation window,
sqrt(E[tau2]), and lmer/glmer MLE),
sim_method_used ("DEFAULT" for the exact-iid engine,
"TWO_BLOCK_GIBBS" for the two-block Gibbs engine, or NULL
when simulated is FALSE; see lmerb's
sim_method),
tau2_overview and tau2_percentiles_overview (posterior mode,
mean, SD on the variance scale, Mean SD, and tau^2 quantiles when
simulated), tau2_sd_percentiles_overview (2.5\
sqrt(tau^2) draws vs lmer/glmer SD), a Residual row
for observation-level when dispersion_ranef is
supplied (same columns as the rows; not used for
per-group gamma_list priors — see summary_sigma2),
and optionally ranef_groups.
lmerb, glmerb, print.lmerb,
summary_sigma2, summary.glmb,
summary.mlmb
ranef.lmerb objectsSummarize ranef.lmerb objects
## S3 method for class 'ranef.lmerb' summary(object, groups = NULL, digits = 4L, ...)## S3 method for class 'ranef.lmerb' summary(object, groups = NULL, digits = 4L, ...)
object |
An object of class |
groups |
Optional character vector of grouping levels to highlight. |
digits |
Number of digits for printing. |
... |
Ignored. |
A data.frame (invisibly) with long-format random effects.