glmbayes

CRAN status CRAN downloads Monthly downloads GitHub release (latest by date) License: GPL-3 GitHub Workflow Status

glmbayes provides independent and identically distributed (iid) samples for Bayesian Generalized Linear Models (GLMs). Its primary interface, glmb(), serves as a Bayesian analogue to R’s glm() function, supporting Gaussian, Poisson, Binomial, and Gamma families under log-concave likelihoods. Sampling for most models is performed using accept-reject methods based on likelihood subgradients (Nygren and Nygren, 2006). For Gaussian models, the package also includes lmb(), a Bayesian counterpart to R’s lm().

The package includes a rich set of supporting tools for prior specification, model diagnostics, and method functions that mirror those for lm() and glm(). Most functions are extensively documented, and a comprehensive set of vignettes are available to guide users through the package’s capabilities.

This repository is 0.9.6 in development. The current CRAN release is version 0.9.5 (CRAN). The GitHub repository holds the source; R-Universe builds binaries from it. See NEWS.md for changes.

Installation

CRAN (release 0.9.5)

install.packages("glmbayes")

GitHub / R-Universe (install from both CRAN and R-Universe repositories if you want R-Universe binaries or faster mirrors):

install.packages("glmbayes",
                 repos = c("https://cloud.r-project.org",
                           "https://knygren.r-universe.dev"))

Prebuilt binaries from CRAN (0.9.5) and R-Universe are built without OpenCL GPU support. For the CRAN release, OpenCL requires installing from source on a system with OpenCL development files available. To set up GPU acceleration, follow

Chapter 16 — Large models: GPU acceleration using OpenCL https://knygren.r-universe.dev/articles/glmbayes/Chapter-16.html

Minimal Working Example

library(glmbayes)

# Dobson (1990), p. 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 glm
glm.D93 <- glm(counts ~ outcome + treatment,
               family = poisson())

## Bayesian glmb
# Step 1: Set up prior
ps <- Prior_Setup(counts ~ outcome + treatment, family = poisson())
mu <- ps$mu
V  <- ps$Sigma

# Step 2: Fit using glmb
glmb.D93 <- glmb(counts ~ outcome + treatment,
                 family = poisson(),
                 pfamily = dNormal(mu = mu, Sigma = V))

summary(glmb.D93)

As with glm(), models are defined by a formula for the linear predictor and a family() describing the likelihood and link. In addition, glmb() requires a pfamily object specifying the prior.

Priors on regression coefficients

The primary table below covers priors on the regression coefficients β. The standard prior for all families is dNormal. The conjugate priors dBeta and dGamma(Inv_Dispersion = FALSE) provide closed-form IID posterior draws for intercept-only models with an identity link.

Likelihood family Link functions Compatible pfamilies (coefficient priors)
Gaussian identity dNormal, dNormal_Gamma, dIndependent_Normal_Gamma
Poisson / Quasi-Poisson log dNormal
Poisson identity (intercept-only) dGamma(Inv_Dispersion = FALSE) — conjugate Gamma–Poisson rate prior
Binomial / Quasi-Binomial logit, probit, cloglog dNormal
Binomial identity (intercept-only) dBeta — conjugate Beta–Binomial probability prior
Gamma log dNormal
Gamma identity (intercept-only) dGamma(Inv_Dispersion = FALSE) — conjugate Gamma–Gamma rate prior

dNormal_Gamma and dIndependent_Normal_Gamma also model precision jointly with the coefficients; see the precision/dispersion table below.

Priors on precision / dispersion

dGamma(Inv_Dispersion = TRUE) (the default when Inv_Dispersion is omitted) places a Gamma prior on the inverse dispersion 1/φ with the regression coefficients β held fixed. This is the precision prior used in Gibbs sampling steps for dispersion estimation.

Likelihood family Link Compatible pfamilies (precision prior)
Gaussian identity dGamma — prior on 1/σ² (precision)
Gamma log dGamma — prior on 1/φ (shape / dispersion)

dNormal_Gamma and dIndependent_Normal_Gamma model β and precision jointly in a single conjugate step, avoiding the need for a separate Gibbs precision update.

Prior_Setup

For a default, data‑aligned prior using the same formula and family as glm(), call Prior_Setup(formula, family, data = ..., ...). The returned list includes default settings for the following:

Optional arguments adjust prior weight, centering, and related settings (see the function help and vignette Chapter 04).

Typical Prior_Setup wiring

Assuming ps <- Prior_Setup(...):

The default priors have limiting behaviors that produce estimates resembling classical estimates as priors get weak (see documentation and vignettes for details).

All supported models have log‑concave likelihoods, enabling efficient iid sampling via enveloping functions and subgradient‑based accept–reject algorithms, especially for models lacking standard iid samplers.

Examples and Demos

Use example() and demo() to explore built-in examples and demos for supported families and links:

## Bayesian linear regression
example("lmb")

## Bayesian generalized linear models
example("glmb")

## Beta-Binomial conjugacy: dBeta() prior; Bechdel test (requires bayesrules)
## See also: vignette("Chapter-02-S03", package = "glmbayes")
demo("Ex_12_BetaBinomial")

## Gamma-Poisson conjugacy: dGamma(Inv_Dispersion=FALSE); bike counts + heart
## transplant mortality (requires bayesrules; Appendix A requires LearnBayes)
## See also: vignette("Chapter-02-S04", package = "glmbayes")
demo("Ex_13_GammaPoisson")

## Predictions for fitted glmb objects (newdata, type, etc.)
example("predict.glmb")

## Deviance residuals and simulate() for posterior predictive checks (menarche)
example("residuals.glmb")

## Two-block Gibbs sampler compared with iid sampling (linear model)
example("rlmb")

## Default prior specification using Prior_Setup
example("Prior_Setup")

## Matrix-input GLM example with an informative prior
example("rglmb")

## Two-step Boston example: estimates and summarizes models with unknown
## dispersion using dGamma priors via rGamma_reg, rglmb, rlmb, glmb, and lmb
example("summary.rGamma_reg")

## High-dimensional Gaussian model (14 predictors) with GPU acceleration (requires OpenCL)
example("Boston_centered")

## High-dimensional binomial model (14 predictors) with GPU acceleration (requires OpenCL)
example("Cleveland")

## Hierarchical linear model (Rubin/Gelman 8-schools) via rlmb
demo("Ex_07_Schools")

## Hierarchical generalized linear model (Poisson BikeSharing) via rglmb
demo("Ex_09_BikeSharingPoisson")

## Detailed simulation pipeline for rNormalGLM models (JASA 2006; Vignette Chapter A05)
example("rNormalGLM_std")

## Detailed simulation pipeline for rIndepNormalGammaReg models (Vignette Chapter A07)
example("rIndepNormalGammaReg_std")

Methodology

For generalized linear models where well known sampling methods are unavailable, sampling follows the framework from Nygren and Nygren (2006), using likelihood subgradients to construct enveloping functions for the posterior distribution. When the posterior is approximately normal, the expected number of draws per acceptance is bounded as per that paper and as discussed in our vignettes. Dispersion can be sampled via rGamma_reg() (standalone) or jointly with coefficients via rNormalGamma_reg() and rindepNormalGamma_reg().

GPU Acceleration Using OpenCL

The implemented algorithms tend to have acceptable performance on CPUs up to around 10-14 dimensions. For larger models, the envelope construction is embarrassingly parallel. To accelerate envelope construction in such cases, the package provides optional GPU acceleration using OpenCL. This requires that users have GPU enabled machines and an OpenCL installation. These features are discussed in more detail in two of our vignettes.

Vignettes

The glmbayes package includes a comprehensive set of vignettes organized into five major parts. These vignettes guide users from introductory material through applied modeling, advanced topics, and the underlying simulation methods that support the package.

Part 1: An Introduction

Overview of the package, its design philosophy, single-parameter conjugate models, and the basic workflow for fitting Bayesian linear and generalized linear models.

Part 2: Bayesian regression models

These chapters focus on Bayesian linear regression (Gaussian family). Topics include lmb() fitting, Prior_Setup(), posterior predictive checks (bayesplot), deviance residuals and model summaries, bayestestR-style summaries, and the bridge to Bayesian GLMs in Part 3.

Part 3: Generalized Linear Models

This part presents Bayesian GLMs across the major likelihood families, including binomial, quasi-binomial, Poisson, quasi-Poisson, and Gamma models. It covers model specification, link functions, log-concavity, diagnostics, interpretation of posterior results, and tooling (bayesplot, bayestestR) for visualization and summaries.

Part 4: Advanced Topics

These chapters explore more complex modeling scenarios and computational strategies, such as informative priors, two-block Gibbs sampling, linear and generalized linear mixed-effects models, models with unknown dispersion parameters, and large-scale model fitting using GPU acceleration using OpenCL.

Part 5: Simulation Methods and Technical Implementation

This part documents the mathematical and algorithmic foundations of the package. Topics include estimation procedures, likelihood subgradient densities, envelope construction, accept-reject sampling, and technical reports on sampler design including implementation aspects for GPU acceleration using OpenCL.

Together, these vignettes form a comprehensive reference that supports users at all levels, from first-time Bayesian GLM users to researchers interested in the mathematical and computational details behind the samplers.

Feature Highlights

Limitations

Future Plans

Further performance and algorithm work: