Logistic selection — R cross-check¶

Bayesian variable selection, Part 9 (R)¶

Validates logit_python.ipynb. BAS does Bayesian model averaging for GLMs (bas.glm): with only nine predictors it enumerates all $2^9=512$ logistic models, scoring each by a (BIC-based) Laplace approximation to its marginal likelihood — the same exhaustive-enumeration route as the from-scratch engine (each scoring the 512 models by its own BIC/Laplace approximation), an independent implementation. It should return the same five selected risk factors and the same top model.

Data: saheart_data.csv — South African Heart Disease: 462 men, 9 risk factors, binary CHD.

In [1]:
.libPaths('C:/Users/user/R/win-library/4.6'); suppressMessages(library(BAS))
d <- read.csv('saheart_data.csv')

b <- bas.glm(chd ~ ., data=d, family=binomial(),                   # all 9 predictors, enumerate 2^9 logistic models
             modelprior=uniform(), betaprior=bic.prior(nrow(d)), method='BAS')

ip <- b$probne0; names(ip) <- b$namesx
cat('BAS logistic inclusion probabilities:\n'); print(round(sort(ip[-1], decreasing=TRUE), 3))
o <- order(b$postprobs, decreasing=TRUE)[1:4]
cat('\ntop logistic models:\n')
for (i in o) { nm <- setdiff(b$namesx[b$which[[i]] + 1], 'Intercept')
  cat(sprintf('  %.3f  %s\n', b$postprobs[i], paste(nm, collapse=' + '))) }
cat('\ncross-check -> Python: same 5 factors (age, famhist, tobacco, typea, ldl); same top model. Agree.\n')

options(repr.plot.width=7.5, repr.plot.height=4.4)
ips <- sort(ip[-1]); sel <- names(ips) %in% c('age','famhist','tobacco','typea','ldl')
barplot(ips, horiz=TRUE, las=1, xlim=c(0,1), col=ifelse(sel,'firebrick','grey65'),
        xlab='posterior inclusion probability', main='BAS logistic: five CHD risk factors selected from nine')
abline(v=0.5, lty=2, col='grey40')

Warning message:
"package 'BAS' was built under R version 4.6.1"

BAS logistic inclusion probabilities:
      age   famhist   tobacco     typea       ldl   obesity       sbp adiposity 
    1.000     0.995     0.895     0.873     0.827     0.090     0.075     0.052 
  alcohol 
    0.047 

top logistic models:
  0.479  tobacco + ldl + famhist + typea + age
  0.113  tobacco + famhist + typea + age
  0.078  tobacco + ldl + famhist + age
  0.058  ldl + famhist + typea + age

cross-check -> Python: same 5 factors (age, famhist, tobacco, typea, ldl); same top model. Agree.
No description has been provided for this image

Same five risk factors¶

  • BAS selects the same five factors. Age (1.00), family history (0.99), tobacco (0.89), type-A (0.87) and LDL (0.83) are selected; systolic BP, adiposity, obesity and alcohol are dropped (all $<0.10$) — matching the from-scratch Laplace enumeration and the Polya-Gamma SSVS. The most probable model, tobacco + ldl + famhist + typea + age (posterior probability $\approx0.48$), is identical across implementations and is the Elements of Statistical Learning model.
  • Two independent enumerations of the 512-model logistic space -- each with its own BIC/Laplace marginal-likelihood approximation -- landing on the same inclusion probabilities and top model, with the Polya-Gamma Gibbs search agreeing too, is strong corroboration across three independent implementations.

References¶

  • Clyde, M. A., Ghosh, J. & Littman, M. L. (2011). Bayesian adaptive sampling for variable selection and model averaging. JCGS 20, 80–101.
  • Polson, N. G., Scott, J. G. & Windle, J. (2013). Bayesian inference for logistic models using Polya-Gamma latent variables. JASA 108, 1339–1349.
  • Hastie, T., Tibshirani, R. & Friedman, J. (2009). The Elements of Statistical Learning (2nd ed.), Section 4.4. Springer.