BayesRLasso: Bayesian Reciprocal Regularization

In himelmallick/BayesRecipe: Bayesian Reciprocal Regularization

Bayesian Reciprocal Regularization

Description

Wrapper for Bayesian Reciprocal LASSO using MCMC

Usage

BayesRLasso(
  x,
  y,
  method = "slice",
  lambda.estimate = "AP",
  update.sigma2 = TRUE,
  max.steps = 11000,
  n.burn = 1000,
  n.thin = 1,
  ridge.CV = TRUE,
  a = 0.001,
  b = 0.001,
  posterior.summary.beta = "mean",
  posterior.summary.lambda = "median",
  beta.ci.level = 0.95,
  lambda.ci.level = 0.95,
  seed = 1234,
  na.rm = TRUE
)

Arguments

x	A numeric matrix with standardized predictors in columns and samples in rows.
y	A mean-centered continuous response variable with matching rows with x.
method	If Gibbs sampling is desired, 'SMU' or 'SMN' must be selected as the data augmentation scheme or mixture representation, although not recommended when the data size is as big as hundreds of covariates. By default it uses an ellipitical slice sampler ('slice') which is scalable for large-scale problems.
lambda.estimate	Estimating lambda by empirical bayes ('EB'), MCMC ('MCMC'), or apriori ('AP'). Default is 'AP' (currently only available option when method = 'slice').
update.sigma2	Whether sigma2 should be updated. Default is TRUE (currently only available option when method = 'slice').
max.steps	Number of MCMC iterations. Default is 11000.
n.burn	Number of burn-in iterations. Default is 1000.
n.thin	Lag at which thinning should be done. Default is 1 (no thinning).
ridge.CV	When method = 'SMU' and X is rank deficient, a ridge parameter is added to the diagonal of the crossproduct (XtX) to allow for proper calculation of the inverse. If TRUE, the ridge parameter is estimated by cross-validation using glmnet. Otherwise, it falls back to adding a small number (1e-05) without the cross-validation. Default is TRUE. This parameter is ignored when method = 'SMN' or 'slice'.
a	If lambda.estimate = 'MCMC', shape hyperparameter for the Gamma prior on lambda. Default is 0.001. This parameter is ignored when method = 'slice'.
b	If lambda.estimate = 'MCMC', rate hyperparameter for the Gamma prior on lambda. Default is 0.001. This parameter is ignored when method = 'slice'.
posterior.summary.beta	Posterior summary measure for beta (mean, median, or mode). Default is 'mean'.
posterior.summary.lambda	Posterior summary measure for lambda (mean, median, or mode). Default is 'median'.
beta.ci.level	Credible interval level for beta. Default is 0.95 (95%).
lambda.ci.level	Credible interval level for lambda. Default is 0.95 (95%).
seed	Seed value for reproducibility. Default is 1234.
na.rm	Logical. Should missing values (including NaN) be omitted from the calculations?

Value

A list containing the following components is returned:

time	Computational time in minutes.
beta	Posterior estimates of beta.
lowerbeta	Lower limit of the credible interval of beta.
upperbeta	Upper limit of the credible interval of beta.
lambda	Posterior estimate of lambda.
lambdaci	Posterior credible interval of lambda.
beta.post	Post-burn-in posterior samples of beta.
sigma2.post	Post-burn-in posterior samples of sigma2.
lambda.post	Post-burn-in posterior samples of lambda.

Examples

## Not run: 

#########################
# Load Prostate dataset #
#########################

library(ElemStatLearn)
prost<-prostate

###########################################
# Scale data and prepare train/test split #
###########################################

prost.std <- data.frame(cbind(scale(prost[,1:8]),prost$lpsa))
names(prost.std)[9] <- 'lpsa'
data.train <- prost.std[prost$train,]
data.test <- prost.std[!prost$train,]

##################################
# Extract standardized variables #
##################################

y.train   = data.train$lpsa - mean(data.train$lpsa)
y.test <- data.test$lpsa - mean(data.test$lpsa)
x.train = scale(as.matrix(data.train[,1:8], ncol=8))
x.test = scale(as.matrix(data.test[,1:8], ncol=8))

#############################
# Reciprocal Bayesian LASSO #
#############################

fit_BayesRLasso<- BayesRLasso(x.train, y.train)
y.pred.BayesRLasso<-x.test%*%fit_BayesRLasso$beta
mean((y.pred.BayesRLasso - y.test)^2) # Performance on test data

######################################
# Visualization of Posterior Samples #
######################################

##############
# Trace Plot #
##############

library(coda)
plot(mcmc(fit_BayesRLasso$beta.post),density=FALSE,smooth=TRUE)

#############
# Histogram #
#############

library(psych)
multi.hist(fit_BayesRLasso$beta.post,density=TRUE,main="")


## End(Not run)

himelmallick/BayesRecipe documentation built on Aug. 4, 2024, 7:21 a.m.