simfuncs: Simulation Functions for Bayesian Generalized Linear Models
In glmbayes: Bayesian Generalized Linear Models (IID Samples)

simfuncs

R Documentation

Simulation Functions for Bayesian Generalized Linear Models

Description

Simulation functions provide a unified interface for generating posterior samples from Bayesian GLMs. These functions are typically used within model fitting routines such as rglmb and rlmb, and are also suitable for use in Block Gibbs sampling and other simulation-based inference techniques.

Usage

simfunction(object, ...)

rNormal_reg(n, y, x, prior_list, offset = NULL, weights = 1,
            family = gaussian(), Gridtype = 2, n_envopt = NULL,
            use_parallel = TRUE, use_opencl = FALSE, verbose = FALSE,progbar=FALSE)

rNormalGamma_reg(n, y, x, prior_list, offset = NULL, weights = 1, family = gaussian(),
                  Gridtype = 2,n_envopt = NULL, 
                  use_parallel = TRUE, use_opencl = FALSE, verbose = FALSE,progbar=FALSE)

rindepNormalGamma_reg(n, y, x, prior_list, offset = NULL, weights = 1,
                             family = gaussian(), Gridtype = 2,n_envopt = NULL,
                              use_parallel = TRUE, use_opencl = FALSE, verbose = FALSE, 
                             progbar = TRUE)

rGamma_reg(n, y, x, prior_list, offset = NULL, weights = 1, family = gaussian(),
           Gridtype = 2,n_envopt = NULL,
            use_parallel = TRUE, use_opencl = FALSE, verbose = FALSE,progbar=FALSE)

## S3 method for class 'rGamma_reg'
print(x, digits = max(3, getOption("digits") - 3), ...)

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

rGamma_Conjugate_reg(n, y, x, prior_list, offset = NULL, weights = 1, family = gaussian(),
           Gridtype = 2,n_envopt = NULL,
            use_parallel = TRUE, use_opencl = FALSE, verbose = FALSE,progbar=FALSE)

rBeta_reg(n, y, x, prior_list, offset = NULL, weights = 1,
                family = gaussian(), Gridtype = 2, n_envopt = NULL,
                use_parallel = TRUE, use_opencl = FALSE,
                verbose = FALSE, progbar = FALSE)

Arguments

`object`	A fitted model object containing a `pfamily` component. The generic function `simfunction()` accesses the simulation metadata stored within such objects.
`n`	Number of draws to generate. If `length(n) > 1`, the length is taken to be the number required.
`y`	A vector of observations of length `m`.
`x`	for the simulation functions a design matrix of dimension `m * p` and for the print functions the object to be printed.
`prior_list`	A list with prior parameters (e.g., shape, rate, beta) used in the simulation.
`offset`	Optional numeric vector of length `m` specifying known components of the linear predictor.
`weights`	Optional numeric vector of prior weights.
`family`	A description of the error distribution and link function (see `family`).
`Gridtype`	Optional integer specifying the method used to construct the envelope function.
`n_envopt`	Effective sample size passed to EnvelopeOpt for grid construction. Defaults to match `n`. Larger values encourage tighter envelopes.
`use_parallel`	Logical. Whether to use parallel processing.
`use_opencl`	Logical. Whether to use OpenCL acceleration.
`verbose`	Logical. Whether to print progress messages.
`progbar`	Logical. Whether to display a progress base during simulation.
`digits`	Number of significant digits to use for printed output.
`...`	Additional arguments passed to or from other methods.

Details

The low-level simulation functions rNormal_reg(), rNormalGamma_reg(), rindepNormalGamma_reg(), and rGamma_reg() generate iid samples from posterior distributions for specific model components. These model functions are used internally by the functions rglmb() and rlmb() to generate samples.

The simfunction() generic extracts metadata from simulation objects, including the function name, call, and arguments used. This is useful for introspection, reproducibility, and diagnostics.

The lower-level simulation functions generate iid samples from posterior distributions for specific model components. These functions are used internally by pfamily constructors and model fitting routines.

Simulation Functions

rNormal_reg(): Produces iid draws for regression coefficients in models with multivariate normal priors and log-concave likelihood functions. For Gaussian likelihoods, these are conjugate priors and standard simulation procedures for multivariate normal distributions are utilized \insertCiteLindleySmith1972,DiaconisYlvisaker1979glmbayes. For all other families/link functions, the likelihood subgradient approach of \insertCiteNygren2006glmbayes is used to generate iid samples.
rNormalGamma_reg(): Produces iid draws for regression coefficients and the dispersion parameter in models with Normal-Gamma priors and Gaussian likelihoods, where this is a conjugate prior distribution. Standard simulation procedures for gamma and multivariate normal distributions are utilized \insertCiteRaiffa1961,LindleySmith1972glmbayes.
rindepNormalGamma_reg(): Produces iid draws for regression coefficients and the dispersion parameter in models with independent Normal and truncated Gamma priors. This is a non-conjugate specification but can still be sampled using accept-reject procedures based on an enveloping approach (see vignette \insertCiteglmbayesIndNormGammaVignetteglmbayes).
rGamma_reg(): Simulates dispersion parameters for Gaussian and Gamma families using either standard gamma sampling or accept-reject methods based on likelihood subgradients \insertCiteChen1979,glmbayesGammaVignetteglmbayes.

Value

simfunction()

An object of class "simfunction" containing:

name: Character string with the name of the simulation function.
call: The matched call used to generate the simulation.
args: A named list of arguments passed to the simulation function.

rNormal_reg()

A list object with classes "rglmb", "glmb", "glm", and "lm". Elements include:

coefficients: Matrix (n * p) of simulated regression coefficients, with column names from x.
coef.mode: Posterior mode of the coefficients. Gaussian: from lm.fit; non-Gaussian: BFGS mode shifted by prior mean.
dispersion: Scalar dispersion used. Poisson/Binomial: 1; otherwise the supplied value. Quasi families: mean residual-based dispersion computed in the wrapper.
Prior: List with mean (prior mean vector) and Precision (prior precision matrix P).
prior.weights: Vector of prior weights used in the simulation (unscaled).
offset: Offset vector passed to the C++ sampler.
offset2: Offset used internally by the wrapper (copy of input or a zero vector).
y: Response vector.
x: Design matrix.
fit: Fitted/diagnostic object. Gaussian: result of lm.fit (class "lm"). Non-Gaussian: result of glmb.wfit(...).
iters: Vector of iteration counts per sample. Gaussian: vector of ones; non-Gaussian: counts from the sampler.
Envelope: Envelope list used for accept-reject sampling (non-Gaussian); NULL for Gaussian.
family: Family object describing distribution and link.
famfunc: Processed family functions used internally (e.g., f2, f3).
call: Matched call to rNormal_reg().
formula: Formula reconstructed from y and x.
model: Model frame corresponding to formula.
data: Data frame combining y and x.

rNormalGamma_reg()

A list with class "rglmb" containing:

coefficients: Matrix (n * p) of simulated regression coefficients; row i equals Btilde + IR %*% rnorm(p) * sqrt(dispersion[i]). Column names are set to colnames(x).
coef.mode: Posterior mean/mode vector Btilde from rNormal_reg.wfit().
dispersion: Numeric vector of length n with draws from the inverse-gamma posterior 1/rgamma(shape = shape + nobs/2, rate = rate + 0.5*S).
Prior: List with mean (as numeric vector mu) and Precision (matrix P).
offset: Offset vector as supplied.
prior.weights: Vector of prior weights wt.
y: Response vector.
x: Design matrix.
fit: Result from rNormal_reg.wfit(), including fields such as Btilde, IR, S, and k.
famfunc: Processed family functions for Gaussian models (from glmbfamfunc(gaussian())).
iters: Numeric vector (length n) of ones indicating per-draw iteration counts.
Envelope: NULL; no envelope is constructed in this conjugate setup.
call: Matched call to rNormalGamma_reg().

rindepNormalGamma_reg()

A list with class "rglmb" containing:

coefficients: Matrix (n * p) of simulated regression coefficients, back-transformed to the original scale; column names set to colnames(x).
coef.mode: Vector with the conditional posterior mode used for envelope anchoring (from the Gaussian fit).
dispersion: Numeric vector of length n with simulated dispersion draws.
Prior: List with prior components: mean (prior mean mu), Sigma (prior covariance), shape and rate (Gamma prior for dispersion), Precision (solve(Sigma)).
family: The gaussian() family object.
prior.weights: Vector of prior weights used in the simulation.
y: Response vector.
x: Design matrix.
call: Matched call to rindepNormalGamma_reg().
famfunc: Processed family functions for Gaussian models (from glmbfamfunc).
iters: Vector with per-draw iteration counts returned by the joint sampler.
Envelope: NULL; envelope diagnostics are not returned by this function.
loglike: NULL; placeholder for log-likelihood values.
weight_out: Numeric vector of per-draw weights returned by the C++ routine.
sim_bounds: List with low and upp, the dispersion bounds used by the shared envelope.
offset2: Offset vector used internally (copy of input or a zero vector).

rGamma_reg()

An object of class "rGamma_reg" containing:

coefficients: A 1 * p matrix of assumed regression coefficients.
coef.mode: Currently NULL; reserved for future use.
dispersion: A vector of simulated dispersion values.
Prior: A list with prior parameters: shape and rate.
prior.weights: Vector of prior weights used in the simulation.
y: The response vector.

Author(s)

The simulation framework was developed by Kjell Nygren as part of the glmbayes package. It builds on the likelihood subgradient approach described in \insertCiteNygren2006glmbayes, and extends classical Bayesian GLM sampling techniques.

References

\insertAllCited

Examples

############################### Start of rNormal_reg examples ####################
## During CRAN checks, run examples sequentially.
use_parallel <- identical(Sys.getenv("NOT_CRAN"), "true")

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))

## Poisson Prior and rNormal_reg call (using Prior_Setup for x, y, and prior values)
ps <- Prior_Setup(counts ~ outcome + treatment, family = poisson(), data = d.AD)

out_pois <- rNormal_reg(
  n = 1000,
  y = ps$y,
  x = ps$x,
  prior_list = list(mu = ps$mu, Sigma = ps$Sigma),
  family = poisson(link = "log"),
  weights = rep(1, nrow(ps$x)),
  use_parallel = use_parallel
)
summary(out_pois)


## Menarche Binomial Data Example
data(menarche, package = "MASS")
menarche$Age2 <- menarche$Age - 13

## Logit Prior and rNormal_reg call (use proportion + trial weights)
ps1 <- Prior_Setup(
  Menarche / Total ~ Age2,
  family = binomial(logit),
  data = menarche,
  weights = menarche$Total
)

out_logit <- rNormal_reg(
  n = 1000,
  y = ps1$y,
  x = ps1$x,
  prior_list = list(mu = ps1$mu, Sigma = ps1$Sigma),
  family = binomial(logit),
  weights = menarche$Total,
  use_parallel = use_parallel
)
summary(out_logit)

## Probit Prior and rNormal_reg call
ps2 <- Prior_Setup(
  Menarche / Total ~ Age2,
  family = binomial(probit),
  data = menarche,
  weights = menarche$Total
)

out_probit <- rNormal_reg(
  n = 1000,
  y = ps2$y,
  x = ps2$x,
  prior_list = list(mu = ps2$mu, Sigma = ps2$Sigma),
  family = binomial(probit),
  weights = menarche$Total,
  use_parallel = use_parallel
)
summary(out_probit)

## clog-log Prior and rNormal_reg call
ps3 <- Prior_Setup(
  Menarche / Total ~ Age2,
  family = binomial(cloglog),
  data = menarche,
  weights = menarche$Total
)

out_cloglog <- rNormal_reg(
  n = 1000,
  y = ps3$y,
  x = ps3$x,
  prior_list = list(mu = ps3$mu, Sigma = ps3$Sigma),
  family = binomial(cloglog),
  weights = menarche$Total,
  use_parallel = use_parallel
)
summary(out_cloglog)


## Gamma regression
data(carinsca)
carinsca$Merit <- ordered(carinsca$Merit)
carinsca$Class <- factor(carinsca$Class)
oldopt <- options(contrasts = c("contr.treatment", "contr.treatment"))

psg <- Prior_Setup(
  Cost / Claims ~ Merit + Class,
  family = Gamma(link = "log"),
  data = carinsca,
  weights = carinsca$Claims
)

out_gamma <- rNormal_reg(
  n = 1000,
  y = psg$y,
  x = psg$x,
  prior_list = list(mu = psg$mu, Sigma = psg$Sigma, dispersion = psg$dispersion),
  family = Gamma(link = "log"),
  weights = carinsca$Claims,
  use_parallel = use_parallel
)
summary(out_gamma)
options(oldopt)
############################### Start of rNormalGamma_reg examples ####################
## During CRAN checks, run examples sequentially.
use_parallel <- identical(Sys.getenv("NOT_CRAN"), "true")

## 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)
mu <- ps$mu
shape <- ps$shape
rate <- ps$rate

y <- ps$y
x <- as.matrix(ps$x)
prior_list <- list(mu = mu, Sigma = ps$Sigma_0, shape = shape, rate = rate)
ngamma.D9 <- rNormalGamma_reg(n = 1000, y = y, x = x,
  prior_list = prior_list, use_parallel = use_parallel)

summary(ngamma.D9)
############################### Start of rindepNormalGamma_reg examples ####################
## During CRAN checks, run examples sequentially.
use_parallel <- identical(Sys.getenv("NOT_CRAN"), "true")

## 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)
p_setup <- Prior_Setup(weight ~ group, family = gaussian())

mu <- p_setup$mu
Sigma_prior <- p_setup$Sigma
dispersion <- p_setup$dispersion
shape <- p_setup$shape
rate <- p_setup$rate
y <- p_setup$y
x <- p_setup$x

prior_list <- list(
  mu = mu,
  Sigma = Sigma_prior,
  dispersion = dispersion,
  shape = shape,
  rate = rate,
  Precision = solve(Sigma_prior),
  max_disp_perc = 0.99
)


set.seed(360)

sim2 <- rindepNormalGamma_reg(n = 1000, y, x, prior_list = prior_list,
  use_parallel = use_parallel)
summary(sim2)
 
 
 
 
 
############################### Start of rGamma_reg examples ####################
## During CRAN checks, run examples sequentially.
use_parallel <- identical(Sys.getenv("NOT_CRAN"), "true")

## 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)

## Set up prior hyperparameters (shape/rate) and model matrix via Prior_Setup
ps <- Prior_Setup(weight ~ group, family = gaussian())
y <- ps$y
x <- as.matrix(ps$x)

## rGamma_reg uses a dGamma-style prior on dispersion with fixed beta.
## Use coefficients from Prior_Setup (full-model GLM MLE by default).
prior_list <- list(beta = ps$coefficients, shape = ps$shape, rate = ps$rate)

out <- rGamma_reg(n = 1000, y = y, x = x, prior_list = prior_list,
  family = gaussian(), use_parallel = use_parallel)
summary(out)

glmbayes documentation built on Aug. 5, 2026, 1:07 a.m.