Sample with or without replacement from a posterior distribution on GLMs
bas.glm(formula, family = binomial(link = "logit"), data, weights, subset,
offset, na.action = "na.omit", n.models = NULL, betaprior = CCH(alpha =
0.5, beta = nrow(data), s = 0), modelprior = beta.binomial(1, 1),
initprobs = "Uniform", method = "MCMC", update = NULL,
bestmodel = NULL, prob.rw = 0.5, MCMC.iterations = NULL,
control = glm.control(), laplace = FALSE, renormalize = FALSE)
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.
a description of the error distribution and link function for exponential family; currently only binomial() with the logitistic link and poission() with the log link are available.
data frame
optional vector of weights to be used in the fitting process. May be missing in which case weights are 1.
subset of data used in fitting
a priori known component to be included in the linear predictor; by default 0.
a function which indicates what should happen when the data contain NAs. The default is "na.omit".
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.
Family of prior distribution on the models. Choices
include uniform
, Bernoulli
,
beta.binomial
, truncated Beta-Binomial,
tr.beta.binomial
, and truncated power family
tr.power.prior
.
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
eplogprob
function to aproximate the Bayes factor using
p-values to find initial marginal inclusion probabilitites and sample
without replacement using these inclusion probabilaties, which may be
updated using estimates of the marginal inclusion probabilites. "eplogp"
assumes that MLEs from the full model exist; for problems where that is not
the case or 'p' is large, initial sampling probabilities may be obtained
using eplogprob.marg
which fits a model to each predictor
seaparately. For variables that should always be included set the
corresponding initprobs to 1. To run a Markov Chain to provide initial
estimates of marginal inclusion probabilities, use method="MCMC+BAS" below.
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 probabilites (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.
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.
optional binary vector representing a model to initialize the sampling. If NULL sampling starts with the null model
For any of the MCMC methods, probability of using the random-walk proposal; otherwise use a random "flip" move to propose a new model.
Number of models to sample when using any of the MCMC options; should be greater than 'n.models'.
a list of parameters that control convergence in the fitting
process. See the documentation for glm.control()
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.
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
bas.glm
returns an object of class basglm
An object of class basglm
is a list containing at least the following
components:
the posterior probabilities of the models selected
the prior probabilities of the models selected
values of the log of the marginal likelihood for the models
total number of independent variables in the full model, including the intercept
the number of independent variables in each of the models, includes the intercept
a list of lists with one list per model with variables that are included in the model
the posterior probability that each variable is non-zero
list of lists with one list per model giving the GLM estimate of each (nonzero) coefficient for each model.
list of lists with one list per model giving the GLM standard error of each coefficient for each model
the GLM deviance for each model
the prior distribution on models that created the BMA object
the Q statistic for each model used in the marginal likelihood approximation
response
matrix of predictors
family object from the original call
family object for prior on coefficients, including hyperparamters
family object for prior on the models
BAS provides several search algorithms to find high probability models for
use in Bayesian Model Averaging or Bayesian model selection. For p less than
20-25, 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 probabilties. BAS
uses different methods to obtain the initprobs
, which may impact the
results in high-dimensional problems. The deterinistic 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 g-priors that provide approximations to the power prior.
Li, Y. and Clyde, M. (2015) Mixtures of g-priors in Generalized Linear Models. http://arxiv.org/abs/1503.06913
Clyde, M. Ghosh, J. and Littman, M. (2010) Bayesian Adaptive Sampling for Variable Selection and Model Averaging. Journal of Computational Graphics and Statistics. 20:80-101 http://dx.doi.org/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.
# NOT RUN {
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)
pima.robust = bas.glm(type ~ ., data=Pima.tr, n.models= 2^7,
method="MCMC", MCMC.iterations=20000,
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=1000,
betaprior=bic.prior(), family=binomial(),
modelprior=uniform())
# }
Run the code above in your browser using DataLab