Bayesian Variable Selection in R

Useful R Packages

• The following R functions are useful in the context of regression and variable selection:
  • MCMCregress from the MCMCpack package
  • bayesreg from the bayesreg package
  • bayeslm from the bayeslm package
  • spikeslab from the spikeslab package for Spike and Slab Prior

Basic Bayesian Regression

MCMCregress function

• This function generates posterior samples for a linear regression model
• Gaussian errors
• Gaussian prior for \(\beta\)
• Inverse-Gamma prior for \(\sigma^2\)
• Gibbs sampler

R Usage

library(MCMCpack)
line   <- list(X = c(-2,-1,0,1,2), Y = c(1,3,3,3,5))  # Using fake data
posterior  <- MCMCregress(Y ~ X, data=line, b0 = 0, B0 = 0.1,
                          c0 = 0.01, d0 = 0.01,
                            burnin = 1000, mcmc = 2000, thin = 2)
# plot(posterior) # plotting the posterior distributions

raftery.diag(posterior) # diagnostics

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95 

You need a sample size of at least 3746 with these values of q, r and s

summary(posterior) # Summary

Iterations = 1001:2999
Thinning interval = 2 
Number of chains = 1 
Sample size per chain = 1000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

              Mean     SD Naive SE Time-series SE
(Intercept) 2.9491 0.4274  0.01352       0.013516
X           0.7881 0.3264  0.01032       0.009737
sigma2      1.1681 1.8806  0.05947       0.067804

2. Quantiles for each variable:

               2.5%    25%    50%    75% 97.5%
(Intercept) 2.08909 2.7161 2.9561 3.1959 3.787
X           0.08524 0.6262 0.7999 0.9594 1.425
sigma2      0.17335 0.3825 0.6655 1.2215 5.285

Spike and Slab Model

The Spike and Slab Model (SSVS)

library(spikeslab)
data("diabetesI")
obj <- spikeslab(Y~.,data=diabetesI,
                 n.iter1 = 1000, n.iter2=2000,
                 intercept = TRUE,
                 fast=TRUE)

## Beta estimates without variable selection
obj$bma.scale

         age          sex          bmi          map           tc          ldl 
   1.0546352 -176.4307456  510.3793679  300.7960221  -22.8506075   -6.8949325 
         hdl          tch          ltg          glu        age.2        bmi.2 
-224.2485177   22.3719867  489.0309185    9.6234135   12.7525865   22.5139545 
       map.2         tc.2        ldl.2        hdl.2        tch.2        ltg.2 
   3.0489192   -0.3355920   -4.4530964    3.0273221    1.3688765   -4.7701312 
       glu.2      age.sex      age.bmi      age.map       age.tc      age.ldl 
  37.8141484  122.7514425    3.1749267   16.9775013    1.9579708   -8.0920417 
     age.hdl      age.tch      age.ltg      age.glu      sex.bmi      sex.map 
   4.3708389    0.5697418   12.8155743   12.0853182   10.8104077   10.4530648 
      sex.tc      sex.ldl      sex.hdl      sex.tch      sex.ltg      sex.glu 
  -0.5107624   -5.3120485    5.6031141   -6.2275763    0.4039139    3.2242668 
     bmi.map       bmi.tc      bmi.ldl      bmi.hdl      bmi.tch      bmi.ltg 
  71.5350155   -1.4557174   -0.1602788   -3.6879989    5.8927310    3.3820885 
     bmi.glu       map.tc      map.ldl      map.hdl      map.tch      map.ltg 
  10.3464233    4.1911610    1.9901894    5.0566182   -2.9279020    4.1697434 
     map.glu       tc.ldl       tc.hdl       tc.tch       tc.ltg       tc.glu 
  -0.3817975    1.1562117    3.6473967  -19.7079914   -4.6475583    5.6606889 
     ldl.hdl      ldl.tch      ldl.ltg      ldl.glu      hdl.tch      hdl.ltg 
  -1.6393168    3.0586068   14.7827413    9.1459831   -3.2686320    6.2649421 
     hdl.glu      tch.ltg      tch.glu      ltg.glu 
   0.6143416   -8.8665670   10.9287598    4.0502587 

Spike and Slab Outputs

## Beta estimates with variable selection
obj$gnet.scale

        age         sex         bmi         map          tc         ldl 
   0.000000 -143.590161  505.724366  271.621915    0.000000    0.000000 
        hdl         tch         ltg         glu       age.2       bmi.2 
-210.177129    0.000000  475.670332   12.313127   13.239734   36.130017 
      map.2        tc.2       ldl.2       hdl.2       tch.2       ltg.2 
   0.000000    0.000000    0.000000    0.000000    0.000000    0.000000 
      glu.2     age.sex     age.bmi     age.map      age.tc     age.ldl 
  71.445803  122.131862    0.000000   25.552380    0.000000    0.000000 
    age.hdl     age.tch     age.ltg     age.glu     sex.bmi     sex.map 
   0.000000    0.000000   11.326147    9.779283    0.000000    5.030626 
     sex.tc     sex.ldl     sex.hdl     sex.tch     sex.ltg     sex.glu 
   0.000000    0.000000    0.000000    0.000000    0.000000    0.000000 
    bmi.map      bmi.tc     bmi.ldl     bmi.hdl     bmi.tch     bmi.ltg 
  93.620681    0.000000    0.000000    0.000000    0.000000    0.000000 
    bmi.glu      map.tc     map.ldl     map.hdl     map.tch     map.ltg 
   0.000000    0.000000    0.000000    0.000000    0.000000    0.000000 
    map.glu      tc.ldl      tc.hdl      tc.tch      tc.ltg      tc.glu 
   0.000000    0.000000    0.000000    0.000000    0.000000    0.000000 
    ldl.hdl     ldl.tch     ldl.ltg     ldl.glu     hdl.tch     hdl.ltg 
   0.000000    0.000000    0.000000    0.000000    0.000000    0.000000 
    hdl.glu     tch.ltg     tch.glu     ltg.glu 
   0.000000    0.000000    0.000000    0.000000 

Spike and Slab Outputs (Contd.)

## Indices of selected variables
obj$gnet.obj.vars

 [1]  3  9  4  7  2 20 37 19  5 12  8 52 22 57 27 11 28 63 29 30 43 10 58 62 24
[26]  6 60 34 41 54 33 32 46 18 53 15 25 44 48 64 40 51 42 59 36 21 56 13 16 47
[51] 45 23 55 38 17 50  1 61 26 31 35 49 14 39

## Summary of Model
obj$summary

                  bma       gnet    bma.scale  gnet.scale
bmi      24.276270994 24.0548551  510.3793679  505.724366
ltg      23.260828671 22.6253304  489.0309185  475.670332
map      14.307407649 12.9197369  300.7960221  271.621915
hdl     -10.666414285 -9.9971066 -224.2485177 -210.177129
sex      -8.391955195 -6.8298878 -176.4307456 -143.590161
age.sex   5.838690696  5.8092203  122.7514425  122.131862
bmi.map   3.402573700  4.4530817   71.5350155   93.620681
glu.2     1.798635637  3.3983303   37.8141484   71.445803
tc       -1.086892561  0.0000000  -22.8506075    0.000000
bmi.2     1.070879610  1.7185297   22.5139545   36.130017
tch       1.064126868  0.0000000   22.3719867    0.000000
tc.tch   -0.937413560  0.0000000  -19.7079914    0.000000
age.map   0.807537386  1.2154028   16.9775013   25.552380
ldl.ltg   0.703143298  0.0000000   14.7827413    0.000000
age.ltg   0.609574709  0.5387299   12.8155743   11.326147
age.2     0.606578675  0.6297499   12.7525865   13.239734
age.glu   0.574839967  0.4651531   12.0853182    9.779283
tch.glu   0.519828094  0.0000000   10.9287598    0.000000
sex.bmi   0.514198658  0.0000000   10.8104077    0.000000
sex.map   0.497201589  0.2392825   10.4530648    5.030626
bmi.glu   0.492129168  0.0000000   10.3464233    0.000000
glu       0.457739102  0.5856757    9.6234135   12.313127
ldl.glu   0.435030053  0.0000000    9.1459831    0.000000
tch.ltg  -0.421739584  0.0000000   -8.8665670    0.000000
age.ldl  -0.384899170  0.0000000   -8.0920417    0.000000
ldl      -0.327958497  0.0000000   -6.8949325    0.000000
hdl.ltg   0.297992907  0.0000000    6.2649421    0.000000
sex.tch  -0.296215598  0.0000000   -6.2275763    0.000000
bmi.tch   0.280288630  0.0000000    5.8927310    0.000000
tc.glu    0.269251514  0.0000000    5.6606889    0.000000
sex.hdl   0.266512959  0.0000000    5.6031141    0.000000
sex.ldl  -0.252668383  0.0000000   -5.3120485    0.000000
map.hdl   0.240518802  0.0000000    5.0566182    0.000000
ltg.2    -0.226892003  0.0000000   -4.7701312    0.000000
tc.ltg   -0.221061806  0.0000000   -4.6475583    0.000000
ldl.2    -0.211812193  0.0000000   -4.4530964    0.000000
age.hdl   0.207899603  0.0000000    4.3708389    0.000000
map.tc    0.199353199  0.0000000    4.1911610    0.000000
map.ltg   0.198334467  0.0000000    4.1697434    0.000000
ltg.glu   0.192651162  0.0000000    4.0502587    0.000000
bmi.hdl  -0.175420221  0.0000000   -3.6879989    0.000000
tc.hdl    0.173488970  0.0000000    3.6473967    0.000000
bmi.ltg   0.160869548  0.0000000    3.3820885    0.000000
hdl.tch  -0.155472969  0.0000000   -3.2686320    0.000000
sex.glu   0.153362733  0.0000000    3.2242668    0.000000
age.bmi   0.151015860  0.0000000    3.1749267    0.000000
ldl.tch   0.145483089  0.0000000    3.0586068    0.000000
map.2     0.145022298  0.0000000    3.0489192    0.000000
hdl.2     0.143995026  0.0000000    3.0273221    0.000000
map.tch  -0.139266098  0.0000000   -2.9279020    0.000000
map.ldl   0.094663658  0.0000000    1.9901894    0.000000
age.tc    0.093131173  0.0000000    1.9579708    0.000000
ldl.hdl  -0.077974348  0.0000000   -1.6393168    0.000000
bmi.tc   -0.069241415  0.0000000   -1.4557174    0.000000
tch.2     0.065110815  0.0000000    1.3688765    0.000000
tc.ldl    0.054995383  0.0000000    1.1562117    0.000000
age       0.050163881  0.0000000    1.0546352    0.000000
hdl.glu   0.029221250  0.0000000    0.6143416    0.000000
age.tch   0.027099855  0.0000000    0.5697418    0.000000
sex.tc   -0.024294491  0.0000000   -0.5107624    0.000000
sex.ltg   0.019212225  0.0000000    0.4039139    0.000000
map.glu  -0.018160257  0.0000000   -0.3817975    0.000000
tc.2     -0.015962484  0.0000000   -0.3355920    0.000000
bmi.ldl  -0.007623687  0.0000000   -0.1602788    0.000000

High Dimensional Usage

# add 2000 noise variables
data(diabetesI, package = "spikeslab")
diabetes.noise <- cbind(diabetesI,
      noise = matrix(rnorm(nrow(diabetesI) * 2000), nrow(diabetesI)))

# example of a big p, small n call
# don't use formula call; make call with x and y arguments
x <- diabetes.noise[, -1]
y <- diabetes.noise[, 1]
obj <- spikeslab(x=x, y=y, bigp.smalln=TRUE, max.var=100,
                 n.iter1=1000, n.iter2=2000)

obj$gnet.obj.vars

  [1]    3    9    7 1544  457  627    8  906  839  693 1177 1114  537  640 1772
 [16] 2008 1522   97 1969  675 1598 1179  892 1661  746  183 1120  289  865  306
 [31]  901 2063  121 1418  151  979 1485  276  907  541  546 1117 1734 1005  851
 [46] 1737  663  506  781  585   74 1859 1583  637 1034  496 1234 1803  980 1634
 [61] 1862 1095  587  334  764 1595 1856 1771  520 1918 1466 1755 1685 1294  735
 [76] 1789 1267  791  417 1490  339  766 1482  432 1310 1993 1891 1770  368 1746
 [91] 1061 1735 1205 1930  402  834 1104 1861 1204 1221

Continuous Shrinkage Priors

R Packages

• Both bayesreg and bayeslm does variable selection using continuous shrinkage priors

• bayesreg allows for Gaussian, Laplace, Student-t errors, as well as GLM for count and binary data

• It allows for ridge regression, LASSO, Horseshoe, and Horseshoe+ priors

• bayeslm only allows Gaussian errors and ridge regression, LASSO (Laplace), Horseshoe and some other priors using elliptical slice sampler

• You can run a non-penalized model with bayeslm as well

bayeslm Usage

n = 100
kappa = 1.25
beta_true = c(c(1,2,3),rnorm(17,0,0.01))
sig_true = kappa*sqrt(sum(beta_true^2))

x = matrix(rnorm(20*n),n,20)
y = x %*% beta_true + sig_true * rnorm(n)
df = data.frame(x = x, y = y)
library(bayeslm)
fita = bayeslm(y ~ x, prior = 'horseshoe', data = df,
        block_vec = rep(1,20), N = 3000, burnin = 1000)

bayeslm Outputs

summary(fita)

Average number of rejections before one acceptance : 
29.72667 
Summary of beta draws 
   based on 2400 valid draws (number of burn in is 600) 
Summary of Posterior draws 
 Moments 
               mean std dev eff sample size    
(Intercept)  1.8082    0.46            2146 ***
x1           0.5466    0.49             247    
x2           2.0771    0.48             471 ***
x3           2.2983    0.50             670 ***
x4           0.0022    0.25            1878    
x5          -0.0592    0.27            1460    
x6           0.1694    0.32             510    
x7          -0.0477    0.25            1398    
x8          -0.3635    0.43             335    
x9          -0.7844    0.53             184   .
x10          0.2089    0.34             484    
x11          0.0710    0.25            1054    
x12          0.2214    0.36             624    
x13         -0.0346    0.25            1937    
x14         -0.1034    0.26             918    
x15          0.0133    0.26            1340    
x16          0.1911    0.31             509    
x17          0.0713    0.26            1164    
x18         -0.0237    0.24            1522    
x19         -0.1088    0.26             927    
x20          0.1328    0.29             683    
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Summary of sigma draws 
Mean of standard deviation is  4.401584 
S.d. of standard deviation samples is  0.3434586 
Effective sample size of s.d. is  913.9362 

bayesreg Usage

X = matrix(rnorm(100*20),100,20)
b = matrix(0,20,1)
b[1:5] = c(5,4,3,2,1)
y = X %*% b + rnorm(100, 0, 1)

df <- data.frame(X,y)

library(bayesreg)
rv.hs <- bayesreg(y~., data=df, model="gaussian", prior="hs",
                  n.samples=2000, burnin=1000, thin=1)