R: Bayesian Variable Selection for linear regression models...

GibbsBvs {BayesVarSel}

R Documentation

Bayesian Variable Selection for linear regression models using Gibbs sampling.

Description

Approximate computation of summaries of the posterior distribution using a Gibbs sampling algorithm to explore the model space and frequency of "visits" to construct the estimates.

Usage

GibbsBvs(
  formula,
  data,
  null.model = paste(as.formula(formula)[[2]], " ~ 1", sep = ""),
  prior.betas = "Robust",
  prior.models = "ScottBerger",
  n.iter = 10000,
  init.model = "Full",
  n.burnin = 500,
  n.thin = 1,
  time.test = TRUE,
  priorprobs = NULL,
  seed = runif(1, 0, 16091956)
)

Arguments

`formula`	Formula defining the most complex regression model in the analysis. See details.
`data`	data frame containing the data.
`null.model`	A formula defining which is the simplest (null) model. It should be nested in the full model. By default, the null model is defined to be the one with just the intercept.
`prior.betas`	Prior distribution for regression parameters within each model (to be literally specified). Possible choices include "Robust", "Robust.G", "Liangetal", "gZellner", "ZellnerSiow", "FLS", "intrinsic.MGC" and "IHG" (see details).
`prior.models`	Prior distribution over the model space (to be literally specified). Possible choices are "Constant", "ScottBerger" and "User" (see details).
`n.iter`	The total number of iterations performed after the burn in process.
`init.model`	The model at which the simulation process starts. Options include "Null" (the model only with the covariates specified in `fixed.cov`), "Full" (the model defined by `formula`), "Random" (a randomly selected model) and a vector with p (the number of covariates to select from) zeros and ones defining a model. When p>n the dimension of the init.model must be smaller than n. Otherwise the function produces an error.
`n.burnin`	Length of burn in, i.e. number of iterations to discard at the beginning.
`n.thin`	Thinning rate. Must be a positive integer. Set 'n.thin' > 1 to save memory and computation time if 'n.iter' is large. Default is 1. This parameter jointly with `n.iter` sets the number of simulations kept and used to construct the estimates so is important to keep in mind that a large value for 'n.thin' can reduce the precision of the results
`time.test`	If TRUE and the number of variables is large (>=21) a preliminary test to estimate computational time is performed.
`priorprobs`	A p+1 dimensional vector defining the prior probabilities Pr(M_i) (should be used in the case where `prior.models`="User"; see the details in `Bvs`.)
`seed`	A seed to initialize the random number generator

Details

This is a heuristic approximation to the function Bvs so the details there apply also here.

The algorithm implemented is a Gibbs sampling-based searching algorithm originally proposed by George and McCulloch (1997). Garcia-Donato and Martinez-Beneito (2013) have shown that this simple sampling strategy in combination with estimates based on frequency of visits (the one here implemented) provides very reliable results.

Value

GibbsBvs returns an object of class Bvs with the following elements:

`time`	The internal time consumed in solving the problem
`lmfull`	The `lm` class object that results when the model defined by `formula` is fitted by `lm`
`lmnull`	The `lm` class object that results when the model defined by `fixed.cov` is fitted by `lm`
`variables`	The name of all the potential explanatory variables
`n`	Number of observations
`p`	Number of explanatory variables to select from
`k`	Number of fixed variables
`HPMbin`	The binary expression of the most probable model found.
`inclprob`	A named vector with the estimates of the inclusion probabilities of all the variables.
`jointinclprob`	A `data.frame` with the estimates of the joint inclusion probabilities of all the variables.
`postprobdim`	Estimates of posterior probabilities of the dimension of the true model.
`modelslogBF`	A matrix with both the binary representation of the visited models after the burning period and the Bayes factor (log scale) of that model to the null model.
`priorprobs`	If `prior.models`="User" then this vector is stored here. Else, the #' type of prior as defined in `prior.models`
`call`	The `call` to the function.
`C`	An estimation of the normalizing constant (C=sum Bi Pr(Mi), for Mi in the model space) using the method in George and McCulloch (1997).
`method`	`gibbs`
`prior.betas`	prior.betas
`prior.models`	prior.models
`priorprobs`	priorprobs

Author(s)

Gonzalo Garcia-Donato and Anabel Forte

References

Garcia-Donato, G. and Martinez-Beneito, M.A. (2013)<DOI:10.1080/01621459.2012.742443> On sampling strategies in Bayesian variable selection problems with large model spaces. Journal of the American Statistical Association, 108: 340-352.

George E. and McCulloch R. (1997) Approaches for Bayesian variable selection. Statistica Sinica, 7, 339:372.

Examples


## Not run: 
#Analysis of Ozone35 data

data(Ozone35)

#We use here the (Zellner) g-prior for
#regression parameters and constant prior
#over the model space
#In this Gibbs sampling scheme, we perform 10100 iterations,
#of which the first 100 are discharged (burnin) and of the remaining
#only one each 10 is kept.
#as initial model we use the Full model
Oz35.GibbsBvs<- GibbsBvs(formula= y ~ ., data=Ozone35, prior.betas="gZellner",
prior.models="Constant", n.iter=10000, init.model="Full", n.burnin=100,
time.test = FALSE)

#Note: this is a heuristic approach and results are estimates
#of the exact answer.

#with the print we can see which is the most probable model
#among the visited
Oz35.GibbsBvs

#The estimation of inclusion probabilities and
#the model-averaged estimation of parameters:
summary(Oz35.GibbsBvs)

#Plots:
plot(Oz35.GibbsBvs, option="conditional")

## End(Not run)

[Package BayesVarSel version 2.2.5 Index]