Illustrative Example: Why Bayesian Variable Selection Doesn't Scale

code.R

## ##########################################################################################################################
## ##########################################################################################################################
## @section: Prep workspace
## ##########################################################################################################################
## ##########################################################################################################################

options(width=200, digits=6, digits.secs=6)
rm(list=ls())




## ##########################################################################################################################
## ##########################################################################################################################
## @section: load data
## ##########################################################################################################################
## ##########################################################################################################################

mat.df.DATA <- read.csv("data.csv")
scl.int.N   <- dim(mat.df.DATA)[1]
scl.flt.LAM <- log(scl.int.N)/scl.int.N

vec.flt.R  <- mat.df.DATA$r
vec.flt.X1 <- mat.df.DATA$x1
vec.flt.X2 <- mat.df.DATA$x2
vec.flt.X3 <- mat.df.DATA$x3




## ##########################################################################################################################
## ##########################################################################################################################
## @section: evaluate models
## ##########################################################################################################################
## ##########################################################################################################################

mat.df.RESULTS <- data.frame(activeSet = c("0", "1", "2", "3", "12", "13", "23", "123"),
                             error     = NA,
                             penalty   = NA
                             )

## @active: none
mat.df.RESULTS[mat.df.RESULTS$activeSet == "0", ]$error   <- var(lm(vec.flt.R ~ 0)$residuals)
mat.df.RESULTS[mat.df.RESULTS$activeSet == "0", ]$penalty <- scl.flt.LAM * 0

## @active: x1
mat.df.RESULTS[mat.df.RESULTS$activeSet == "1", ]$error   <- var(lm(vec.flt.R ~ 0 + vec.flt.X1)$residuals)
mat.df.RESULTS[mat.df.RESULTS$activeSet == "1", ]$penalty <- scl.flt.LAM * 1

## @active: x2
mat.df.RESULTS[mat.df.RESULTS$activeSet == "2", ]$error   <- var(lm(vec.flt.R ~ 0 + vec.flt.X2)$residuals)
mat.df.RESULTS[mat.df.RESULTS$activeSet == "2", ]$penalty <- scl.flt.LAM * 1

## @active: x3
mat.df.RESULTS[mat.df.RESULTS$activeSet == "3", ]$error   <- var(lm(vec.flt.R ~ 0 + vec.flt.X3)$residuals)
mat.df.RESULTS[mat.df.RESULTS$activeSet == "3", ]$penalty <- scl.flt.LAM * 1

## @active: x1, x2
mat.df.RESULTS[mat.df.RESULTS$activeSet == "12", ]$error   <- var(lm(vec.flt.R ~ 0 + vec.flt.X1 + vec.flt.X2)$residuals)
mat.df.RESULTS[mat.df.RESULTS$activeSet == "12", ]$penalty <- scl.flt.LAM * 2

## @active: x1, x3
mat.df.RESULTS[mat.df.RESULTS$activeSet == "13", ]$error   <- var(lm(vec.flt.R ~ 0 + vec.flt.X1 + vec.flt.X3)$residuals)
mat.df.RESULTS[mat.df.RESULTS$activeSet == "13", ]$penalty <- scl.flt.LAM * 2

## @active: x2, x3
mat.df.RESULTS[mat.df.RESULTS$activeSet == "23", ]$error   <- var(lm(vec.flt.R ~ 0 + vec.flt.X2 + vec.flt.X3)$residuals)
mat.df.RESULTS[mat.df.RESULTS$activeSet == "23", ]$penalty <- scl.flt.LAM * 2

## @active: x1, x2, x3
mat.df.RESULTS[mat.df.RESULTS$activeSet == "123", ]$error   <- var(lm(vec.flt.R ~ 0 + vec.flt.X1 + vec.flt.X2 + vec.flt.X3)$residuals)
mat.df.RESULTS[mat.df.RESULTS$activeSet == "123", ]$penalty <- scl.flt.LAM * 3

mat.df.RESULTS$objective <- mat.df.RESULTS$error + mat.df.RESULTS$penalty




## ##########################################################################################################################
## ##########################################################################################################################
## @section: report results
## ##########################################################################################################################
## ##########################################################################################################################

print(round(scl.flt.LAM, 3))
print(round(c(mean(vec.flt.X1), mean(vec.flt.X2), mean(vec.flt.X3)), 3))
print(round(var(cbind(vec.flt.X1, vec.flt.X2, vec.flt.X3)), 3))
print(round(c(cor(vec.flt.R, vec.flt.X1), cor(vec.flt.R, vec.flt.X2), cor(vec.flt.R, vec.flt.X3)), 3))

print(mat.df.RESULTS)

data.csv