Exact enumeration & BMA — R cross-check

Bayesian variable selection, Part 2 (R)

Validates bvs_python.ipynb. The from-scratch enumerator scored all $2^9=512$ models under a $g$-prior with $g=n$. BAS computes exactly the same marginal likelihood, so with prior='g-prior', alpha=n it must reproduce the from-scratch posterior inclusion probabilities and top models to the decimal — the tightest kind of cross-check, two independent implementations of the same exact calculation. It also gives the classic image() picture of the model space.

Data: cystfibr_data.csv — cystic-fibrosis PEmax data, 25 patients, 9 predictors (ISwR::cystfibr; Congdon Example 3.3).

.libPaths('C:/Users/user/R/win-library/4.6'); suppressMessages(library(BAS))
d <- read.csv('cystfibr_data.csv')
b <- bas.lm(pemax ~ age + sex + height + weight + bmp + fev1 + rv + frc + tlc,   # all 9 predictors
            data=d, prior='g-prior', alpha=nrow(d), modelprior=uniform())        # g = n = 25; enumerate 2^9

ip <- b$probne0; names(ip) <- b$namesx
cat('marginal inclusion probabilities:\n'); print(round(sort(ip[-1], decreasing=TRUE), 3))
o <- order(b$postprobs, decreasing=TRUE)[1:5]
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 enumeration: weight 0.641, fev1 0.596, bmp 0.492; best model weight+bmp+fev1 (0.051). Exact match.\n')

options(repr.plot.width=7.5, repr.plot.height=4.2)
ips <- sort(ip[-1])
cols <- ifelse(names(ips) %in% c('weight','fev1','bmp'), 'firebrick', 'grey60')
barplot(ips, horiz=TRUE, las=1, xlim=c(0,1), col=cols, xlab='posterior inclusion probability',
        main='BAS exact enumeration (g-prior): weight, FEV1, BMP lead')
abline(v=0.5, lty=2, col='grey40')

marginal inclusion probabilities:
weight   fev1    bmp    age     rv height    tlc    frc    sex 
 0.641  0.596  0.492  0.376  0.303  0.296  0.235  0.225  0.213 

top models by posterior probability:
  0.051  weight + bmp + fev1
  0.035  weight
  0.033  weight + bmp + fev1 + rv
  0.029  weight + bmp
  0.022  weight + bmp + fev1 + frc

cross-check -> Python enumeration: weight 0.641, fev1 0.596, bmp 0.492; best model weight+bmp+fev1 (0.051). Exact match.

Identical to the decimal

• BAS reproduces the from-scratch enumeration exactly. Same inclusion probabilities (weight 0.641, FEV1 0.596, BMP 0.492, …) and the same ranked models (best weight + bmp + fev1 at 0.051). Because both evaluate the same closed-form $g$-prior marginal likelihood over the same 512 models, agreement is exact rather than approximate — the strongest validation available.
• Same story, too: three influential predictors, but a diffuse posterior over models — the evidence for any single model is weak, so model averaging is the honest summary.

References

• Clyde, M. A., Ghosh, J. & Littman, M. L. (2011). Bayesian adaptive sampling for variable selection and model averaging. JCGS 20, 80–101.
• Liang, F. et al. (2008). Mixtures of g-priors for Bayesian variable selection. JASA 103, 410–423.