bas.glm {BAS}  R Documentation 
Bayesian Adaptive Sampling Without Replacement for Variable Selection in Generalized Linear Models
Description
Sample with or without replacement from a posterior distribution on GLMs
Usage
bas.glm(
formula,
family = binomial(link = "logit"),
data,
weights,
subset,
contrasts = NULL,
offset,
na.action = "na.omit",
n.models = NULL,
betaprior = CCH(alpha = 0.5, beta = as.numeric(nrow(data)), s = 0),
modelprior = beta.binomial(1, 1),
initprobs = "Uniform",
include.always = ~1,
method = "MCMC",
update = NULL,
bestmodel = NULL,
prob.rw = 0.5,
MCMC.iterations = NULL,
thin = 1,
control = glm.control(),
laplace = FALSE,
renormalize = FALSE,
force.heredity = FALSE,
bigmem = FALSE
)
Arguments
formula 
generalized linear model formula for the full model with all predictors, Y ~ X. All code assumes that an intercept will be included in each model. 
family 
a description of the error distribution and link function for exponential family; currently only 'binomial()' with the logistic link and 'poisson()' and 'Gamma()'with the log link are available. 
data 
data frame 
weights 
optional vector of weights to be used in the fitting process. May be missing in which case weights are 1. 
subset 
subset of data used in fitting 
contrasts 
an optional list. See the contrasts.arg of 'model.matrix.default()'. 
offset 
a priori known component to be included in the linear predictor; by default 0. 
na.action 
a function which indicates what should happen when the data contain NAs. The default is "na.omit". 
n.models 
number of unique models to keep. If NULL, BAS will attempt to enumerate unless p > 35 or method="MCMC". For any of methods using MCMC algorithms that sample with replacement, sampling will stop when the number of iterations exceeds the min of 'n.models' or 'MCMC.iterations' and on exit 'n.models' is updated to reflect the unique number of models that have been sampled. 
betaprior 
Prior on coefficients for model coefficients (except
intercept). Options include

modelprior 
Family of prior distribution on the models. Choices
include 
initprobs 
vector of length p with the initial inclusion probabilities
used for sampling without replacement (the intercept will be included with
probability one and does not need to be added here) or a character string
giving the method used to construct the sampling probabilities if "Uniform"
each predictor variable is equally likely to be sampled (equivalent to
random sampling without replacement). If "eplogp", use the

include.always 
A formula with terms that should always be included in the model with probability one. By default this is '~ 1' meaning that the intercept is always included. This will also override any of the values in 'initprobs' above by setting them to 1. 
method 
A character variable indicating which sampling method to use: method="BAS" uses Bayesian Adaptive Sampling (without replacement) using the sampling probabilities given in initprobs and updates using the marginal inclusion probabilities to direct the search/sample; method="MCMC" combines a random walk Metropolis Hastings (as in MC3 of Raftery et al 1997) with a random swap of a variable included with a variable that is currently excluded (see Clyde, Ghosh, and Littman (2010) for details); method="MCMC+BAS" runs an initial MCMC as above to calculate marginal inclusion probabilities and then samples without replacement as in BAS; method = "deterministic" runs an deterministic sampling using the initial probabilities (no updating); this is recommended for fast enumeration or if a model of independence is a good approximation to the joint posterior distribution of the model indicators. For BAS, the sampling probabilities can be updated as more models are sampled. (see 'update' below). We recommend "MCMC+BAS" or "MCMC" for high dimensional problems. 
update 
number of iterations between potential updates of the sampling probabilities in the "BAS" method. If NULL do not update, otherwise the algorithm will update using the marginal inclusion probabilities as they change while sampling takes place. For large model spaces, updating is recommended. If the model space will be enumerated, leave at the default. 
bestmodel 
optional binary vector representing a model to initialize the sampling. If NULL sampling starts with the null model 
prob.rw 
For any of the MCMC methods, probability of using the randomwalk proposal; otherwise use a random "flip" move to propose a new model. 
MCMC.iterations 
Number of models to sample when using any of the MCMC options; should be greater than 'n.models'. By default 10*n.models. 
thin 
oFr "MCMC", thin the MCMC chain every "thin" iterations; default is no thinning. For large p, thinning can be used to significantly reduce memory requirements as models and associated summaries are saved only every thin iterations. For thin = p, the model and associated output are recorded every p iterations,similar to the Gibbs sampler in SSVS. 
control 
a list of parameters that control convergence in the fitting
process. See the documentation for 
laplace 
logical variable for whether to use a Laplace approximate for integration with respect to g to obtain the marginal likelihood. If FALSE the Cephes library is used which may be inaccurate for large n or large values of the Wald Chisquared statistic. 
renormalize 
logical variable for whether posterior probabilities should be based on renormalizing marginal likelihoods times prior probabilities or use Monte Carlo frequencies. Applies only to MCMC sampling. 
force.heredity 
Logical variable to force all levels of a factor to be included together and to include higher order interactions only if lower order terms are included. Currently only supported with ‘method=’MCMC'' and ‘method=’BAS'' (experimental) on nonSolaris platforms. Default is FALSE. 
bigmem 
Logical variable to indicate that there is access to large amounts of memory (physical or virtual) for enumeration with large model spaces, e.g. > 2^25. 
Details
BAS provides several search algorithms to find high probability models for
use in Bayesian Model Averaging or Bayesian model selection. For p less than
2025, BAS can enumerate all models depending on memory availability, for
larger p, BAS samples without replacement using random or deterministic
sampling. The Bayesian Adaptive Sampling algorithm of Clyde, Ghosh, Littman
(2010) samples models without replacement using the initial sampling
probabilities, and will optionally update the sampling probabilities every
"update" models using the estimated marginal inclusion probabilities. BAS
uses different methods to obtain the initprobs
, which may impact the
results in highdimensional problems. The deterministic sampler provides a
list of the top models in order of an approximation of independence using
the provided initprobs
. This may be effective after running the
other algorithms to identify high probability models and works well if the
correlations of variables are small to modest. The priors on coefficients
are mixtures of gpriors that provide approximations to the power prior.
Value
bas.glm
returns an object of class basglm
An object of class basglm
is a list containing at least the following
components:
postprobs 
the posterior probabilities of the models selected 
priorprobs 
the prior probabilities of the models selected 
logmarg 
values of the log of the marginal likelihood for the models 
n.vars 
total number of independent variables in the full model, including the intercept 
size 
the number of independent variables in each of the models, includes the intercept 
which 
a list of lists with one list per model with variables that are included in the model 
probne0 
the posterior probability that each variable is nonzero 
mle 
list of lists with one list per model giving the GLM estimate of each (nonzero) coefficient for each model. 
mle.se 
list of lists with one list per model giving the GLM standard error of each coefficient for each model 
deviance 
the GLM deviance for each model 
modelprior 
the prior distribution on models that created the BMA object 
Q 
the Q statistic for each model used in the marginal likelihood approximation 
Y 
response 
X 
matrix of predictors 
family 
family object from the original call 
betaprior 
family object for prior on coefficients, including hyperparameters 
modelprior 
family object for prior on the models 
include.always 
indices of variables that are forced into the model 
Author(s)
Merlise Clyde (clyde@duke.edu), Quanli Wang and Yingbo Li
References
Li, Y. and Clyde, M. (2018) Mixtures of gpriors in Generalized
Linear Models.
Journal of the American Statistical Association. 113:18281845
doi:10.1080/01621459.2018.1469992
Clyde, M. Ghosh, J. and Littman, M. (2010) Bayesian Adaptive Sampling for
Variable Selection and Model Averaging. Journal of Computational Graphics
and Statistics. 20:80101
doi:10.1198/jcgs.2010.09049
Raftery, A.E, Madigan, D. and Hoeting, J.A. (1997) Bayesian Model Averaging
for Linear Regression Models. Journal of the American Statistical
Association.
Examples
library(MASS)
data(Pima.tr)
# enumeration with default method="BAS"
pima.cch = bas.glm(type ~ ., data=Pima.tr, n.models= 2^7,
method="BAS",
betaprior=CCH(a=1, b=532/2, s=0), family=binomial(),
modelprior=beta.binomial(1,1))
summary(pima.cch)
image(pima.cch)
# Note MCMC.iterations are set to 2500 for illustration purposes due to time
# limitations for running examples on CRAN servers.
# Please check convergence diagnostics and run longer in practice
pima.robust = bas.glm(type ~ ., data=Pima.tr, n.models= 2^7,
method="MCMC", MCMC.iterations=2500,
betaprior=robust(), family=binomial(),
modelprior=beta.binomial(1,1))
pima.BIC = bas.glm(type ~ ., data=Pima.tr, n.models= 2^7,
method="BAS+MCMC", MCMC.iterations=2500,
betaprior=bic.prior(), family=binomial(),
modelprior=uniform())
# Poisson example
if(requireNamespace("glmbb", quietly=TRUE)) {
data(crabs, package='glmbb')
#short run for illustration
crabs.bas = bas.glm(satell ~ color*spine*width + weight, data=crabs,
family=poisson(),
betaprior=EB.local(), modelprior=uniform(),
method='MCMC', n.models=2^10, MCMC.iterations=2500,
prob.rw=.95)
# Gamma example
if(requireNamespace("faraway", quietly=TRUE)) {
data(wafer, package='faraway')
wafer_bas = bas.glm(resist~ ., data=wafer, include.always = ~ .,
betaprior = bic.prior() ,
family = Gamma(link = "log"))
}
}