BayesianLasso

Introduction

The ‘BayesianLasso’ package provides efficient sampling algorithms for Bayesian Lasso regression, modified Hans and PC samplers. The modified Hans sampler is based on a newly defined Lasso distribution.

This vignette introduces the main functionality of the package and demonstrates how to apply the samplers to a simulated data.

Simulated Example

library(BayesianLasso)
# Simulate data
set.seed(123)
Ns <- 2000
ns <- 100
ps <- 10
X <- matrix(rnorm(ns * ps), nrow = ns)
beta <- c(rep(2, 3), rep(0, ps - 3))
y <- X %*% beta + rnorm(ns)

vtime_val_Hans = c()
results_Hans <- NULL 

# Run the modified Hans sampler
for(i in 1:5){
  time_val <- system.time({
    res_Hans <- Modified_Hans_Gibbs(
      X = X, y = y, beta_init= rep(1,10), a1=2, b1=1, u1=2, v1=1,
      nsamples=Ns, lambda_init=1, sigma2_init=1, verbose=0, tune_lambda2 = TRUE,
      rao_blackwellization = FALSE)
  })[3]
  
  vtime_val_Hans[i] <- time_val
  
  # Initialize accumulators after first run
  if (is.null(results_Hans)) {
    results_Hans <- list(
      mBeta_Hans = res_Hans$mBeta,
      vsigma2_Hans = res_Hans$vsigma2,
      vlambda2_Hans = res_Hans$vlambda2
    )
  } else {
    results_Hans$mBeta_Hans    <- results_Hans$mBeta_Hans + res_Hans$mBeta
    results_Hans$vsigma2_Hans  <- results_Hans$vsigma2_Hans + res_Hans$vsigma2
    results_Hans$vlambda2_Hans <- results_Hans$vlambda2_Hans + res_Hans$vlambda2
  }
}

# Take averages
mBeta_Hans    <- (results_Hans$mBeta_Hans) / 5
vsigma2_Hans  <- (results_Hans$vsigma2_Hans) / 5
vlambda2_Hans <- (results_Hans$vlambda2_Hans) / 5
time_val_Hans = mean(vtime_val_Hans)




# ======================== PC sampler =======================================

# Run the modified PC sampler

vtime_val_PC = c()
results_PC <- NULL 

for(i in 1:5){
  time_val <- system.time({
    res_PC <- Modified_PC_Gibbs(
      X = X, y = y, a1=2, b1=1, u1=2, v1=1,
      nsamples=Ns, lambda_init=1, sigma2_init=1, verbose=0)
  })[3]
  
  vtime_val_PC[i] <- time_val
  
  # Initialize accumulators after first run
  if (is.null(results_PC)) {
    results_PC <- list(
      mBeta = res_PC$mBeta,
      vsigma2 = res_PC$vsigma2,
      vlambda2 = res_PC$vlambda2
    )
  } else {
    results_PC$mBeta    <- results_PC$mBeta + res_PC$mBeta
    results_PC$vsigma2  <- results_PC$vsigma2 + res_PC$vsigma2
    results_PC$vlambda2 <- results_PC$vlambda2 + res_PC$vlambda2
  }
}

# Take averages
mBeta = (results_PC$mBeta)/5
vsigma2 = (results_PC$vsigma2)/5
vlambda2 = (results_PC$vlambda2)/5
time_val_PC = mean(vtime_val_PC)

Convergence Diagnostics

Assessing convergence of the MCMC chains is an important step in Bayesian analysis. We use the mcmc_stats() and mcmc_diagnostics() helper functions to evaluate the mixing and convergence of the Modified Hans Gibbs sampler.

stats_Hans <- mcmc_stats(
  mBeta_Hans, vsigma2_Hans, vlambda2_Hans,
  time_val = time_val_Hans, inds_use = 200:Ns
)

print(stats_Hans)
#>    mix_beta    eff_beta  mix_sigma2  eff_sigma2 mix_lambda2 eff_lambda2 
#>    82.68013 62044.54438    87.51104 65669.73972    72.05850 54073.90298 
#>        time 
#>     0.02400

mcmc_diagnostics(
  mBeta_Hans, vsigma2_Hans, vlambda2_Hans,
  beta_inds = 1:3, mStat = stats_Hans, doplots = TRUE
)

#> $ldens_beta
#> $ldens_beta[[1]]
#> 
#> Call:
#>  density.default(x = mBeta[, j])
#> 
#> Data: mBeta[, j] (2000 obs.);    Bandwidth 'bw' = 0.01009
#> 
#>        x               y           
#>  Min.   :1.838   Min.   :0.000222  
#>  1st Qu.:1.949   1st Qu.:0.082179  
#>  Median :2.061   Median :0.833078  
#>  Mean   :2.061   Mean   :2.229172  
#>  3rd Qu.:2.173   3rd Qu.:4.316045  
#>  Max.   :2.285   Max.   :7.579985  
#> 
#> $ldens_beta[[2]]
#> 
#> Call:
#>  density.default(x = mBeta[, j])
#> 
#> Data: mBeta[, j] (2000 obs.);    Bandwidth 'bw' = 0.009996
#> 
#>        x               y           
#>  Min.   :1.757   Min.   :0.000233  
#>  1st Qu.:1.860   1st Qu.:0.103123  
#>  Median :1.963   Median :0.877695  
#>  Mean   :1.963   Mean   :2.422321  
#>  3rd Qu.:2.066   3rd Qu.:4.713308  
#>  Max.   :2.169   Max.   :7.850303  
#> 
#> $ldens_beta[[3]]
#> 
#> Call:
#>  density.default(x = mBeta[, j])
#> 
#> Data: mBeta[, j] (2000 obs.);    Bandwidth 'bw' = 0.009598
#> 
#>        x               y          
#>  Min.   :1.670   Min.   :0.00000  
#>  1st Qu.:1.797   1st Qu.:0.01344  
#>  Median :1.924   Median :0.24579  
#>  Mean   :1.924   Mean   :1.96578  
#>  3rd Qu.:2.051   3rd Qu.:3.55550  
#>  Max.   :2.178   Max.   :8.09667  
#> 
#> 
#> $dens_sigma2
#> 
#> Call:
#>  density.default(x = vsigma2)
#> 
#> Data: vsigma2 (2000 obs.);   Bandwidth 'bw' = 0.01358
#> 
#>        x                y           
#>  Min.   :0.8083   Min.   :0.000165  
#>  1st Qu.:0.9552   1st Qu.:0.042315  
#>  Median :1.1021   Median :0.626729  
#>  Mean   :1.1021   Mean   :1.698430  
#>  3rd Qu.:1.2490   3rd Qu.:3.318473  
#>  Max.   :1.3959   Max.   :5.842377  
#> 
#> $dens_lambda2
#> 
#> Call:
#>  density.default(x = vlambda2)
#> 
#> Data: vlambda2 (2000 obs.);  Bandwidth 'bw' = 0.09214
#> 
#>        x                y            
#>  Min.   :0.7931   Min.   :0.0000243  
#>  1st Qu.:1.6775   1st Qu.:0.0163150  
#>  Median :2.5619   Median :0.1569984  
#>  Mean   :2.5619   Mean   :0.2821238  
#>  3rd Qu.:3.4463   3rd Qu.:0.5399681  
#>  Max.   :4.3308   Max.   :0.8281255  
#> 
#> $mStat
#>    mix_beta    eff_beta  mix_sigma2  eff_sigma2 mix_lambda2 eff_lambda2 
#>    82.68013 62044.54438    87.51104 65669.73972    72.05850 54073.90298 
#>        time 
#>     0.02400 
#> 
#> $beta_inds
#> [1] 1 2 3
#> 
#> $rhat_sigma2
#> [1] 1.000383
#> 
#> $rhat_lambda2
#> [1] 0.9997187