High-dimensional selection — R cross-check

Bayesian variable selection, Part 8 (R)

Validates hshd_python.ipynb. The standard frequentist tool for $p\gg n$ selection is the lasso (glmnet): an $L_1$ penalty that sets most coefficients exactly to zero, with the penalty strength chosen by cross-validation. It is the natural comparison for the Bayesian horseshoe — both deliver sparse selection, one by penalizing, the other by shrinking. It should put the same gene, YOAB_at, at the top.

Data: ribo_data.csv — riboflavin production in B. subtilis: n = 71, p = 4088 genes (hdi::riboflavin).

What the lasso is

Ordinary least squares chooses coefficients to minimize the squared error $\sum_i (y_i - x_i'\beta)^2$. With $p=4088$ predictors and only $n=71$ samples that is hopeless — infinitely many coefficient vectors fit perfectly. The lasso (Tibshirani 1996) adds a penalty on the sum of absolute coefficients: $$ \min_{\beta}\ \underbrace{\sum_{i}(y_i - x_i'\beta)^2}_{\text{fit}}\ +\ \lambda\underbrace{\sum_{j}|\beta_j|}_{L_1\text{ penalty}}. $$

The entire trick is the absolute value. Because $|\beta_j|$ has a sharp corner at zero, the optimum lands exactly on zero for most coefficients — so the lasso does not merely shrink, it selects: most genes get coefficient 0 and drop out, a handful stay nonzero. (Ridge regression penalizes $\beta_j^2$ instead, which shrinks smoothly but never gives exact zeros — no selection.)

• $\lambda$ is the dial. $\lambda=0$ is ordinary regression (no selection); larger $\lambda$ → harsher penalty → fewer genes survive.
• Cross-validation picks $\lambda$. cv.glmnet splits the samples into 10 folds, refits at a grid of $\lambda$ values, and keeps the $\lambda$ with the best held-out prediction error (lambda.min).

How this relates to the Bayesian horseshoe. The lasso is exactly the posterior mode under a Laplace (double-exponential) prior on $\beta$ — a shrinkage method, but a single frequentist point estimate rather than a full posterior. Both give sparse $p\gg n$ selection, but differently: the lasso returns a hard in/out list under one global $\lambda$ and tends to over-select at the CV-optimal value; the horseshoe gives graded, per-gene evidence (the shrinkage weight $\kappa_j$) with heavy tails that let strong signals escape while crushing noise harder — hence sparser. Their agreeing on the same top gene despite these different philosophies is what makes this a convincing cross-check. (One caveat shared by both: with correlated predictors the lasso arbitrarily keeps one of a correlated group.)

.libPaths('C:/Users/user/R/win-library/4.6'); suppressMessages(library(glmnet))
d <- read.csv('ribo_data.csv', check.names=FALSE)
y <- d$y; X <- as.matrix(d[,-1]); colnames(X) <- sub('^X\\.','',colnames(X))

set.seed(1); cv <- cv.glmnet(X, y, alpha=1)                        # 10-fold cross-validated lasso
co <- coef(cv, s='lambda.min'); co <- co[co[,1]!=0, ]; co <- co[names(co)!='(Intercept)']
co <- co[order(-abs(co))]
cat(sprintf('lasso (lambda.min) selects %d of %d genes\n', length(co), ncol(X)))
cat('top genes by |coefficient|:\n'); for (i in 1:8) cat(sprintf('  %-12s %+.3f\n', names(co)[i], co[i]))
cat('\ncross-check -> Python horseshoe: YOAB_at is the #1 gene in both selectors; their secondary genes differ, as expected of two different selectors.\n')
cat('(the lasso keeps ~41 genes; the horseshoe is sparser -- lambda.min lasso is known to over-select.)\n')

options(repr.plot.width=11, repr.plot.height=4.2); par(mfrow=c(1,2))
plot(cv)                                                            # CV error vs penalty
barplot(rev(abs(co[1:10])), horiz=TRUE, las=1, col='#2c6fbb', xlab='|lasso coefficient|',
        names.arg=rev(names(co)[1:10]), main='Top genes selected by the lasso', cex.names=0.8)

lasso (lambda.min) selects 41 of 4088 genes
top genes by |coefficient|:
  YOAB_at      -0.812
  YEBC_at      -0.535
  LYSC_at      -0.298
  SPOVAA_at    +0.265
  YQJU_at      +0.232
  YXLD_at      -0.201
  YCLB_at      +0.199
  ARGF_at      -0.191

cross-check -> Python horseshoe: YOAB_at is the #1 gene in both selectors; their secondary genes differ, as expected of two different selectors.
(the lasso keeps ~41 genes; the horseshoe is sparser -- lambda.min lasso is known to over-select.)

Same top gene, a different penalty

• The lasso puts YOAB_at first, matching the horseshoe's strongest signal — an $L_1$ penalty and a global-local Bayesian prior agreeing on the driver of riboflavin production is a strong cross-check across two very different philosophies.
• Sparsity differs, as expected. Cross-validated lasso at lambda.min keeps ~41 genes — it is known to over-select (the penalty that predicts best is not the one that recovers the sparsest true model), whereas the horseshoe concentrates its posterior on the two clear signals. The horseshoe's shrinkage weights give a graded measure of evidence per gene, where the lasso gives a hard in/out list.

References

• Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. JRSS B 58, 267–288.
• Friedman, J., Hastie, T. & Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. J. Statistical Software 33(1).
• Bühlmann, P., Kalisch, M. & Meier, L. (2014). High-dimensional statistics with a view toward applications in biology. Annual Review of Statistics 1, 255–278.