Package {IBGS}

Iterated Block Gibbs Sampler for Ultrahigh-Dimensional Variable Selection and Model Averaging

Description

Variable selection for generalized linear models, the Cox proportional-hazards model and the linear mixed model in ultrahigh dimensions via the iterated block Gibbs sampler (IBGS). The sampler is implemented in C with parallel block screening through OpenMP, and supports the gaussian, binomial and poisson families (fitted by least squares or iteratively reweighted least squares), the Cox model for survival analysis (fitted by its Efron partial likelihood), and the fixed-effect part of a linear mixed model (fitted on whitened data). Candidate models are compared by the AIC, BIC, AICc or extended BIC criterion.

Details

Each model type comes with two samplers: a block sampler that splits the predictors into small groups for ultrahigh dimensions, and a standalone (non-block) sampler for problems where the full design is manageable.

Every sampler returns an object of class "IBGS" with print, summary, coef and plot methods. predict and fitted form model-averaged predictions, and plotICtrace, plotModelFreq and plotMargProb provide diagnostic plots.

Author(s)

Lizhong Chen <chen.l@wehi.edu.au>

Maintainer: Lizhong Chen <chen.l@wehi.edu.au>

References

Chen, L. (2022). Model selection and averaging by Gibbs sampler with a tropical cyclone seasonal forecasting application. PhD thesis, University of Melbourne. https://hdl.handle.net/11343/311691

See Also

glmIBGS, coxIBGS, lmeIBGS, predict.IBGS

Examples

x <- matrix(rnorm(100 * 50), 100, 50)
y <- rowSums(x[, 1:3]) + rnorm(100)
fit <- glmIBGS(y, x, criterion = "BIC")
fit
predict(fit, x[1:5, ])

Coefficients of an IBGS fit

Description

Returns the coefficient vector of one of the retained top models, or the model-averaged coefficient vector (smooth-SIC weights over the best n.models, the same weighting predict.IBGS uses). The vector is named by the predictors (with "(Intercept)" first for the families that carry an intercept).

Usage

## S3 method for class 'IBGS'
coef(object, model = 1, average = FALSE, n.models = object$n.models, ...)

Arguments

Value

a named numeric vector of coefficients

The restricted Gibbs sampler for the Cox model

Description

The standalone (non-block) Metropolis-within-Gibbs sampler for the Cox proportional-hazards model that searches the model space directly over all predictors, with the model size capped at max.size. This is the survival-analysis counterpart of glmGibbs: every candidate Cox model is fitted by its Efron partial likelihood in a fast C routine (cox_gibbs_glm), and the summary is assembled by refitting the top models in C. Intended for problems where the full design is manageable; use coxIBGS in ultrahigh dimensions.

Usage

coxGibbs(
  y,
  status,
  x,
  max.size = ncol(x),
  permute = TRUE,
  n.models = 10,
  threshold = 0.9,
  n.draws = 1000,
  inv.temp = 0.5,
  ebic.gamma = 0.5,
  criterion = c("AIC", "BIC", "AICc", "exBIC"),
  weights = NULL,
  cor.check = NULL,
  start = c("null", "full")
)

Arguments

Details

coxGibbs runs a single random-permutation Metropolis-within-Gibbs walk over the inclusion indicators of all p predictors, with the model size capped at max.size and no block screening, just as glmGibbs does for the GLM families. The acceptance rule min(1, exp(inv.temp * (IC.current - IC.proposed))), the four information criteria and the model averaging over the top n.models models are as described under glmIBGS; each candidate is fitted by its Efron partial likelihood, the model has no intercept, and the criteria use the number of events as the effective sample size. Use coxIBGS in ultrahigh dimensions.

Value

an object of class "IBGS" summarising the search, with components marginal.prob (the marginal inclusion probability of each predictor), selected.vars (the names of the predictors above threshold), var.names, coef (the coefficient matrix of the top n.models models), model.ic and model.freq (the criterion and visit frequency of each top model), ic.trace (the criterion at every generation) and criterion (the criterion name). Has print, summary, coef and plot methods.

Examples

x <- matrix(rnorm(2000), ncol = 10)
eta <- x[, 1:3] %*% c(1, -1, 0.5)
time <- rexp(200, exp(eta))
status <- rbinom(200, 1, 0.8)
m.s <- coxGibbs(time, status, x, criterion = "BIC")

The iterated block Gibbs sampler for the Cox model

Description

Variable selection for the Cox proportional-hazards model in ultrahigh dimensions. This is the survival-analysis counterpart of glmIBGS: the whole search (refinement iterations, random block screening and the final long run) is performed in a single multicore C routine (cox_ibgs_glm), and this R function is a thin wrapper that assembles the summary (refitting the top models in C). Each candidate model is fitted by its Efron partial likelihood; the Cox model has no intercept, and the information criteria use the number of events as the effective sample size.

Usage

coxIBGS(
  y,
  status,
  x,
  n.refine = 3,
  n.models = 10,
  block.size = 30,
  n.keep = 20,
  threshold = 0.9,
  permute = TRUE,
  n.draws = 250,
  inv.temp = 0.5,
  ebic.gamma = 0.5,
  criterion = c("AIC", "BIC", "AICc", "exBIC"),
  weights = NULL,
  n.cores = 1L,
  cor.check = NULL,
  start = c("null", "full")
)

Arguments

Details

coxIBGS applies the same iterated block Gibbs algorithm as glmIBGS – the screen / select / threshold refinement loop, the final long Gibbs run, the Metropolis-within-Gibbs acceptance rule min(1, exp(inv.temp * (IC.current - IC.proposed))), the four information criteria and the model averaging over the top n.models models – but for the Cox proportional-hazards model rather than a GLM. Each candidate model is fitted in C by its Efron partial likelihood; the Cox model carries no intercept, and the information criteria use the number of events (not the sample size n) as the effective sample size. See glmIBGS for the screening, sampler and criteria details.

Value

an object of class "IBGS" summarising the search, with components marginal.prob (the marginal inclusion probability of each predictor), selected.vars (the names of the predictors above threshold), var.names, coef (the coefficient matrix of the top n.models models), model.ic and model.freq (the criterion and visit frequency of each top model), ic.trace (the criterion at every generation) and criterion (the criterion name). Has print, summary, coef and plot methods.

Examples

x <- matrix(rnorm(60 * 20), 60, 20)
time <- rexp(60, exp(rowSums(x[, 1:3])))
status <- rbinom(60, 1, 0.7)
fit <- coxIBGS(time, status, x, criterion = "BIC")
fit

Fitted values of an IBGS fit

Description

The model-averaged fitted values on the training data, on the response (default) or link scale. The averaged training linear predictor (the linear.predictors component) is computed once at fit time and stored on the object, so this is cheap even in high dimensions. Equivalent to predict(object) with no new data.

Usage

## S3 method for class 'IBGS'
fitted(object, type = c("response", "link"), ...)

Arguments

Value

a numeric vector of model-averaged fitted values on the training data (length equal to the number of training observations).

See Also

predict.IBGS, coef.IBGS

Examples

x <- matrix(rnorm(1000), ncol = 10)
y <- rowSums(x[, 1:5]) + rnorm(100)
fit <- glmGibbs(y, x, criterion = "BIC")
head(fitted(fit))

The restricted Gibbs sampler (linear model)

Description

The standalone (non-block) Metropolis-within-Gibbs sampler that searches the model space directly over all predictors, with the model size capped at max.size. Every candidate linear model is fitted by ordinary least squares in a fast C routine (gibbs_sampler_gaussian); the summary is assembled by refitting the top models in C to recover their coefficients. Unlike glmIBGS, this sampler does not split the predictors into blocks and is intended for problems where the full design is manageable.

Usage

glmGibbs(
  y,
  x,
  max.size = ncol(x),
  permute = TRUE,
  n.models = 10,
  threshold = 0.9,
  n.draws = 1000,
  inv.temp = 0.5,
  ebic.gamma = 0.5,
  criterion = c("AIC", "BIC", "AICc", "exBIC"),
  family = c("gaussian", "binomial", "poisson"),
  weights = NULL,
  cor.check = NULL,
  start = c("null", "full"),
  fast = TRUE
)

Arguments

Details

glmGibbs runs a single random-permutation Metropolis-within-Gibbs walk directly over the inclusion indicators of all p predictors, with the model size capped at max.size; there is no block screening or refinement. The flip acceptance rule min(1, exp(inv.temp * (IC.current - IC.proposed))), the four information criteria (AIC, BIC, AICc and exBIC) and the model averaging over the top n.models models are exactly as described under glmIBGS, and each candidate is fitted in the same way (gaussian by ordinary least squares, binomial and poisson by IRLS). Because every candidate fit may use all p predictors, this sampler is intended for problems where the full design is manageable (p not much larger than n); use glmIBGS in ultrahigh dimensions.

Value

an object of class "IBGS" summarising the search, with components marginal.prob (the marginal inclusion probability of each predictor), selected.vars (the names of the predictors above threshold), var.names, coef (the coefficient matrix of the top n.models models), model.ic and model.freq (the criterion and visit frequency of each top model), ic.trace (the criterion at every generation) and criterion (the criterion name). Has print, summary, coef and plot methods.

Examples

x <- matrix(rnorm(1000), ncol = 10)
y <- rowSums(x[,1:5]) + rnorm(100)
m.s <- glmGibbs(y, x, criterion = "BIC")

The iterated block Gibbs sampler algorithm

Description

Variable selection for the generalized linear model in ultrahigh dimensions. The whole search – the refinement iterations, the random block screening, and the final long run – is performed in a single multicore C routine (ibgs_glm); this R function is a thin wrapper that assembles the result summary (refitting the top models in C to recover their coefficients). The gaussian family is fitted by ordinary least squares; the binomial and poisson families by IRLS.

Usage

glmIBGS(
  y,
  x,
  n.refine = 3,
  n.models = 10,
  block.size = 30,
  n.keep = 20,
  threshold = 0.9,
  permute = TRUE,
  n.draws = 250,
  inv.temp = 0.5,
  ebic.gamma = 0.5,
  criterion = c("AIC", "BIC", "AICc", "exBIC"),
  family = c("gaussian", "binomial", "poisson"),
  weights = NULL,
  n.cores = 1L,
  cor.check = NULL,
  start = c("null", "full"),
  fast = TRUE
)

Arguments

Details

In ultrahigh dimensions the number of predictors p can far exceed the sample size n, so a single Gibbs sampler over all p inclusion indicators is both slow (each candidate fit can use at most n columns) and noisy. IBGS instead alternates screening and refinement to grow a small, stable set of important predictors:

• Screen. Randomly partition the not-yet-important predictors into blocks of size block.size (capped at n minus the current important set, so every block design is estimable). An independent within-block Gibbs sampler is run on each block – always including the important set – giving a marginal inclusion probability for every predictor in the block.

• Select. Keep the top n.keep predictors by inclusion probability, union them with the important set, and run one combined Gibbs sampler over that union.

• Threshold. Promote every predictor whose combined inclusion probability exceeds threshold into the important set.

These steps are repeated for n.refine refinement rounds; a final screen is then followed by a long Gibbs run over the final candidate set, which produces the recorded indicator samples, their information criteria, and the marginal inclusion probabilities returned to R.

Each Gibbs sampler is a random-permutation Metropolis-within-Gibbs walk over the model indicator vector: each sweep visits every coordinate once in a fresh random order, flipping it and accepting the flip with probability min(1, exp(inv.temp * (IC.current - IC.proposed))). Because a lower criterion is better, a flip that lowers the criterion is always accepted and a worsening flip is accepted with a probability that decays in the increase; the tuning parameter inv.temp controls how greedy the walk is. The chain's stationary distribution favours low-criterion models and its visit frequencies estimate the marginal inclusion probabilities. Every candidate model is fitted in C – the gaussian family by ordinary least squares (a Cholesky solve on the active sub-block of a precomputed Gram matrix), the binomial and poisson families by iteratively reweighted least squares.

Models are compared by one of four information criteria, each of the form -2logLik + penalty: AIC = -2logLik + 2 npar, BIC = -2logLik + log(n) npar, AICc = AIC + 2 npar (npar + 1) / (n - npar - 1) (a small-sample AIC), and exBIC = BIC + 2 ebic.gamma npred log(p0) (the extended BIC of Chen and Chen, 2008, whose extra term penalises the size p0 of the candidate pool and is essential in ultrahigh dimensions where BIC under-penalises). Here npar counts the estimated parameters and npred the predictors excluding the intercept.

The returned object summarises the search: the top n.models distinct models are refit in C to recover their coefficients, and averaged coefficients and predictions combine them with smooth-criterion weights (see coef.IBGS and predict.IBGS). The block screening is parallelised with n.cores OpenMP threads (one thread per block); the random block assignment and the per-block seeds are drawn serially from R's generator, so the result is reproducible and does not depend on the number of threads.

Value

an object of class "IBGS" summarising the search, with components marginal.prob (the marginal inclusion probability of each predictor), selected.vars (the names of the predictors above threshold), var.names, coef (the coefficient matrix of the top n.models models), model.ic and model.freq (the criterion and visit frequency of each top model), ic.trace (the criterion at every generation) and criterion (the criterion name). Has print, summary, coef and plot methods.

Examples

x <- matrix(rnorm(100 * 50), 100, 50)
y <- rowSums(x[, 1:3]) + rnorm(100)
fit <- glmIBGS(y, x, criterion = "BIC")
fit

The restricted Gibbs sampler for the linear mixed model (fixed-effect part)

Description

The standalone (non-block) Metropolis-within-Gibbs sampler that selects the fixed effects of a linear mixed model while holding the random part fixed. With the random structure held, the marginal model is y ~ N(X beta, V) with V fixed; whitening by V = L Lt reduces each candidate fit to ordinary least squares on the whitened data, performed in a fast C routine. The variance component(s) are estimated once (random-intercept case) or supplied, and held fixed across all candidate models. This is the mixed-model counterpart of glmGibbs; use lmeIBGS in ultrahigh dimensions.

Usage

lmeGibbs(
  y,
  x,
  group = NULL,
  Z = NULL,
  varcomp = NULL,
  V = NULL,
  max.size = ncol(x),
  permute = TRUE,
  n.models = 10,
  threshold = 0.9,
  n.draws = 1000,
  inv.temp = 0.5,
  ebic.gamma = 0.5,
  criterion = c("AIC", "BIC", "AICc", "exBIC"),
  cor.check = NULL,
  start = c("null", "full")
)

Arguments

Details

lmeGibbs runs a single random-permutation Metropolis-within-Gibbs walk over the inclusion indicators of all p fixed-effect predictors, with the model size capped at max.size and no block screening, just as glmGibbs does for the GLM families – but on whitened data. With the random part held fixed the marginal model is y ~ N(X beta, V) with V fixed; whitening by V = L Lt reduces each candidate fit to ordinary least squares. The acceptance rule min(1, exp(inv.temp * (IC.current - IC.proposed))), the four information criteria and the model averaging over the top n.models models are as described under glmIBGS; the variance component(s) are estimated once (random-intercept case, by REML) or supplied, and held fixed. Use lmeIBGS in ultrahigh dimensions.

Specify the random part with exactly one of: group (a grouping factor for a random intercept; its variance ratio is estimated once by REML), Z + varcomp, or a marginal covariance V.

Value

an object of class "IBGS" summarising the search, with components marginal.prob (the marginal inclusion probability of each predictor), selected.vars (the names of the predictors above threshold), var.names, coef (the coefficient matrix of the top n.models models), model.ic and model.freq (the criterion and visit frequency of each top model), ic.trace (the criterion at every generation), criterion (the criterion name) and, when the random part is given via group or Z + varcomp, re (the random-effect BLUPs used for conditional prediction; a directly supplied V gives the fixed-effect fit without re). Has print, summary, coef and plot methods.

Examples

x <- matrix(rnorm(1000), ncol = 10)
g <- gl(10, 10)
y <- rowSums(x[, 1:5]) + rnorm(10)[g] + rnorm(100)
fit <- lmeGibbs(y, x, group = g, criterion = "BIC")

The iterated block Gibbs sampler for the linear mixed model (fixed-effect part)

Description

Fixed-effect variable selection for a linear mixed model in ultrahigh dimensions, holding the random part fixed. This is the mixed-model counterpart of glmIBGS: the whole search (refinement iterations, random block screening and the final long run) runs in a single multicore C routine on the whitened data. With the random structure held the marginal model is y ~ N(X beta, V) with V fixed; whitening by V = L Lt reduces each candidate fit to ordinary least squares. The variance component(s) are estimated once (random-intercept case, by REML from the largest estimable fixed model) or supplied, and held fixed.

Usage

lmeIBGS(
  y,
  x,
  group = NULL,
  Z = NULL,
  varcomp = NULL,
  V = NULL,
  n.refine = 3,
  n.models = 10,
  block.size = 30,
  n.keep = 20,
  threshold = 0.9,
  permute = TRUE,
  n.draws = 250,
  inv.temp = 0.5,
  ebic.gamma = 0.5,
  criterion = c("AIC", "BIC", "AICc", "exBIC"),
  n.cores = 1L,
  cor.check = NULL,
  start = c("null", "full")
)

Arguments

Details

lmeIBGS runs the same iterated block Gibbs algorithm as glmIBGS – the screen / select / threshold refinement loop, the final long Gibbs run, the acceptance rule min(1, exp(inv.temp * (IC.current - IC.proposed))), the four information criteria and the model averaging over the top n.models models – but on whitened data. With the random part held fixed the marginal model is y ~ N(X beta, V) with V fixed; whitening by V = L Lt (left-multiplying y and X by the inverse root) reduces each candidate fit to ordinary least squares, so the C sampler is the gaussian one applied to the whitened design. The variance component(s) are estimated once – for a random intercept by REML from the largest estimable fixed model – or supplied directly, and then held fixed across all candidate models. See glmIBGS for the screening, sampler and criteria details.

Specify the random part with exactly one of: group (a grouping factor for a random intercept), Z + varcomp, or a marginal covariance V.

Value

an object of class "IBGS" summarising the search, with components marginal.prob (the marginal inclusion probability of each predictor), selected.vars (the names of the predictors above threshold), var.names, coef (the coefficient matrix of the top n.models models), model.ic and model.freq (the criterion and visit frequency of each top model), ic.trace (the criterion at every generation), criterion (the criterion name) and, when the random part is given via group or Z + varcomp, re (the random-effect BLUPs used for conditional prediction; a directly supplied V gives the fixed-effect fit without re). Has print, summary, coef and plot methods.

Examples

x <- matrix(rnorm(100 * 50), 100, 50)
g <- gl(10, 10)
y <- rowSums(x[, 1:3]) + rnorm(10)[g] + rnorm(100)
fit <- lmeIBGS(y, x, group = g, criterion = "BIC")
fit

Diagnostic plots for an IBGS fit

Description

Dispatches to the diagnostic plots: the trace of the criterion sequence (plotICtrace), the marginal inclusion probabilities (plotMargProb), the top-model visit frequencies (plotModelFreq), the Gelman-Rubin shrink factor (plotGelman), and the criterion-trace autocorrelation (plotAutocorr). Several may be requested at once, in which case they are drawn in a multi-panel layout.

Usage

## S3 method for class 'IBGS'
plot(
  x,
  which = c("ictrace", "margprob", "modelfreq", "gelman", "autocorr"),
  ...
)

Arguments

Value

x, invisibly

Autocorrelation plot for the criterion trace

Description

Draws the lagged autocorrelation of the criterion trace (result$convergence$autocorr) as a stick plot. A trace that mixes well decays to zero quickly; persistent high autocorrelation signals a slowly-mixing chain and a small effective sample size.

Usage

plotAutocorr(result, col = NULL, ...)

Arguments

Details

The autocorrelations are computed in C with no dependency on the coda package, matching stats::acf applied to the criterion trace.

Value

invisibly, the named vector of autocorrelations (named by lag)

Gelman-Rubin shrink-factor plot for the criterion trace

Description

Plots the evolving Gelman-Rubin potential scale reduction factor (R-hat) of the criterion trace as the run lengthens, computed from the single chain split into equal contiguous segments (result$convergence$shrink). The median and the 97.5% upper limit are drawn together with a horizontal guide at 1 (perfect mixing) and at 1.1 (the usual convergence rule of thumb): once both curves settle below about 1.1 the sampler has stabilised.

Usage

plotGelman(result, col = NULL, legend = TRUE, ...)

Arguments

Details

The diagnostic is computed in C with no dependency on the coda package, but follows the same definition as coda::gelman.diag / coda::gelman.plot applied to the criterion trace split into segments (split-R-hat).

Value

invisibly, the data.frame of shrink-factor values

Trace plot of the information-criterion sequence

Description

Plots the model selection criterion at each Gibbs generation. The best (lowest) criterion value reached is highlighted and a horizontal guide marks the best retained model's criterion, so the minimal information criterion stands out.

Usage

plotICtrace(
  result,
  col = NULL,
  highlight = NULL,
  running.min = FALSE,
  legend = TRUE,
  ...
)

Arguments

Value

invisibly, a list with the best value and the generation that reached it

Plot the marginal inclusion probability of the top covariates

Description

Draws a dot-and-whisker plot of the n.vars covariates with the highest marginal inclusion probability, sorted descending. Each point carries an error bar of plus/minus one Monte-Carlo standard error of the inclusion proportion, se = \sqrt{p (1 - p) / ESS}, where ESS is the effective sample size of the criterion trace (result$convergence$ess); the bars are clamped to [0, 1]. The selected predictors (result$selected.vars) are highlighted and the rest are muted, so the selected set stands out. With horizontal = TRUE the predictors run down the y-axis (names read horizontally on the left); otherwise they run along the x-axis (the default).

Usage

plotMargProb(
  result,
  n.vars = 20,
  col = NULL,
  horizontal = FALSE,
  lwd = 2,
  side = 1,
  line = 0.25,
  las = 2,
  cex = 1
)

Arguments

Value

invisibly, the named vector of plotted inclusion probabilities

Plot the visit frequency of the top selected models

Description

Draws a bar for each of the top n.models models, showing how often the sampler settled on it: the relative frequency with which the Gibbs chain visited that model (result$model.freq). The bars are ordered best-first (the model with the lowest information criterion on the left, highlighted), and each bar is annotated with its criterion value. Tall leading bars indicate a search that concentrated sharply on a few models.

Usage

plotModelFreq(
  result,
  n.models = result$n.models,
  col = NULL,
  annotate = TRUE,
  cumulative = FALSE,
  horizontal = FALSE,
  ...
)

Arguments

Value

invisibly, the vector of plotted relative model frequencies

Predict method for an IBGS fit

Description

Model-averaged prediction (the default) or single-model prediction on new data. The top n.models models are averaged on the linear-predictor (link) scale with smooth-SIC weights (proportional to exp(-inv.temp * (IC - min IC))) and the family inverse link is applied once; averaging before the inverse link keeps binomial predictions in the unit interval and poisson/cox predictions positive. With newdata = NULL the model-averaged fitted values stored on the object at fit time are returned.

Usage

## S3 method for class 'IBGS'
predict(object, newdata = NULL, type = c("response", "link"),
  average = TRUE, n.models = object$n.models, model = 1,
  group.new = NULL, Z.new = NULL, ...)

Arguments

Value

a numeric vector of predictions: the response scale for the GLM families (probabilities in [0, 1] for binomial), the relative risk exp(x*beta) for cox, and the fixed-effect mean – optionally plus the random-effect BLUP – for lme.

See Also

fitted.IBGS, coef.IBGS

Examples

x <- matrix(rnorm(1000), ncol = 10)
y <- rowSums(x[, 1:5]) + rnorm(100)
fit <- glmGibbs(y, x, criterion = "BIC")
predict(fit, x[1:5, ])

Print a concise summary of an IBGS fit

Description

Print a concise summary of an IBGS fit

Usage

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

Arguments

Value

x, invisibly

Tabular summary of an IBGS fit

Description

Tabular summary of an IBGS fit

Usage

## S3 method for class 'IBGS'
summary(object, ...)

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

Arguments

Value

invisibly, an object of class "summary.IBGS" with the selected variable table, the top-model table and the convergence diagnostics of the criterion trace (Gelman-Rubin R-hat, Geweke z, effective sample size and lag-1 autocorrelation), which print displays as a "Convergence diagnostics" block