R/BayesRLassoSlice.R

In himelmallick/BayesRecipe: Bayesian Reciprocal Regularization

#' Bayesian Reciprocal Regularization
#'
#' Elliptical slice sampler for Bayesian Reciprocal LASSO using inverse Laplace prior 
#' 
#' @param x A numeric matrix with standardized predictors in columns and samples in rows.
#' @param y A mean-centered continuous response variable with matching rows with x.
#' @param lambda.estimate Estimating lambda by empirical bayes ('EB'), MCMC ('MCMC'), or apriori ('AP'). Default is 'AP'.
#' @param update.sigma2 Whether sigma2 should be updated. Default is TRUE.
#' @param max.steps Number of MCMC iterations. Default is 11000.
#' @param n.burn Number of burn-in iterations. Default is 1000.
#' @param n.thin Lag at which thinning should be done. Default is 1 (no thinning).
#' @param ridge.CV If 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.
#' @param a If lambda.estimate = 'MCMC', shape hyperparameter for the Gamma prior on lambda. Default is 0.001.
#' @param b If lambda.estimate = 'MCMC', rate hyperparameter for the Gamma prior on lambda. Default is 0.001.
#' @param posterior.summary.beta Posterior summary measure for beta (mean, median, or mode). Default is 'mean'. 
#' @param posterior.summary.lambda Posterior summary measure for lambda (mean, median, or mode). Default is 'median'. 
#' @param beta.ci.level Credible interval level for beta. Default is 0.95 (95\%).
#' @param lambda.ci.level Credible interval level for lambda. Default is 0.95 (95\%).
#' @param seed Seed value for reproducibility. Default is 1234.
#' @param na.rm Logical. Should missing values (including NaN) be omitted from the calculations?
#' 
#' @return A list containing the following components is returned:
#' \item{time}{Computational time in minutes.}
#' \item{beta}{Posterior estimates of beta.}
#' \item{lowerbeta}{Lower limit of the credible interval of beta.}
#' \item{upperbeta}{Upper limit of the credible interval of beta.}
#' \item{lambda}{Posterior estimate of lambda.}
#' \item{lambdaci}{Posterior credible interval of lambda.}
#' \item{beta.post}{Post-burn-in posterior samples of beta.}
#' \item{sigma2.post}{Post-burn-in posterior samples of sigma2.}
#' \item{lambda.post}{Post-burn-in posterior samples of lambda.}
#' 
#' @examples
#' \dontrun{
#' 
#' #########################
#' # 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 (slice) #
#' ###################################
#' 
#' fit_BayesRLasso_slice<- BayesRLasso(x.train, y.train, method = 'slice')
#' y.pred.BayesRLasso_slice<-x.test%*%fit_BayesRLasso_slice$beta
#' mean((y.pred.BayesRLasso_slice - y.test)^2) # Performance on test data
#'
#' ######################################
#' # Visualization of Posterior Samples #
#' ######################################
#' 
#' ##############
#' # Trace Plot #
#' ##############
#' 
#' library(coda)
#' plot(mcmc(fit_BayesRLasso_slice$beta.post),density=FALSE,smooth=TRUE)
#' 
#' #############
#' # Histogram #
#' #############
#' 
#' library(psych)
#' multi.hist(fit_BayesRLasso_slice$beta.post,density=TRUE,main="")
#' 
#' }
#' 
#' @export
#' @keywords Ellipitical slice sampler, MCMC, Bayesian regularization, Reciprocal LASSO

BayesRLassoSlice<-function(x,
                         y, 
                         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) {
  # Set random seed
  set.seed(seed)
  
  # Print parameters ignored by the slice Sampler
  cat("The following parameters will be ignored by the slice sampler:\n")
      cat("update.sigma2, ridge.CV, a, b\n")
  
  # Set parameter values and extract relevant quantities
	x <- as.matrix(x)
 	n <- nrow(x)
	p <- ncol(x)
	# Calculate cross-product quantities for OLS
	XtX <- t(x) %*% x	

	# Estimate lambda (remains the same throughout the sampler if estimated beforehand)
	# EB and MCMC not supported for the slice sampler)
	if (lambda.estimate == 'AP') {
	  cat("Finding the best lambda using the apriori method. The slice sampler will start soon...\n")
	  lambda<-hyper_par_BayesRLasso(x, y)
	} else{
	  stop("The selected lambda.estimate is not valid for the slice sampler. Use lambda.estimate = 'AP'")
	}
	lambdaSamples<-rep(lambda, max.steps)
	
	# Track time
	cat(c("Job started at:",date()),fill=TRUE)
	start.time <- Sys.time()

	#####################
	# Run slice Sampler #
	#####################
	
	# Check if X is rank-deficient and correct option accordingly
	v <- try(solve(XtX), silent = TRUE)
	if(inherits(v, "try-error")){
	  fit = bayeslm::bayeslm(y, 
	                         x, 
	                         prior = 'inverselaplace', 
	                         lambda = lambda,  
	                         icept = FALSE, 
	                         N = max.steps,
	                         burnin = 0,
	                         standardize = FALSE, 
	                         singular = TRUE)
	} else{
	  fit = bayeslm::bayeslm(y, 
	                         x, 
	                         prior = 'inverselaplace', 
	                         lambda = lambda,  
	                         icept = FALSE, 
	                         N = max.steps,
	                         burnin = 0,
	                         standardize = FALSE, 
	                         singular = FALSE)
	}
	

	# Collect all quantities and prepare output
	beta.post=fit$beta[seq(n.burn+1, max.steps, n.thin),]
	sigma2.post=fit$sigma[seq(n.burn+1, max.steps, n.thin)]
	lambda.post=lambdaSamples[seq(n.burn+1, max.steps, n.thin)]

	# Posterior estimates of beta
	if (posterior.summary.beta == 'mean') {
	  beta <-apply(beta.post, 2, mean, na.rm = na.rm)
	} else if (posterior.summary.beta == 'median') {
	  beta <-apply(beta.post, 2, median, na.rm = na.rm)
	} else if (posterior.summary.beta == 'mode') {
	  beta <-posterior.mode(mcmc(beta.post), na.rm = na.rm)
	} else{
	  stop('posterior.summary.beta must be either mean, median, or mode.')
	}
	
	# Posterior estimates of lambda
	if (posterior.summary.lambda == 'mean') {
	  lambda<-mean(lambda.post, na.rm = na.rm)
	} else if (posterior.summary.lambda == 'median') {
	  lambda<-median(lambda.post, na.rm = na.rm)
	} else if (posterior.summary.lambda == 'mode') {
	  lambda <-posterior.mode(mcmc(lambda.post), na.rm = na.rm)
	} else{
	  stop('posterior.summary.lambda must be either mean, median, or mode.')
	}
	
	# Credible intervals for beta
	lowerbeta<-colQuantiles(beta.post,prob=(1 - beta.ci.level)/2, na.rm = na.rm)
	upperbeta<-colQuantiles(beta.post,prob=(1 + beta.ci.level)/2, na.rm = na.rm)
	
	# Assign names
	names(beta)<-names(lowerbeta)<-names(upperbeta)<-colnames(beta.post)<-colnames(x)
	
	# Credible intervals for lambda
	lambdaci<-credInt(lambda.post, cdf = NULL,conf=lambda.ci.level, type="twosided")
	
	# Track stopping time
	stop.time <- Sys.time()
	time<-round(difftime(stop.time, start.time, units="min"),3)
	cat(c("Job finished at:",date()),fill=TRUE)
	cat("Computational time:", time, "minutes \n")
	
	# Return output
	return(list(time = format(time),
	            beta = beta,
	            lowerbeta = lowerbeta,
	            upperbeta = upperbeta,
	            lambda = lambda,
	            lambdaci = lambdaci,
	            beta.post = beta.post,
	            sigma2.post = sigma2.post,
	            lambda.post = lambda.post))
}