Prior_Setup: Setup Prior Objects
In glmbayes: Bayesian Generalized Linear Models (IID Samples)

Prior_Setup

R Documentation

Setup Prior Objects

Description

Helper function to facilitate the Setup of Prior Distributions for glm models.

Usage

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

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

Arguments

`formula`	an object of class `"formula"` (or one that can be coerced to that class): a symbolic description of the model to be fitted. The details of model specification are given under ‘Details’.
`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 `as.data.frame` to a data frame) containing the variables in the model. If not found in `data`, the variables are taken from `environment(formula)`, typically the environment from which `glm` is called.
`weights`	an optional vector of ‘prior weights’ to be used in the fitting process. Should be `NULL` or a numeric vector.
`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 `NAs` are treated. The default is first, any `na.action` attribute of data, second a `na.action` setting of options, and third `na.fail` if that is unset. The `factory-fresh` default is `na.omit`. Another possible value is `NULL`.
`offset`	this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be `NULL` or a numeric vector of length equal to the number of cases. One or more `offset` terms can be included in the formula instead or as well, and if more than one is specified their sum is used. See `model.offset`.
`contrasts`	an optional list. See the `contrasts.arg` of `model.matrix.default`.
`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 `n_prior` is provided and `pwt` is still a scalar and `sd` was not supplied, `pwt` is set to `n_prior / (n_prior + n_effective)`. If `length(pwt) > 1` (including from `sd`) or `sd` was supplied, `n_prior` does not overwrite `pwt`; it is used only as a scalar for Gamma / S_marg steps. If `sd` is provided, `pwt` is computed from the prior standard deviations. If none of these are supplied, `pwt` defaults to `pwt_default_low` for models with fewer than 14 coefficients, and `pwt_default_high` otherwise.
`pwt_default_low`	Default prior weight used when `pwt` is not supplied and the model dimension is below 14. Defaults to 0.01.
`pwt_default_high`	Default prior weight used when `pwt` is not supplied and the model dimension is 14 or greater. Defaults to 0.05.
`n_prior`	Optional scalar effective prior sample size (on the `n_effective` scale). If provided with scalar `pwt` and without `sd`, `pwt` is recomputed from `n_prior`. With vector `pwt` or with `sd`, `pwt` is left unchanged and `n_prior` is used for the Gamma prior on precision and related Gaussian calibration only. If missing and `pwt` is scalar, `n_prior = (pwt/(1-pwt))*n_effective`.
`sd`	Optional vector argument with the prior standard deviations for the coefficients
`dispersion`	Optional scalar dispersion override (default `NULL`). For now, this is documented as an optional argument used to scale the `Sigma` (variance-covariance) matrix; see Details for additional context.
`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 `1`), non-negative (`k \geq 0`), with `k + p \geq 2` where `p` is the number of coefficients (columns of the model matrix). `k` controls the tail behavior and effective degrees of freedom of the variance prior. It does not change the posterior mean of `\sigma^2` or the covariance of `\beta`, but larger `k` makes the prior and posterior for `\sigma^2` more concentrated and less heavy-tailed. Not yet used in calibration; passed through to `compute_gaussian_prior` for future use.
`...`	For `glm`: arguments to be used to form the default `control` argument if it is not supplied directly. For `weights`: further arguments passed to or from other methods.
`x`	An object of class `"PriorSetup"`

Details

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 V_0 = \mathrm{vcov}(\hat\beta) be the covariance matrix of the full-model GLM coefficients. For non-Gaussian families, the prior covariance is:

\Sigma = \begin{cases} \dfrac{1 - \mathrm{pwt}}{\mathrm{pwt}} V_0, & \text{scalar pwt},\\[4pt] V_0 \circ \left[\sqrt{\dfrac{1 - \mathrm{pwt}_i}{\mathrm{pwt}_i}} \sqrt{\dfrac{1 - \mathrm{pwt}_j}{\mathrm{pwt}_j}}\right], & \text{vector pwt}, \end{cases}

where \circ 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 n_{\mathrm{prior}} already satisfies n_{\mathrm{prior}} / (n_{\mathrm{prior}} + n_{\mathrm{effective}}) = \mathrm{pwt} (otherwise conj_poisson remains NULL with a warning). Writing the weighted mean \bar{y}_w = \sum_i w_i y_i / \sum_i w_i, the output list component conj_poisson stores \texttt{shape} = n_{\mathrm{prior}} \bar{y}_w and \texttt{rate} = n_{\mathrm{prior}}, so the prior mean for the rate matches \bar{y}_w. Omitting optional mu resets the surrogate Normal summaries mu and the sole diagonal element of Sigma to these Gamma moments (Sigma_{11} = \bar{y}_w / n_{\mathrm{prior}}).

For Gaussian families, Prior_Setup() also constructs the dispersion-free covariance

\Sigma_0 = \Sigma / \texttt{dispersion},

which under scalar pwt and the default calibration reduces to

\Sigma_0 = \frac{1 - \mathrm{pwt}}{\mathrm{pwt}} (X^\top W X)^{-1}.

Gaussian Normal–Gamma calibration and S_{\mathrm{marg}}

For family = gaussian(), the function performs the Normal–Gamma calibration described in \insertCiteglmbayesChapterA12glmbayes. Let:

p = \texttt{ncol}(x),
n_{\mathrm{effective}} = \sum_i w_i,
\hat\beta the weighted least-squares estimator,
\Sigma_0 the dispersion-free prior covariance.

The marginal quadratic term is

S_{\mathrm{marg}} = \mathrm{RSS}_w + (\hat\beta - \mu)^\top \left(\Sigma_0 + (X^\top W X)^{-1}\right)^{-1} (\hat\beta - \mu),

where \mathrm{RSS}_w is the weighted residual sum of squares at \hat\beta. Under the default scalar-pwt Zellner mapping \Sigma_0 = \frac{1 - \mathrm{pwt}}{\mathrm{pwt}}(X^\top W X)^{-1}, this simplifies to

S_{\mathrm{marg}} = \mathrm{RSS}_w + \mathrm{pwt}\, (\hat\beta - \mu)^\top (X^\top W X)(\hat\beta - \mu),

which makes the limiting behavior as \mathrm{pwt} \to 0 transparent.

The calibrated dispersion is

\texttt{dispersion} = \frac{S_{\mathrm{marg}}}{n_{\mathrm{effective}} - p},

and the Normal–Gamma hyperparameters are

\text{shape} = \frac{n_{\mathrm{prior}} + k}{2},\qquad \text{rate} = \frac{1}{2} S_{\mathrm{marg}} \frac{n_{\mathrm{prior}} + k + p - 2}{n_{\mathrm{effective}} - p}.

The independent Normal–Gamma shape is

\text{shape}_{ING} = \text{shape} + \frac{p}{2}.

Posterior summaries for the conjugate Normal–Gamma prior

Under the conjugate Normal–Gamma prior (used by dNormal_Gamma()), the posterior has:

Posterior mean

E[\beta \mid y] = (1 - \mathrm{pwt})\,\hat\beta + \mathrm{pwt}\,\mu.
Posterior expectation of \sigma^2

E[\sigma^2 \mid y] = \frac{S_{\mathrm{marg}}}{n_{\mathrm{effective}} - p}.
Posterior covariance

\mathrm{Cov}(\beta \mid y) = E[\sigma^2 \mid y]\, \left(\Sigma_0^{-1} + X^\top W X\right)^{-1}.

Weak-prior limits (Theorems 2 and 3)

As n_{\mathrm{prior}} \to 0^+ (equivalently \mathrm{pwt} \to 0), S_{\mathrm{marg}} \to \mathrm{RSS}_w, and the conjugate Normal–Gamma posterior converges to the classical weighted least-squares limit:

E[\beta \mid y] \to \hat\beta,\qquad E[\sigma^2 \mid y] \to \frac{\mathrm{RSS}_w}{n_{\mathrm{effective}} - p},\qquad \mathrm{Cov}(\beta \mid y) \to \frac{\mathrm{RSS}_w}{n_{\mathrm{effective}} - p} (X^\top W X)^{-1}.

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 \insertCiteglmbayesChapterA12glmbayes shows that the ING posterior has the same weak-prior limit as the conjugate Normal–Gamma posterior:

E[\beta \mid y] \to \hat\beta,\qquad \mathrm{Cov}(\beta \mid y) \to \frac{\mathrm{RSS}_w}{n_{\mathrm{effective}} - p} (X^\top W X)^{-1}.

Value

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 `family = gaussian()` only: dispersion-independent prior covariance on the precision-weighted coefficient scale (the `Sigma_0` passed to `compute_gaussian_prior`). Under scalar `pwt`, `\Sigma_0^{-1} = \frac{\mathrm{pwt}}{1-\mathrm{pwt}} X^\top W X`.
`dispersion`	Calibrated dispersion (Gaussian models only), equal to `S_{\mathrm{marg}}/(n_{\mathrm{effective}} - p)` under the default calibration.
`shape`	Derived prior Gamma shape parameter for the Normal–Gamma prior on precision (Gaussian only), `(n_{\mathrm{prior}} + k)/2`.
`shape_ING`	For `gaussian()` only when `shape` is available: dedicated shape parameter for `dIndependent_Normal_Gamma()`, `\texttt{shape} + p/2`.
`rate`	Derived prior Gamma rate parameter (Gaussian only), using the calibrated `S_{\mathrm{marg}}`.
`rate_gamma`	For `gaussian()` only, when Gaussian calibration runs: prior Gamma rate for `dGamma()` / fixed-`\beta` use, based on `\mathrm{RSS}_w(\beta_\star)` at the Zellner blend.
`coefficients`	Named numeric vector of returned coefficient values. For `gaussian()` with scalar or vector `pwt`, this is the closed-form posterior-mean blend `(1-\mathrm{pwt})\hat\beta + \mathrm{pwt}\mu` when inputs are valid; otherwise it falls back to the full-model GLM coefficients.
`model`	The model frame used to construct the design matrix (if `model = TRUE`).
`x`	The model matrix used (if `x = TRUE`).
`y`	The response vector used (if `y = TRUE`).
`call`	The matched call to `Prior_Setup()`.
`PriorSettings`	A list containing prior configuration details, including `pwt`, `n_prior`, `n_effective`, `n_likelihood`, `intercept_source`, and `effects_source`.
`conj_poisson`	`NULL` unless `family = poisson(link = "identity")` with a single-column design (intercept-only), scalar `pwt`, finite positive scalar `n_prior`, and negligible `offset`. Then a named list with conjugate `Gamma(shape, rate)` hyperparameters on the Poisson rate, with `\texttt{shape} = n_{\mathrm{prior}}\,\bar{y}_w`, `\texttt{rate} = n_{\mathrm{prior}}`, and `beta` centered at the weighted sample mean `\bar{y}_w`; pass to `dGamma` with `Inv_Dispersion = FALSE`. See Details.

References

\insertAllCited

Examples

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

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

## Set up Prior for Poisson Model
ps <- Prior_Setup(counts ~ outcome + treatment, family = poisson())
ps

## Normal prior for glmb
glmb.D93 <- glmb(
  counts ~ outcome + treatment,
  family = poisson(),
  pfamily = dNormal(mu = ps$mu, Sigma = ps$Sigma),
  use_parallel = use_parallel
)

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

## Conjugate Normal Prior (fixed dispersion)
lmb.D9 <- lmb(
  weight ~ group,
  pfamily = dNormal(mu = ps2$mu, ps2$Sigma, dispersion = ps2$dispersion),
  use_parallel = use_parallel
)

## Conjugate Normal_Gamma Prior
lmb.D9_v2 <- lmb(
  weight ~ group,
  pfamily = dNormal_Gamma(
    ps2$mu,
    Sigma_0 = ps2$Sigma_0,
    shape = ps2$shape,
    rate = ps2$rate
  ),
  use_parallel = use_parallel
)

## Independent_Normal_Gamma_Prior
lmb.D9_v3 <- lmb(
  weight ~ group,
  dIndependent_Normal_Gamma(
    ps2$mu,
    ps2$Sigma,
    shape = ps2$shape_ING,
    rate = ps2$rate
  ),
  use_parallel = use_parallel
)

## -------------------------------------------------------------------------
## Matrix-input bridge: use Prior_Setup outputs with rglmb() and rlmb()
## -------------------------------------------------------------------------
y <- ps2$y
x <- as.matrix(ps2$x)
wt <- rep(1, length(y))

rglmb.D9 <- rglmb(
  n = 1000,
  y = y,
  x = x,
  pfamily = dIndependent_Normal_Gamma(
    ps2$mu,
    ps2$Sigma,
    shape = ps2$shape_ING,
    rate = ps2$rate
  ),
  weights = wt,
  family = gaussian(),
  use_parallel = use_parallel
)

rlmb.D9 <- rlmb(
  n = 1000,
  y = y,
  x = x,
  pfamily = dIndependent_Normal_Gamma(
    ps2$mu,
    ps2$Sigma,
    shape = ps2$shape_ING,
    rate = ps2$rate
  ),
  weights = wt,
  use_parallel = use_parallel
)

## -------------------------------------------------------------------------
## Prior-list templates for lower-level samplers
## -------------------------------------------------------------------------
prior_list_rNormalGamma <- list(
  mu = ps2$mu,
  Sigma = ps2$Sigma_0,
  shape = ps2$shape,
  rate = ps2$rate
)

prior_list_rindepNormalGamma <- list(
  mu = ps2$mu,
  Sigma = ps2$Sigma,
  dispersion = ps2$dispersion,
  shape = ps2$shape_ING,
  rate = ps2$rate,
  Precision = solve(ps2$Sigma),
  max_disp_perc = 0.99
)

rate_dg <- if (!is.null(ps2$rate_gamma)) ps2$rate_gamma else ps2$rate
prior_list_rGamma <- list(
  beta = ps2$coefficients,
  shape = ps2$shape,
  rate = rate_dg
)

## Note: for a full dGamma run across rGamma_reg/rglmb/rlmb/glmb/lmb, see:
## example("summary.rGamma_reg")

## -------------------------------------------------------------------------
## dGamma prior illustration: Prior_Setup(shape, rate_gamma or rate) + fixed beta
## -------------------------------------------------------------------------
## For gaussian models, Prior_Setup() provides Gamma hyperparameters for the
## precision: tau = 1/dispersion ~ Gamma(shape, rate). The lmb() / rGamma_reg()
## constructors consume these via dGamma(pfamily) and rGamma_reg(prior_list),
## respectively.

## dGamma via pfamily for lmb()
lmb.D9_dGamma <- lmb(
  weight ~ group,
  pfamily = dGamma(shape = ps2$shape, rate = rate_dg, beta = ps2$coefficients),
  use_parallel = use_parallel
)

## dGamma via prior_list for rGamma_reg()
out.rGamma_reg <- rGamma_reg(
  n = 1000,
  y = y,
  x = x,
  prior_list = prior_list_rGamma,
  offset = rep(0, length(y)),
  weights = wt,
  family = gaussian(),
  use_parallel = use_parallel
)

## -------------------------------------------------------------------------
## Poisson(link = "identity"), intercept-only: `conj_poisson` + dGamma(Inv_Dispersion=FALSE)
## -------------------------------------------------------------------------
y_p <- c(rep(1L, 3L), rep(0L, 6L))
df_p <- data.frame(y = y_p)
ps_p <- Prior_Setup(
  y ~ 1,
  family = poisson(link = "identity"),
  data = df_p,
  pwt = 0.4
)
if (!is.null(ps_p$conj_poisson)) {
  cp <- ps_p$conj_poisson
  pf_conj <- dGamma(shape = cp$shape, rate = cp$rate, beta = cp$beta, Inv_Dispersion = FALSE)
  ## glmb(n = 500, formula = y ~ 1, data = df_p,
  ##      family = poisson(link = "identity"), pfamily = pf_conj)
}

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