SSVS / Bayesian variable selection — R cross-check¶

Bayesian variable selection, Part 1 (R)¶

Validates ssvs_python.ipynb. The Python notebook searched the model space by MCMC (spike-and-slab Gibbs, and PyMC). Here we take the complementary exact route: with only $p=6$ predictors there are just $2^6=64$ models, so BAS (Bayesian Adaptive Sampling, Clyde et al. 2011) can enumerate every one, compute its marginal likelihood in closed form, and return exact posterior model probabilities and marginal inclusion probabilities. If the MCMC search is working, its inclusion probabilities should line up with this exact enumeration — and the most-probable model should again be "soil-moisture tension only".

Data: ssvs_data.csv — the turnip-greens vitamin-B2 experiment of Draper & Smith (1966), *Applied Regression Analysis (Congdon Example 3.10): riboflavin content of turnip greens vs radiation, soil-moisture tension and temperature over 27 plots.*

In [2]:
.libPaths('C:/Users/user/R/win-library/4.6'); suppressMessages(library(BAS))
d <- read.csv('ssvs_data.csv')
names(d) <- c('radiation','moisture','temp','riboflavin')          # X1,X2,X3,Y -> real names (Draper-Smith)
d$rad_moist <- d$radiation*d$moisture                              # the 3 pairwise interactions
d$rad_temp  <- d$radiation*d$temp
d$moist_temp<- d$moisture*d$temp
f <- riboflavin ~ radiation + moisture + temp + rad_moist + rad_temp + moist_temp

# enumerate ALL 2^6 = 64 models; g-prior on coefficients, uniform prior over models
b <- bas.lm(f, data=d, prior='g-prior', alpha=nrow(d), modelprior=uniform())

ip <- b$probne0; names(ip) <- b$namesx
cat('marginal inclusion probabilities  P(gamma_j = 1 | data):\n'); print(round(ip[-1], 3))

o <- order(b$postprobs, decreasing=TRUE)[1:5]                       # most probable models
cat('\ntop models by posterior probability:\n')
for (i in o) { nm <- setdiff(b$namesx[b$which[[i]] + 1], 'Intercept')
  cat(sprintf('  %.3f  %s\n', b$postprobs[i], if(length(nm)) paste(nm, collapse=' + ') else '(intercept only)')) }
cat('\ncross-check -> Python SSVS: same selection (moisture wins, MAP = moisture only); inclusion 0.61 here vs ~0.69 there differs by prior, not MC error.\n')

options(repr.plot.width=7, repr.plot.height=4)
cols <- ifelse(names(ip)=='moisture', 'firebrick', 'grey60')
bp <- barplot(ip[-1], col=cols[-1], ylim=c(0,1), ylab='posterior inclusion probability', las=2,
              main='BAS exact enumeration: soil-moisture tension is the predictor that matters')
abline(h=0.5, lty=2, col='grey40'); text(bp[2], ip['moisture']+0.05, sprintf('%.2f', ip['moisture']), col='firebrick')

marginal inclusion probabilities  P(gamma_j = 1 | data):
 radiation   moisture       temp  rad_moist   rad_temp moist_temp 
     0.231      0.614      0.232      0.201      0.264      0.449 

top models by posterior probability:
  0.185  moisture
  0.105  moist_temp
  0.073  moisture + rad_temp
  0.062  radiation + moisture
  0.060  moisture + temp

cross-check -> Python SSVS: same selection (moisture wins, MAP = moisture only); inclusion 0.61 here vs ~0.69 there differs by prior, not MC error.
No description has been provided for this image

Same selection, by a different route¶

  • Exact enumeration confirms the search. Across all 64 models, soil-moisture tension carries the largest marginal inclusion probability (~0.61) and the single most-probable model is "moisture only" (~0.19) — the same qualitative selection as the from-scratch spike-and-slab Gibbs and PyMC (moisture wins, "moisture only" is the MAP model, moisture:temperature is the runner-up). The exact inclusion probabilities differ modestly — 0.61 here vs ~0.69 there — because the two routes use different priors (a g-prior with exact enumeration vs a spike-and-slab MCMC), not because of Monte-Carlo error.
  • Two philosophies, one answer. SSVS explores model space stochastically (essential when $p$ is large and $2^p$ is astronomical); BAS here evaluates it exhaustively (possible only because $p$ is small). That they coincide is the validation: the MCMC search is faithfully sampling the true posterior over models.

References¶

  • Clyde, M. A., Ghosh, J. & Littman, M. L. (2011). Bayesian adaptive sampling for variable selection and model averaging. J. Computational and Graphical Statistics 20, 80–101.
  • George, E. I. & McCulloch, R. E. (1993). Variable selection via Gibbs sampling. JASA 88, 881–889.