| SOCR ≫ | DSPA ≫ | DSPA3 Topics ≫ |
DSPA3: Data Science and Predictive Analytics (UMich HS650)
Chapter 11: Variable Importance and Controlled Feature
Selection
Chapter 11: Variable Importance and Controlled Feature Selection
SOCR/MCAIM (Ivo Dinov)
SOCR/MCAIM (Ivo Dinov)
August 28, 2026
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(plotly) # interactive figures and ALL 3-D graphics
library(glmnet)
library(rsample)How this chapter uses graphics
Every two-dimensional figure is drawn with
ggplot2and rendered statically. Immediately after each one, the equivalentplot_ly()code appears in a chunk markedeval=FALSE, echo=TRUE.Every three-dimensional figure is drawn with
plot_ly()and evaluated. Selection is governed throughout by two-parameter families — shrinkage over (coefficient, \(\lambda\)), recovery over (sample size, sparsity), error control over (target FDR, signal amplitude), and reading one of these from a fixed viewpoint discards the interaction that is the whole subject.
1 Learning objectives
After completing this chapter you will be able to:
- Distinguish feature selection from dimension reduction, and filter, wrapper, and embedded strategies from one another.
- Recognize and avoid selection bias, and demonstrate that selecting on all the data can produce excellent cross-validated accuracy on pure noise.
- Derive the ridge and LASSO closed forms under orthonormal design and state the correct soft-threshold level.
- Explain the elastic net’s grouping effect and when it is preferable to the LASSO.
- Implement cyclic coordinate descent and validate it
against
glmnet. - State the irrepresentable condition and explain what it does and does not guarantee about support recovery.
- Use the exact degrees-of-freedom result for the LASSO to connect the path to information criteria.
- Apply Benjamini–Hochberg correctly as a step-up procedure, and explain why the rejection set is contiguous by construction.
- Run stability selection and state its bound on the expected number of false selections.
- Construct knockoffs, compute a valid \(W_j\) statistic, choose the threshold, and verify FDR control by replication.
Estimated time: 12–15 hours including exercises. Prerequisites: Chapter 3 (SVD, conditioning), Chapter 4 (dimension reduction, spurious correlation when \(p\gg n\)), Chapter 9 (optimism, nested resampling, effective degrees of freedom), and Chapter 2 (multiplicity and FDR).
2 PART I: THE SELECTION PROBLEM
3 Selection versus dimension reduction
Both reduce the number of columns entering a model, and they are not the same operation.
| Feature selection | Dimension reduction | |
|---|---|---|
| Process | Discrete — a subset is kept | Continuous — new axes are formed |
| Output | A subset \(S\subseteq\{1,\dots,p\}\) | \(k\) linear or nonlinear combinations |
| Interpretability | Original variables retained | Components mix all variables |
| Estimator variance | Higher — the subset itself is estimated | Lower — no combinatorial choice |
| Typical use | Biomarker discovery, cost reduction | Denoising, visualization, collinearity |
The variance row is the one most often overlooked. Because the selected set \(\hat S\) is itself a random function of the data, resampling produces different subsets, and any downstream quantity computed conditional on \(\hat S\) inherits that variability. Chapter 4 developed the continuous side; this chapter develops the discrete one, and much of it is about controlling exactly this extra variance.
4 Three families
\[ \begin{aligned} \textbf{Filter: }\quad &\text{score each feature (or set) by a statistic computed \emph{before} any model is fitted}\\ \textbf{Wrapper: }\quad &\text{search over subsets, scoring each by a fitted model's performance}\\ \textbf{Embedded: }\quad &\text{selection happens \emph{inside} the fitting procedure} \end{aligned} \]
| Family | Examples | Cost | Sees interactions? | Model-dependent? |
|---|---|---|---|---|
| Filter, univariate | correlation, \(t\)/\(\chi^2\), mutual information, SIS | \(O(np)\) | No | No |
| Filter, multivariate | CFS, mRMR, Markov blanket | \(O(np^2)\) | Partly | No |
| Wrapper, deterministic | forward/backward, RFE | \(O(p)\)–\(O(p^2)\) model fits | Yes | Yes |
| Wrapper, randomized | genetic algorithms, simulated annealing | Many fits | Yes | Yes |
| Embedded | LASSO, elastic net, tree importance, Boruta | One fit or path | Yes | Yes |
5 Selection bias: the cardinal error
Common misconception: “I cross-validated my model, so my accuracy estimate is honest.” If the features were selected using all the data, the cross-validation is contaminated before it begins. Every fold’s training set was chosen with knowledge of every fold’s held-out labels, so the selected features are already tuned to the assessment data.
The consequence is not subtle. With \(p\gg n\), selecting the most correlated features on the whole dataset and then cross-validating a classifier built from them yields excellent apparent accuracy on data with no signal whatsoever (Ambroise & McLachlan, 2002).
The correct protocol re-runs the entire selection procedure inside every training fold, so the assessment set is never consulted, the nested design of Chapter 9, §9.16, applied to selection rather than tuning.
set.seed(11)
# Data with NO relationship between X and y, by construction
simulate_null <- function(n = 100, p = 5000) {
list(X = matrix(rnorm(n * p), n, p),
y = factor(sample(rep(c("a", "b"), length.out = n))))
}
# WRONG: screen on all the data, then cross-validate the model
biased_cv <- function(d, k_keep = 20, folds = 5) {
scores <- abs(apply(d$X, 2, function(col) {
t.test(col ~ d$y)$statistic
}))
keep <- order(scores, decreasing = TRUE)[1:k_keep] # uses ALL labels
Xs <- d$X[, keep, drop = FALSE]
f <- sample(rep(seq_len(folds), length.out = nrow(Xs)))
mean(vapply(seq_len(folds), function(i) {
fit <- MASS::lda(Xs[f != i, , drop = FALSE], grouping = d$y[f != i])
mean(predict(fit, Xs[f == i, , drop = FALSE])$class == d$y[f == i])
}, numeric(1)))
}
# RIGHT: screen inside every training fold
honest_cv <- function(d, k_keep = 20, folds = 5) {
f <- sample(rep(seq_len(folds), length.out = nrow(d$X)))
mean(vapply(seq_len(folds), function(i) {
tr <- f != i
scores <- abs(apply(d$X[tr, , drop = FALSE], 2, function(col)
t.test(col ~ d$y[tr])$statistic))
keep <- order(scores, decreasing = TRUE)[1:k_keep] # TRAINING labels only
fit <- MASS::lda(d$X[tr, keep, drop = FALSE], grouping = d$y[tr])
mean(predict(fit, d$X[!tr, keep, drop = FALSE])$class == d$y[!tr])
}, numeric(1)))
}
set.seed(13)
res_bias <- replicate(20, { d <- simulate_null(); c(biased_cv(d), honest_cv(d)) })
c(truth = 0.5,
biased_screening = round(mean(res_bias[1, ]), 4),
honest_nested_screening = round(mean(res_bias[2, ]), 4),
inflation = round(mean(res_bias[1, ]) - mean(res_bias[2, ]), 4))#> truth biased_screening honest_nested_screening
#> 0.5000 0.8635 0.4985
#> inflation
#> 0.3650There is no signal in these data, the labels are a random permutation. The biased protocol reports accuracy far above chance; the nested protocol reports chance, which is correct.
ggplot(data.frame(acc = c(res_bias[1, ], res_bias[2, ]),
protocol = rep(c("Screen on all data, then CV",
"Screen inside each fold"), each = 20)),
aes(protocol, acc, fill = protocol)) +
geom_hline(yintercept = 0.5, linetype = "dashed", color = "firebrick") +
geom_boxplot(width = 0.5, alpha = 0.85, show.legend = FALSE) +
geom_jitter(width = 0.08, alpha = 0.5, size = 1.2, show.legend = FALSE) +
scale_fill_manual(values = c("#D8433B", "#3B7DD8")) +
coord_flip() +
labs(title = "Cross-validated accuracy on data containing no signal",
subtitle = "Dashed line: the truth. Only the nested protocol recovers it",
x = NULL, y = "CV accuracy") +
theme_dspa()
The inflation grows with the number of candidate features and shrinks with the sample size, which makes it a surface:
p_grid <- c(100, 500, 2000, 8000)
n_grid <- c(40, 60, 100, 160, 260)
bias_at <- function(n, p, reps = 6) {
set.seed(n * 7 + p)
mean(vapply(seq_len(reps), function(r) biased_cv(simulate_null(n, p)), numeric(1)))
}
Zbias <- outer(p_grid, n_grid, Vectorize(bias_at))
plot_ly(x = n_grid, y = p_grid, z = Zbias, type = "surface",
colorscale = "Inferno", reversescale = TRUE,
colorbar = list(title = "Apparent\nCV accuracy")) |>
add_trace(x = rep(n_grid, each = length(p_grid)),
y = rep(p_grid, times = length(n_grid)),
z = rep(0.5, length(n_grid) * length(p_grid)),
type = "scatter3d", mode = "markers", name = "Truth = 0.5",
marker = list(size = 2, color = "black")) |>
layout(title = "Selection bias: apparent accuracy on pure noise",
scene = list(xaxis = list(title = "Sample size n"),
yaxis = list(title = "Candidate features p", type = "log"),
zaxis = list(title = "Apparent CV accuracy")))The black plane marks the truth. The surface rises above it toward the large-\(p\), small-\(n\) corner, exactly the regime of genomics, imaging, and most biomarker discovery.
6 Univariate filters, and what they miss
Filters score each feature independently of any model. They are \(O(np)\), need no fitting, and scale to \(p\) in the millions, which makes them the standard first stage when \(p\) is enormous.
Their weakness follows from the same independence: a feature that is useless alone but informative in combination is invisible to them.
set.seed(21)
n_f <- 400
x1 <- rnorm(n_f); x2 <- rnorm(n_f)
noise <- matrix(rnorm(n_f * 10), n_f, 10)
# The signal is a pure INTERACTION: neither margin carries it
y_f <- factor(ifelse(x1 * x2 + rnorm(n_f, sd = 0.4) > 0, "a", "b"))
Xf <- cbind(x1 = x1, x2 = x2, noise)
colnames(Xf)[3:12] <- paste0("noise", 1:10)
univ <- apply(Xf, 2, function(col) abs(t.test(col ~ y_f)$statistic))
head(sort(univ, decreasing = TRUE), 5)#> noise10 noise1 noise2 noise6 noise4
#> 1.775270 1.104837 0.885433 0.867957 0.857797c(rank_of_x1 = which(names(sort(univ, decreasing = TRUE)) == "x1"),
rank_of_x2 = which(names(sort(univ, decreasing = TRUE)) == "x2"),
interaction_t = round(abs(t.test(x1 * x2 ~ y_f)$statistic), 2))#> rank_of_x1 rank_of_x2 interaction_t.t
#> 6.0 12.0 13.7data.frame(x1 = x1, x2 = x2, y = y_f) |>
ggplot(aes(x1, x2, color = y)) +
geom_point(size = 1.4, alpha = 0.8) +
scale_color_manual(values = c("#3B7DD8", "#D8433B")) +
coord_fixed() +
labs(title = "A signal no univariate filter can see",
subtitle = "The classes separate on the product x1*x2. Each margin alone is uninformative",
x = expression(x[1]), y = expression(x[2]), color = NULL) +
theme_dspa()
Both true predictors rank alongside pure noise, while the product has an enormous \(t\)-statistic. A univariate screen would discard exactly the two features that matter.
6.1 Sure independence screening
When \(p\) is far too large for any wrapper or path algorithm, screening is the only feasible first step. Sure independence screening (Fan & Lv, 2008) ranks by marginal correlation and retains the top \(d\):
\[\hat{\mathcal M}_\gamma=\Big\{1\le j\le p:\ |\hat\omega_j|\ \text{is among the largest } d\Big\}, \qquad \hat\omega_j=\mathbf{x}_j^\top\mathbf{y},\]
typically with \(d=\lfloor n/\log n\rfloor\).
The sure screening property. Under regularity conditions, most importantly that no truly important predictor is marginally uncorrelated with the response — \[P\big(\mathcal M_\star\subseteq\hat{\mathcal M}_\gamma\big)\to1\quad\text{as }n\to\infty,\] so the retained set contains all the true predictors with probability tending to 1. SIS is a screening step, not a selection step: it reduces \(p\) to a manageable size, after which a proper selector runs on the survivors.
set.seed(23)
n_s <- 200; p_s <- 20000
beta_s <- numeric(p_s); true_set <- sample(p_s, 8); beta_s[true_set] <- 2
Xs <- matrix(rnorm(n_s * p_s), n_s, p_s)
ys <- as.vector(Xs %*% beta_s) + rnorm(n_s)
d_sis <- floor(n_s / log(n_s))
omega <- abs(crossprod(Xs, ys))
screened <- order(omega, decreasing = TRUE)[1:d_sis]
c(p_original = p_s, d_retained = d_sis,
reduction_factor = round(p_s / d_sis),
true_predictors = length(true_set),
true_predictors_retained = sum(true_set %in% screened),
sure_screening_achieved = all(true_set %in% screened))#> p_original d_retained reduction_factor
#> 20000 37 541
#> true_predictors true_predictors_retained sure_screening_achieved
#> 8 8 1All eight true predictors survive a 200-fold reduction. That is what screening is for, it does not select, it makes selection possible.
The condition matters. SIS retains a predictor only if it is marginally associated with the outcome. The interaction example above is precisely the case where it fails: \(x_1\) and \(x_2\) have zero marginal correlation with \(y\), so no amount of data will screen them in. Iterated SIS and model-based screening exist for this reason.
7 PART II — REGULARIZED LINEAR MODELING
8 The framework
With \(\mathbf y\in\mathbb R^n\) and design \(X\in\mathbb R^{n\times p}\), the penalized estimator is
\[\boxed{\;\hat{\boldsymbol\beta}(\lambda)=\arg\min_{\boldsymbol\beta}\left\{\underbrace{\frac{1}{2n}\big\|\mathbf y-X\boldsymbol\beta\big\|_2^2}_{\text{fidelity}}+\underbrace{\lambda\,J(\boldsymbol\beta)}_{\text{regularizer}}\right\}\;}\]
with \(\lambda\ge0\) the regularization parameter and \(J\) a penalty. The intercept is conventionally left unpenalized.
Every choice of \(J\) gives a different estimator:
\[ \begin{aligned} J(\boldsymbol\beta)=\|\boldsymbol\beta\|_2^2 &\;\Longrightarrow\; \textbf{ridge} && \text{shrinks, never zeroes}\\ J(\boldsymbol\beta)=\|\boldsymbol\beta\|_1 &\;\Longrightarrow\; \textbf{LASSO} && \text{shrinks and zeroes}\\ J(\boldsymbol\beta)=\tfrac{1-\alpha}{2}\|\boldsymbol\beta\|_2^2+\alpha\|\boldsymbol\beta\|_1 &\;\Longrightarrow\; \textbf{elastic net} && \text{both}\\ J(\boldsymbol\beta)=\|\boldsymbol\beta\|_0 &\;\Longrightarrow\; \textbf{best subset} && \text{NP-hard} \end{aligned} \]
Equivalently, in constrained form. By Lagrangian duality each penalized problem corresponds to a constrained one,
\[\min_{\boldsymbol\beta}\ \|\mathbf y-X\boldsymbol\beta\|_2^2 \quad\text{subject to}\quad J(\boldsymbol\beta)\le t,\]
with a one-to-one decreasing correspondence between \(\lambda\) and \(t\). The constrained form is what makes the geometry visible (§11.8.2).
9 Ridge regression
\[\hat{\boldsymbol\beta}^{\text{ridge}}(\lambda)=\big(X^\top X+n\lambda I\big)^{-1}X^\top\mathbf y\]
The added \(n\lambda I\) makes the system invertible even when \(p>n\), where \(X^\top X\) is singular, this is the ill-posedness repair of Chapter 3, §3.7.
Under orthonormal design (\(X^\top X=I\)), with the \(\frac{1}{2n}\) fidelity convention, differentiating \(\frac{1}{2n}\|\mathbf y-X\boldsymbol\beta\|^2+\lambda\|\boldsymbol\beta\|_2^2\) coordinatewise gives
\[-\frac{1}{n}\big(\hat\beta^{OLS}_j-\beta_j\big)+2\lambda\beta_j=0 \quad\Longrightarrow\quad \boxed{\;\hat\beta_j^{\text{ridge}}=\frac{\hat\beta_j^{OLS}}{1+2n\lambda}\;}\]
every coefficient is scaled by the same factor, so none is ever exactly zero. Ridge does not select.
Effective degrees of freedom. With the SVD \(X=UDV^\top\) and singular values \(d_j\),
\[\mathrm{df}(\lambda)=\operatorname{tr}\big(X(X^\top X+n\lambda I)^{-1}X^\top\big)=\sum_{j=1}^{p}\frac{d_j^2}{d_j^2+n\lambda},\]
a continuous function decreasing from \(p\) to 0, the quantity developed in Chapter 9, §9.3.
set.seed(31)
n_o <- 200; p_o <- 5
Q <- qr.Q(qr(matrix(rnorm(n_o * p_o), n_o, p_o))) # orthonormal columns
beta_true <- c(3, -2, 1.5, 0, 0)
y_o <- as.vector(Q %*% beta_true) + rnorm(n_o, sd = 0.5)
beta_ols <- as.vector(crossprod(Q, y_o))
lam <- 0.02
fit_r <- glmnet(Q, y_o, alpha = 0, lambda = lam, standardize = FALSE,
intercept = FALSE)
data.frame(
OLS = round(beta_ols, 5),
glmnet_ridge = round(as.vector(coef(fit_r))[-1], 5),
closed_form = round(beta_ols / (1 + 2 * n_o * lam), 5))10 The LASSO
\[\hat{\boldsymbol\beta}^{L}(\lambda)=\arg\min_{\boldsymbol\beta}\left\{\frac{1}{2n}\|\mathbf y-X\boldsymbol\beta\|_2^2+\lambda\|\boldsymbol\beta\|_1\right\}\]
10.1 The soft threshold, derived
\(|\beta_j|\) is not differentiable at 0, so the stationarity condition is a subgradient condition. Under orthonormal design the objective separates across coordinates, and for coordinate \(j\),
\[0\in-\frac1n\big(\hat\beta_j^{OLS}-\beta_j\big)+\lambda\,\partial|\beta_j|, \qquad \partial|\beta_j|=\begin{cases}\{\operatorname{sign}(\beta_j)\}&\beta_j\ne0\\ [-1,1]&\beta_j=0.\end{cases}\]
For \(\beta_j>0\) this gives \(\beta_j=\hat\beta_j^{OLS}-n\lambda\), valid while the right side is positive. For \(\beta_j<0\), symmetrically. And \(\beta_j=0\) is optimal precisely when \(|\hat\beta_j^{OLS}|\le n\lambda\). Combining,
\[\boxed{\;\hat\beta_j^{L}=\mathcal S_{n\lambda}\big(\hat\beta_j^{OLS}\big)=\operatorname{sign}\big(\hat\beta_j^{OLS}\big)\Big(\big|\hat\beta_j^{OLS}\big|-n\lambda\Big)_{+}\;}\]
The threshold level depends on the fidelity convention, and the two must match. With \(\frac{1}{2n}\|\cdot\|^2\) the threshold is \(n\lambda\); with \(\frac{1}{n}\|\cdot\|^2\) it is \(\frac{n\lambda}{2}\); with an unscaled \(\|\cdot\|^2\) it is \(\lambda/2\). Quoting a ridge shrinkage factor derived under one convention beside a LASSO threshold derived under another gives two results that cannot both be right.
soft_threshold <- function(z, t) sign(z) * pmax(abs(z) - t, 0)
lam_l <- 0.01
fit_l <- glmnet(Q, y_o, alpha = 1, lambda = lam_l, standardize = FALSE,
intercept = FALSE)
data.frame(
OLS = round(beta_ols, 5),
glmnet_lasso = round(as.vector(coef(fit_l))[-1], 5),
closed_form = round(soft_threshold(beta_ols, n_o * lam_l), 5))10.2 Soft versus hard thresholding
\[ \begin{aligned} \textbf{Soft (}\ell_1\textbf{): }\quad &\mathcal S_t(z)=\operatorname{sign}(z)(|z|-t)_+ && \text{continuous; shrinks survivors}\\ \textbf{Hard (}\ell_0\textbf{): }\quad &\mathcal H_t(z)=z\,\mathbb 1\{|z|>t\} && \text{discontinuous; leaves survivors unchanged} \end{aligned} \]
Best subset and orthogonal matching pursuit are different algorithms. Best subset solves the \(\ell_0\) problem exactly, evaluating up to \(2^p\) subsets, NP-hard in general, feasible to about \(p\approx40\) with modern branch-and-bound. OMP is a greedy procedure: at each step add the predictor most correlated with the current residual, then refit. It runs in polynomial time and is not guaranteed to reach the \(\ell_0\) optimum. They coincide under orthogonal design, where hard thresholding solves both.
z_seq <- seq(-3, 3, length.out = 600); t_thr <- 1
data.frame(z = rep(z_seq, 3),
out = c(z_seq, soft_threshold(z_seq, t_thr),
z_seq * (abs(z_seq) > t_thr)),
op = rep(c("Identity (OLS)", "Soft (LASSO)", "Hard (best subset)"),
each = length(z_seq))) |>
ggplot(aes(z, out, color = op)) +
geom_hline(yintercept = 0, color = "grey75") +
geom_line(linewidth = 1) +
scale_color_manual(values = c("Identity (OLS)" = "grey45",
"Soft (LASSO)" = "#3B7DD8",
"Hard (best subset)" = "#D8433B")) +
coord_fixed() +
labs(title = "Thresholding operators at level t = 1",
subtitle = "Soft is continuous and shrinks survivors; hard jumps and leaves them untouched",
x = expression(hat(beta)^{OLS}), y = expression(hat(beta)), color = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(x = z_seq, y = z_seq, type = "scatter", mode = "lines",
name = "Identity") |>
add_lines(y = soft_threshold(z_seq, t_thr), name = "Soft (LASSO)") |>
add_lines(y = z_seq * (abs(z_seq) > t_thr), name = "Hard (best subset)") |>
layout(title = "Thresholding operators",
xaxis = list(title = "OLS coefficient", scaleanchor = "y"),
yaxis = list(title = "Penalized coefficient"))Shrinkage is a two-variable phenomenon, it depends on the OLS coefficient and on \(\lambda\):
b_grid <- seq(-3, 3, length.out = 80)
lam_grid <- seq(0.001, 1.2, length.out = 80)
Z_ridge <- outer(lam_grid, b_grid, function(l, b) b / (1 + 2 * l))
Z_lasso <- outer(lam_grid, b_grid, function(l, b) sign(b) * pmax(abs(b) - l, 0))
plot_ly() |>
add_surface(x = b_grid, y = lam_grid, z = Z_ridge, opacity = 0.85,
showscale = FALSE, colorscale = "Blues", name = "Ridge") |>
add_surface(x = b_grid, y = lam_grid, z = Z_lasso - 7, opacity = 0.85,
showscale = FALSE, colorscale = "Reds", name = "LASSO (offset)") |>
layout(title = "Shrinkage surfaces: ridge (upper, blue) and LASSO (lower, red, offset)",
scene = list(xaxis = list(title = "OLS coefficient"),
yaxis = list(title = "lambda"),
zaxis = list(title = "Penalized coefficient")))Rotate to look along the \(\lambda\) axis. The ridge surface is a smooth fan, every slice is a line through the origin with shrinking slope, so coefficients approach zero without reaching it. The LASSO surface has a flat plateau of exact zeros that widens with \(\lambda\), and the surviving coefficients sit on planes translated toward zero. That plateau is the selection.
11 Elastic net
\[\hat{\boldsymbol\beta}^{EN}=\arg\min_{\boldsymbol\beta}\left\{\frac{1}{2n}\|\mathbf y-X\boldsymbol\beta\|_2^2+\lambda\Big[\frac{1-\alpha}{2}\|\boldsymbol\beta\|_2^2+\alpha\|\boldsymbol\beta\|_1\Big]\right\}\]
Both parameters are needed and they do different jobs. \(\alpha\in[0,1]\) sets the mix, \(\alpha=0\) is ridge, \(\alpha=1\) is LASSO, while \(\lambda\) sets the total weight of the penalty. Varying \(\alpha\) with no \(\lambda\) changes the shape of the constraint region while leaving its size undetermined.
11.1 The grouping effect
Common misconception: “LASSO picks the important variables.” With a group of highly correlated predictors that are all genuinely associated with the outcome, LASSO tends to select one arbitrarily and zero the rest. Which one it picks is unstable, it changes under resampling, and even under a small perturbation of the data.
The elastic net’s \(\ell_2\) component creates a grouping effect: strongly correlated predictors receive similar coefficients and are selected or dropped together. When the scientific question is “which pathway matters” rather than “which single probe matters,” that is the behavior you want.
set.seed(41)
n_g <- 120
z1 <- rnorm(n_g)
# Three near-copies of the same underlying signal
Xg <- cbind(g1 = z1 + rnorm(n_g, sd = 0.05),
g2 = z1 + rnorm(n_g, sd = 0.05),
g3 = z1 + rnorm(n_g, sd = 0.05),
matrix(rnorm(n_g * 6), n_g, 6,
dimnames = list(NULL, paste0("noise", 1:6))))
yg <- 3 * z1 + rnorm(n_g, sd = 1)
round(cor(Xg[, 1:3]), 3)#> g1 g2 g3
#> g1 1.000 0.998 0.998
#> g2 0.998 1.000 0.998
#> g3 0.998 0.998 1.000alphas <- c(1, 0.5, 0.1)
grp <- bind_rows(lapply(alphas, function(a) {
set.seed(43)
cv <- cv.glmnet(Xg, yg, alpha = a, nfolds = 10)
cf <- as.vector(coef(cv, s = "lambda.min"))[-1]
data.frame(feature = colnames(Xg), coef = cf,
alpha = sprintf("alpha = %.1f%s", a,
ifelse(a == 1, " (LASSO)", "")))
}))
grp |> filter(grepl("^g", feature)) |>
pivot_wider(names_from = alpha, values_from = coef) |>
mutate(across(where(is.numeric), \(z) round(z, 4)))ggplot(grp, aes(feature, coef, fill = alpha)) +
geom_col(position = "dodge") +
scale_fill_manual(values = c("#D8433B", "#7FB069", "#3B7DD8")) +
labs(title = "The grouping effect",
subtitle = "g1, g2, g3 are correlated at r > 0.99. LASSO keeps one; lowering alpha spreads the weight",
x = NULL, y = "Coefficient at lambda.min", fill = NULL) +
theme_dspa(10) +
theme(axis.text.x = element_text(angle = 30, hjust = 1))
11.2 The geometry
In constrained form the estimator is where the smallest SSE contour first touches the feasible region \(\{J(\boldsymbol\beta)\le t\}\). The \(\ell_1\) ball has corners on the axes; a contour is far more likely to touch a corner than a smooth arc, and a corner has a coordinate exactly zero. The \(\ell_2\) ball has no corners, so ridge solutions are almost surely interior to every axis.
if (ONLINE) {
mlb <- dspa_read("https://umich.instructure.com/files/330381/download?download_frd=1",
"mlb.txt", reader = utils::read.table, header = TRUE)
} else mlb <- NULL
if (is.null(mlb)) {
set.seed(45); n_m <- 1034
mlb <- data.frame(Height = rnorm(n_m, 73, 2.3), Weight = rnorm(n_m, 202, 21),
Age = rnorm(n_m, 28, 4.3))
mlb$Height <- 60 + 0.05 * mlb$Weight + 0.02 * mlb$Age + rnorm(n_m, sd = 1.8)
}
Ym <- as.vector(scale(mlb$Height))
Xm <- scale(as.matrix(mlb[, c("Weight", "Age")]))
b_ols <- as.vector(solve(crossprod(Xm), crossprod(Xm, Ym)))
sse_grid <- expand.grid(b1 = seq(-0.3, 0.9, length.out = 220),
b2 = seq(-0.5, 0.4, length.out = 220))
sse_grid$sse <- apply(as.matrix(sse_grid[, 1:2]), 1,
function(b) sum((Ym - Xm %*% b)^2))
# Constraint boundaries at a COMMON radius t, for several alpha
t_rad <- 0.30
theta <- seq(0, 2 * pi, length.out = 800)
en_boundary <- function(alpha, t) {
# Solve (1-alpha)/2 * r^2 * m2 + alpha * r * m1 = t along each ray
ux <- cos(theta); uy <- sin(theta)
m1 <- abs(ux) + abs(uy); m2 <- ux^2 + uy^2
r <- if (alpha == 1) t / m1 else
(-alpha * m1 + sqrt((alpha * m1)^2 + 2 * (1 - alpha) * m2 * t)) /
((1 - alpha) * m2)
data.frame(b1 = r * ux, b2 = r * uy, alpha = sprintf("alpha = %.2f", alpha))
}
bounds <- bind_rows(lapply(c(0, 0.25, 0.5, 0.75, 1), en_boundary, t = t_rad))
### OLD, buggs
# ggplot() +
# geom_contour(data = sse_grid, aes(b1, b2, z = sse),
# breaks = quantile(sse_grid$sse, seq(0.002, 0.25, length.out = 9)),
# color = "grey55", linewidth = 0.35) +
# geom_path(data = bounds, aes(b1, b2), color = "#D8433B", linewidth = 0.8) +
# geom_point(aes(b_ols[1], b_ols[2]), size = 3, color = "black") +
# annotate("text", x = b_ols[1], y = b_ols[2] - 0.06, label = "OLS", size = 3.2) +
# geom_hline(yintercept = 0, color = "grey70", linewidth = 0.3) +
# geom_vline(xintercept = 0, color = "grey70", linewidth = 0.3) +
# facet_wrap(~ alpha, nrow = 2) + coord_fixed() +
# labs(title = "Elastic-net feasible regions at a common radius, with SSE contours",
# subtitle = "Grey: SSE contours centered at OLS. Red: the constraint boundary. Corners appear as alpha increases",
# x = expression(beta[1]~"(Weight)"), y = expression(beta[2]~"(Age)")) +
# theme_dspa(9)
alpha_vals <- c(0, 0.25, 0.5, 0.75, 1)
alpha_labs <- sprintf("alpha = %.2f", alpha_vals)
# Repeat SSE grid for each facet
sse_grid_facet <- do.call(rbind, lapply(alpha_labs, function(a) {
d <- sse_grid
d$alpha <- a
d
}))
# OLS point in every facet
ols_pt <- data.frame(
b1 = b_ols[1],
b2 = b_ols[2],
alpha = alpha_labs
)
# For the text label, shift y slightly
ols_text <- transform(ols_pt, b2 = b2 - 0.06)
ggplot() +
geom_contour(data = sse_grid_facet, aes(b1, b2, z = sse),
breaks = quantile(sse_grid$sse, seq(0.002, 0.25, length.out = 9)),
color = "grey55", linewidth = 0.35) +
geom_path(data = bounds, aes(b1, b2), color = "#D8433B", linewidth = 0.8) +
geom_point(data = ols_pt, aes(b1, b2), size = 3, color = "black") +
geom_text(data = ols_text, aes(b1, b2), label = "OLS", size = 3.2) +
geom_hline(yintercept = 0, color = "grey70", linewidth = 0.3) +
geom_vline(xintercept = 0, color = "grey70", linewidth = 0.3) +
facet_wrap(~ alpha, nrow = 2) +
coord_fixed() +
labs(
title = "Elastic-net feasible regions at a common radius, with SSE contours",
subtitle = "Grey: SSE contours centered at OLS. Red: the constraint boundary. Corners appear as alpha increases",
x = expression(beta[1] ~ "(Weight)"),
y = expression(beta[2] ~ "(Age)")
) +
theme_dspa(9)
All five panels use the same constraint radius, so the only thing changing is the boundary’s shape. At \(\alpha=0\) it is a circle with no corners; by \(\alpha=1\) it is a diamond whose vertices sit on the axes, and a contour touching a vertex sets a coefficient to exactly zero.
a_grid <- seq(0, 1, length.out = 40)
th <- seq(0, 2 * pi, length.out = 120)
Rad <- outer(a_grid, th, Vectorize(function(a, tt) {
ux <- cos(tt); uy <- sin(tt)
m1 <- abs(ux) + abs(uy); m2 <- ux^2 + uy^2
if (a == 1) t_rad / m1 else
(-a * m1 + sqrt((a * m1)^2 + 2 * (1 - a) * m2 * t_rad)) / ((1 - a) * m2)
}))
Xs_ <- Rad * matrix(cos(th), nrow(Rad), ncol(Rad), byrow = TRUE)
Ys_ <- Rad * matrix(sin(th), nrow(Rad), ncol(Rad), byrow = TRUE)
Zs_ <- matrix(a_grid, nrow(Rad), ncol(Rad))
plot_ly(x = Xs_, y = Ys_, z = Zs_, type = "surface",
surfacecolor = Zs_, colorscale = "Portland",
colorbar = list(title = "alpha")) |>
layout(title = "The elastic-net constraint boundary as alpha sweeps 0 (circle) to 1 (diamond)",
scene = list(xaxis = list(title = "beta1"),
yaxis = list(title = "beta2"),
zaxis = list(title = "alpha")))Rotate to look down the \(\alpha\) axis: the cross-section morphs continuously from a circle to a diamond, and the corners emerge rather than appearing abruptly. That continuous emergence is why intermediate \(\alpha\) gives partial selection with grouping.
12 Coordinate descent
The LASSO objective is convex but non-smooth, so gradient descent does not apply directly. Cyclic coordinate descent exploits the fact that each coordinate’s subproblem has the closed-form soft-threshold solution:
\[\beta_j\ \leftarrow\ \frac{1}{\|\mathbf x_j\|^2}\,\mathcal S_{n\lambda}\!\left(\mathbf x_j^\top\mathbf r^{(-j)}\right), \qquad \mathbf r^{(-j)}=\mathbf y-\sum_{k\ne j}\mathbf x_k\beta_k,\]
cycling over \(j=1,\dots,p\) until convergence. The partial residual \(\mathbf r^{(-j)}\) is updated incrementally rather than recomputed, which makes each sweep \(O(np)\).
Common misconception: “the LASSO is just soft-thresholded least squares.” Applying \(\mathcal S_{n\lambda}\) once to the OLS solution gives the LASSO only under orthonormal design, where the objective separates across coordinates and each one can be solved in isolation.
With correlated predictors the coordinates are coupled: shrinking \(\beta_j\) changes the partial residual seen by \(\beta_k\), which changes its update, which feeds back. The estimate is reached only by cycling to convergence. A single pass produces a different estimator, one with no optimality property and no name.
lasso_cd <- function(X, y, lambda, tol = 1e-10, max_iter = 10000) {
X <- scale(X, center = TRUE, scale = FALSE)
y <- y - mean(y)
n <- nrow(X); p <- ncol(X)
xx <- colSums(X^2)
beta <- numeric(p)
r <- y - as.vector(X %*% beta) # current residual
for (it in seq_len(max_iter)) {
beta_old <- beta
for (j in seq_len(p)) {
r <- r + X[, j] * beta[j] # add coordinate j back in
beta[j] <- soft_threshold(sum(X[, j] * r), n * lambda) / xx[j]
r <- r - X[, j] * beta[j] # remove the updated coordinate
}
if (max(abs(beta - beta_old)) < tol) break
}
list(beta = beta, iterations = it)
}
set.seed(51)
n_c <- 150; p_c <- 12
Sig <- 0.7^abs(outer(1:p_c, 1:p_c, "-")) # strongly correlated design
Xc <- MASS::mvrnorm(n_c, rep(0, p_c), Sig)
bc <- c(2, -1.5, 1, rep(0, p_c - 3))
yc <- as.vector(Xc %*% bc) + rnorm(n_c)
lam_c <- 0.08
manual <- lasso_cd(Xc, yc, lam_c)
gl <- as.vector(coef(glmnet(Xc, yc, alpha = 1, lambda = lam_c,
standardize = FALSE, thresh = 1e-14))[-1])
one_pass <- soft_threshold(
as.vector(solve(crossprod(Xc), crossprod(Xc, yc))), n_c * lam_c)
data.frame(
coordinate_descent = round(manual$beta, 5),
glmnet = round(gl, 5),
single_threshold_pass = round(one_pass, 5)) |> head(6)c(iterations_to_converge = manual$iterations,
max_abs_diff_vs_glmnet = signif(max(abs(manual$beta - gl)), 3),
max_abs_diff_one_pass = signif(max(abs(one_pass - gl)), 3))#> iterations_to_converge max_abs_diff_vs_glmnet max_abs_diff_one_pass
#> 4.80e+01 1.88e-07 1.58e+00Coordinate descent matches glmnet to eleven decimal
places. The single thresholding pass does not, on a
correlated design it is a different estimator entirely.
12.1 Solution paths and complexity
glmnet computes the whole path over a decreasing grid of
\(\lambda\), warm-starting each fit
from the previous solution. Two devices make it fast: warm
starts mean few sweeps per \(\lambda\), and active-set
iteration restricts the inner loop to currently non-zero
coordinates.
| Algorithm | Cost | Note |
|---|---|---|
| Coordinate descent (path of \(K\) values) | \(O(npK)\) with warm starts | glmnet’s default |
| LARS | \(O(p^3+np^2)\) | Exact piecewise-linear path |
| Best subset | \(O(2^p)\) | Exact \(\ell_0\); branch-and-bound to \(p\approx40\) |
| OMP | \(O(nps)\) for \(s\) steps | Greedy approximation |
| Ridge via SVD | \(O(np\min(n,p))\) | One decomposition serves all \(\lambda\) |
fit_path_l <- glmnet(Xc, yc, alpha = 1)
fit_path_r <- glmnet(Xc, yc, alpha = 0)
path_df <- function(fit, label) {
as.data.frame(as.matrix(t(fit$beta))) |>
mutate(loglambda = log(fit$lambda)) |>
pivot_longer(-loglambda, names_to = "variable", values_to = "coef") |>
mutate(model = label)
}
bind_rows(path_df(fit_path_l, "LASSO"), path_df(fit_path_r, "Ridge")) |>
ggplot(aes(loglambda, coef, color = variable)) +
geom_hline(yintercept = 0, color = "grey70") +
geom_line(linewidth = 0.7, show.legend = FALSE) +
facet_wrap(~ model) +
labs(title = "Regularization paths on a correlated design",
subtitle = "LASSO coefficients hit zero and stay there; ridge coefficients approach zero without arriving",
x = expression(log(lambda)), y = "Coefficient") +
theme_dspa(10)
# --- Interactive equivalent ------------------------------------------------
d <- path_df(fit_path_l, "LASSO")
plot_ly(d, x = ~loglambda, y = ~coef, color = ~variable,
type = "scatter", mode = "lines") |>
layout(title = "LASSO solution path",
xaxis = list(title = "log(lambda)"),
yaxis = list(title = "Coefficient"))12.2 Choosing \(\lambda\)
# One helper, used everywhere. The y-axis label follows type.measure, and the
# nzero annotations are attached as a proper column.
plot_cv_glmnet <- function(cvfit, title = "") {
d <- data.frame(loglam = log(cvfit$lambda), m = cvfit$cvm,
lo = cvfit$cvlo, hi = cvfit$cvup, nzero = cvfit$nzero)
ggplot(d, aes(loglam, m)) +
geom_ribbon(aes(ymin = lo, ymax = hi), fill = "grey86") +
geom_line(linewidth = 0.9, color = "steelblue") +
geom_point(size = 1.1) +
geom_vline(xintercept = log(cvfit$lambda.min), linetype = "dashed",
color = "firebrick") +
geom_vline(xintercept = log(cvfit$lambda.1se), linetype = "dotted",
color = "grey25") +
labs(title = paste0("Cross-validated ", cvfit$name,
if (nzchar(title)) paste0(" (", title, ")") else ""),
subtitle = sprintf("Dashed: lambda.min (%d active). Dotted: lambda.1se (%d active)",
cvfit$nzero[which(cvfit$lambda == cvfit$lambda.min)],
cvfit$nzero[which(cvfit$lambda == cvfit$lambda.1se)]),
x = expression(log(lambda)), y = cvfit$name) +
theme_dspa(10)
}
c(lambda_min = signif(cv_l$lambda.min, 4),
active_at_min = cv_l$nzero[which(cv_l$lambda == cv_l$lambda.min)],
lambda_1se = signif(cv_l$lambda.1se, 4),
active_at_1se = cv_l$nzero[which(cv_l$lambda == cv_l$lambda.1se)],
true_nonzero = sum(bc != 0))#> lambda_min active_at_min.s42 lambda_1se active_at_1se.s29
#> 0.02715 8.00000 0.09099 4.00000
#> true_nonzero
#> 3.00000lambda.min minimizes the CV error;
lambda.1se is the largest \(\lambda\) whose CV error is within one
standard error of that minimum, the 1-SE rule of Chapter 9, §9.21.2, trading a statistically
indistinguishable amount of accuracy for a sparser, more stable model.
For selection purposes lambda.1se is
usually the better choice; for prediction,
lambda.min.
13 PART III — WHEN DOES SELECTION SUCCEED?
The LASSO produces a subset. Whether that subset is the true support is a separate question with a precise answer.
14 The irrepresentable condition
Let \(S=\{j:\beta_j^\star\ne0\}\) be the true support and \(S^c\) its complement. Write \(X_S\) and \(X_{S^c}\) for the corresponding column blocks.
Irrepresentable condition (Zhao & Yu, 2006; Zou, 2006). The LASSO is sign-consistent, it recovers \(S\) with the correct signs, with probability tending to 1, essentially if and only if \[\boxed{\;\Big\|\,X_{S^c}^\top X_S\big(X_S^\top X_S\big)^{-1}\operatorname{sign}\big(\boldsymbol\beta^\star_S\big)\,\Big\|_\infty\ <\ 1\;}\]
The interpretation is direct: the quantity inside is the vector of coefficients obtained by regressing each irrelevant predictor on the relevant ones, combined according to the true signs. The condition requires that no irrelevant predictor be too well “represented” by the signal variables. If some noise variable is nearly a linear combination of the true predictors, with the wrong sign pattern, the LASSO will select it no matter how much data you collect.
Common misconception: “with enough data the LASSO finds the right variables.” Consistency in prediction holds under mild conditions. Consistency in selection requires the irrepresentable condition, which is a property of the design matrix, not of the sample size. When it fails, more data does not help: the LASSO converges to the wrong support.
This is what motivates the adaptive LASSO (reweight the penalty by \(1/|\hat\beta_j^{\text{init}}|^\gamma\)), the relaxed LASSO, and stability selection, each of which weakens the requirement in a different way.
irrep <- function(X, support, signs) {
XS <- X[, support, drop = FALSE]; XSc <- X[, -support, drop = FALSE]
max(abs(crossprod(XSc, XS) %*% solve(crossprod(XS)) %*% signs))
}
set.seed(61)
n_i <- 400; p_i <- 20; s_i <- 3
supp <- 1:s_i; sgn <- rep(1, s_i)
# Two designs: benign (weakly correlated) and adversarial (a noise column
# nearly equal to the sum of the signal columns)
X_ok <- matrix(rnorm(n_i * p_i), n_i, p_i)
X_bad <- X_ok
X_bad[, s_i + 1] <- rowSums(X_ok[, supp]) / s_i + rnorm(n_i, sd = 0.05)
c(irrepresentable_benign = round(irrep(X_ok, supp, sgn), 4),
irrepresentable_adversarial = round(irrep(X_bad, supp, sgn), 4),
condition_holds_benign = irrep(X_ok, supp, sgn) < 1,
condition_holds_adversarial = irrep(X_bad, supp, sgn) < 1)#> irrepresentable_benign irrepresentable_adversarial
#> 0.1475 0.9976
#> condition_holds_benign condition_holds_adversarial
#> 1.0000 1.0000recover_rate <- function(X, n_rep = 40, seed = 63) {
set.seed(seed)
b <- numeric(ncol(X)); b[supp] <- 2
vapply(seq_len(n_rep), function(r) {
y <- as.vector(X %*% b) + rnorm(nrow(X))
cv <- cv.glmnet(X, y, alpha = 1, nfolds = 10)
sel <- which(as.vector(coef(cv, s = "lambda.min"))[-1] != 0)
identical(sort(sel), sort(supp))
}, logical(1)) |> mean()
}
data.frame(
design = c("Benign", "Adversarial"),
irrepresentable_value = round(c(irrep(X_ok, supp, sgn),
irrep(X_bad, supp, sgn)), 3),
exact_recovery_rate = c(recover_rate(X_ok), recover_rate(X_bad)))The adversarial design violates the condition, and the LASSO recovers the exact support far less often, on the same sample size, with the same signal strength. The failure is structural.
14.1 The phase transition
Support recovery has a sharp threshold in \((n, s)\): for a random Gaussian design, exact recovery requires roughly \(n\gtrsim 2s\log p\) samples.
n_pt <- c(40, 70, 110, 160, 220, 300)
s_pt <- c(2, 4, 6, 9, 13, 18)
p_pt <- 200
## OLD: outer(s_pt, n_pt, Vectorize(recovery_at)) passes the first argument from s_pt and the second from n_pt. Original function recovery_at expects the first argument to be n (sample size) and the second to be s (sparsity). Therefore, inside the function, n gets the small values (2, 4, …) and s gets the larger ones (40, 70, …). With n = 2, cv.glmnet has too few observations, and the response y can appear constant (or the standardization fails), leading to the error.
# recovery_at <- function(n, s, reps = 12) {
# set.seed(n * 31 + s)
# mean(vapply(seq_len(reps), function(r) {
# X <- matrix(rnorm(n * p_pt), n, p_pt)
# S <- sample(p_pt, s); b <- numeric(p_pt); b[S] <- 2.5
# y <- as.vector(X %*% b) + rnorm(n)
# cv <- cv.glmnet(X, y, alpha = 1, nfolds = 5)
# sel <- which(as.vector(coef(cv, s = "lambda.1se"))[-1] != 0)
# length(intersect(sel, S)) / s # proportion of true support recovered
# }, numeric(1)))
# }
recovery_at <- function(s, n, reps = 12) {
set.seed(n * 31 + s)
mean(vapply(seq_len(reps), function(r) {
X <- matrix(rnorm(n * p_pt), n, p_pt)
S <- sample(p_pt, s); b <- numeric(p_pt); b[S] <- 2.5
y <- as.vector(X %*% b) + rnorm(n)
cv <- cv.glmnet(X, y, alpha = 1, nfolds = 5)
sel <- which(as.vector(coef(cv, s = "lambda.1se"))[-1] != 0)
length(intersect(sel, S)) / s
}, numeric(1)))
}
Zpt <- outer(s_pt, n_pt, Vectorize(recovery_at))
plot_ly(x = n_pt, y = s_pt, z = Zpt, type = "surface",
colorscale = "Viridis",
colorbar = list(title = "Fraction of\ntrue support\nrecovered")) |>
layout(title = sprintf("LASSO support recovery over sample size and sparsity (p = %d)", p_pt),
scene = list(xaxis = list(title = "Sample size n"),
yaxis = list(title = "Sparsity s"),
zaxis = list(title = "Recovery fraction")))Rotate to see the ridge running diagonally: the boundary between recovery and failure tracks \(n\propto s\log p\), not \(n\propto s\). Doubling the number of true signals requires more than twice the data, because each additional signal must be distinguished from all \(p-s\) noise variables.
pt_long <- expand.grid(n = n_pt, s = s_pt)
pt_long$recovery <- as.vector(t(Zpt))
pt_long$threshold <- 2 * pt_long$s * log(p_pt)
ggplot(pt_long, aes(n, recovery, color = factor(s))) +
geom_line(linewidth = 0.9) + geom_point(size = 1.8) +
scale_color_viridis_d(option = "plasma", end = 0.9, name = "Sparsity s") +
labs(title = "Recovery against sample size, by sparsity",
subtitle = sprintf("The curves shift right as s grows, at roughly the 2s*log(p) rate (log p = %.1f)",
log(p_pt)),
x = "Sample size n", y = "Fraction of true support recovered") +
theme_dspa()
15 Degrees of freedom for the LASSO
Chapter 9, §9.3 defined effective degrees of freedom as \(\mathrm{df}=\frac{1}{\sigma^2}\sum_i\operatorname{Cov}(\hat y_i,y_i)\). For ridge this is a smooth function of \(\lambda\). For the LASSO there is a remarkable exact result.
Theorem (Zou, Hastie & Tibshirani, 2007). For the LASSO with a fixed \(\lambda\), \[\boxed{\;\mathrm{df}\big(\hat{\boldsymbol\beta}^{L}(\lambda)\big)=\mathbb E\Big[\#\{j:\hat\beta_j(\lambda)\ne0\}\Big]\;},\] the expected number of non-zero coefficients, exactly and without approximation.
This is not obvious. Selection is a discrete, data-dependent operation, so one might expect the degrees of freedom to exceed the number of selected variables, as it does for best subset selection, where the search cost adds degrees of freedom. For the LASSO the shrinkage applied to the survivors exactly compensates for the cost of having chosen them.
set.seed(71)
n_df <- 100; p_df <- 30
X_df <- scale(matrix(rnorm(n_df * p_df), n_df, p_df))
b_df <- c(rep(1.5, 5), rep(0, p_df - 5))
mu_df <- as.vector(X_df %*% b_df); sig_df <- 1
lam_seq <- exp(seq(log(0.6), log(0.01), length.out = 12))
reps_df <- 400
df_study <- do.call(rbind, lapply(lam_seq, function(l) {
out <- replicate(reps_df, {
y <- mu_df + rnorm(n_df, sd = sig_df)
f <- glmnet(X_df, y, alpha = 1, lambda = l, standardize = FALSE)
yh <- as.vector(predict(f, X_df))
c(nz = sum(as.vector(coef(f))[-1] != 0), y = y, yh = yh)
})
nz <- out[1, ]
ys <- out[2:(n_df + 1), ]; yhs <- out[(n_df + 2):(2 * n_df + 1), ]
cov_sum <- sum(vapply(seq_len(n_df), \(i) cov(yhs[i, ], ys[i, ]), numeric(1)))
data.frame(lambda = l,
df_covariance = cov_sum / sig_df^2,
df_expected_nonzero = mean(nz))
}))
df_study |> mutate(across(where(is.numeric), \(z) round(z, 3)))ggplot(df_study, aes(df_expected_nonzero, df_covariance)) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "firebrick") +
geom_point(size = 2.8, color = "steelblue") +
coord_fixed() +
labs(title = "Degrees of freedom equals the expected number of non-zero coefficients",
subtitle = "Dashed: the identity line. Vertical axis from the covariance definition; horizontal from counting",
x = "E[# non-zero]", y = expression(sum(Cov(hat(y)[i], y[i]))/sigma^2)) +
theme_dspa()
The points sit on the identity line. This makes \(C_p\), AIC, and BIC directly computable
along the LASSO path by substituting the number of active
coefficients, which is why glmnet can report them.
16 Post-selection inference
Common misconception: “the LASSO selected these variables, so I can refit OLS on them and report the usual p-values.” Those p-values are invalid. They are computed as though the model had been specified in advance, when in fact it was chosen because it fit these data well. The selection event conditions the distribution of the estimates, and ignoring it produces anti-conservative inference, often dramatically so.
set.seed(81)
naive_p <- replicate(300, {
n <- 60; p <- 40
X <- matrix(rnorm(n * p), n, p)
y <- rnorm(n) # NO signal whatsoever
cv <- cv.glmnet(X, y, alpha = 1, nfolds = 5)
sel <- which(as.vector(coef(cv, s = "lambda.1se"))[-1] != 0)
if (length(sel) == 0 || length(sel) >= n - 2) return(NA_real_)
min(summary(lm(y ~ X[, sel, drop = FALSE]))$coefficients[-1, 4])
}) |> na.omit()
c(selections_with_at_least_one_variable = length(naive_p),
fraction_with_p_below_0.05 = round(mean(naive_p < 0.05), 4),
nominal_rate = 0.05)#> selections_with_at_least_one_variable fraction_with_p_below_0.05
#> 6.00 1.00
#> nominal_rate
#> 0.05ggplot(data.frame(p = as.numeric(naive_p)), aes(p)) +
geom_histogram(bins = 25, fill = "steelblue", color = "white") +
geom_vline(xintercept = 0.05, linetype = "dashed", color = "firebrick") +
labs(title = "Naive p-values after LASSO selection, on data with no signal",
subtitle = "Under a valid procedure these would be uniform. They are concentrated near zero",
x = "Smallest p-value among selected variables", y = "Replications") +
theme_dspa()
There is no signal in these simulations, so a valid procedure would produce uniform p-values and reject at 5%. The observed rejection rate is far higher.
Three sound alternatives:
Selective inference. Condition on the selection
event. For the LASSO the event is a polyhedral constraint \(\{Ay\le b\}\), and the conditional
distribution of a linear contrast is a truncated
Gaussian with computable truncation limits (Lee et al., 2016).
selectiveInference::fixedLassoInf() implements this.
Data splitting. Select on one half, infer on the other. Simple, always valid, and costs statistical efficiency, the inference uses only half the data.
Stability selection or knockoffs. Control an error rate over the selection itself rather than attaching p-values to the selected coefficients (§11.16, §11.17).
SI_OK <- requireNamespace("selectiveInference", quietly = TRUE)
if (SI_OK) {
set.seed(83)
n_si <- 100; p_si <- 25
X_si <- scale(matrix(rnorm(n_si * p_si), n_si, p_si))
b_si <- c(rep(2, 3), rep(0, p_si - 3))
y_si <- as.vector(X_si %*% b_si) + rnorm(n_si)
lam_si <- cv.glmnet(X_si, y_si, alpha = 1)$lambda.1se
fit_si <- glmnet(X_si, y_si, alpha = 1, standardize = FALSE, thresh = 1e-12)
beta_si <- as.vector(coef(fit_si, s = lam_si, exact = TRUE,
x = X_si, y = y_si))[-1]
si <- selectiveInference::fixedLassoInf(
X_si, y_si, beta_si, lam_si * n_si, sigma = 1)
data.frame(variable = si$vars,
coefficient = round(si$coef0, 4),
selective_pvalue = signif(si$pv, 4),
ci_low = round(si$ci[, 1], 3), ci_high = round(si$ci[, 2], 3),
truly_nonzero = si$vars <= 3)
} else {
message("selectiveInference not installed; see the data-splitting alternative below.")
}# The assumption-light alternative: select on one half, infer on the other
set.seed(85)
n_ds <- 200; p_ds <- 30
X_ds <- scale(matrix(rnorm(n_ds * p_ds), n_ds, p_ds))
b_ds <- c(rep(1.8, 4), rep(0, p_ds - 4))
y_ds <- as.vector(X_ds %*% b_ds) + rnorm(n_ds)
half <- sample(n_ds, n_ds / 2)
cv_ds <- cv.glmnet(X_ds[half, ], y_ds[half], alpha = 1)
sel_ds <- which(as.vector(coef(cv_ds, s = "lambda.1se"))[-1] != 0)
fit_inf <- lm(y_ds[-half] ~ X_ds[-half, sel_ds, drop = FALSE])
data.frame(variable = sel_ds,
estimate = round(coef(fit_inf)[-1], 4),
p_value = signif(summary(fit_inf)$coefficients[-1, 4], 4),
truly_nonzero = sel_ds <= 4)Selection used only the first half; inference uses only the second. The p-values are valid because the second half played no part in choosing the model, at the cost of halving the data available for each task.
17 PART IV: ERROR-CONTROLLED SELECTION
Selecting features is easy; selecting them with a guarantee is the harder and more useful problem. Three frameworks, in increasing generality.
18 False discovery rate and Benjamini–Hochberg
For a selected set \(\hat S\), the false discovery proportion and its expectation are
\[\mathrm{FDP}(\hat S)=\frac{\#\{j\in\hat S:\ j\ \text{is null}\}}{\max(1,|\hat S|)}, \qquad \mathrm{FDR}(\hat S)=\mathbb E\big[\mathrm{FDP}(\hat S)\big].\]
Controlling FDR at \(q\) means that on average a fraction \(q\) of the selections are false, a far more useful guarantee than family-wise error control when hundreds or thousands of hypotheses are in play.
Common misconception: “reject every hypothesis whose p-value falls below its own threshold.” That is a step-down rule, and it is not Benjamini–Hochberg. BH is a step-up procedure: sort the p-values \(p_{(1)}\le\cdots\le p_{(m)}\), find \[\hat k=\max\Big\{i:\ p_{(i)}\le \frac{q\,i}{m}\Big\},\] and reject \(H_{(1)},\dots,H_{(\hat k)}\), all of them, including any \(p_{(i)}\) with \(i<\hat k\) that failed its own threshold.
The two rules coincide only when the passing set happens to be contiguous. When it is not, the step-down version rejects strictly fewer hypotheses and is therefore needlessly conservative, it discards real discoveries while claiming the same error guarantee.
pvals <- sort(c(0.9, 0.35, 0.01, 0.013, 0.014, 0.19, 0.35, 0.5,
0.63, 0.67, 0.75, 0.81, 0.01, 0.051))
m <- length(pvals); q <- 0.05
thresh <- q * seq_len(m) / m
passes <- pvals <= thresh
k_hat <- if (any(passes)) max(which(passes)) else 0L
data.frame(i = seq_len(m), p = pvals, threshold = round(thresh, 5),
below_own_threshold = passes,
rejected_by_BH = seq_len(m) <= k_hat)c(k_hat = k_hat,
rejected = k_hat,
bonferroni_threshold = round(q / m, 5),
rejected_by_bonferroni = sum(pvals <= q / m))#> k_hat rejected bonferroni_threshold
#> 4.00000 4.00000 0.00357
#> rejected_by_bonferroni
#> 0.00000Bonferroni controls the probability of any false positive and rejects far fewer hypotheses; BH controls the expected proportion among rejections and is correspondingly more powerful.
ggplot(data.frame(i = seq_len(m), p = pvals,
rejected = seq_len(m) <= k_hat),
aes(i, p)) +
geom_abline(slope = q / m, intercept = 0, color = "#3B7DD8", linewidth = 0.9) +
geom_hline(yintercept = q, linetype = "dashed", color = "grey40") +
geom_hline(yintercept = q / m, linetype = "dotted", color = "firebrick") +
geom_point(aes(color = rejected), size = 3) +
scale_color_manual(values = c(`FALSE` = "grey65", `TRUE` = "#D8433B"),
labels = c("not rejected", "rejected by BH"), name = NULL) +
annotate("text", x = m, y = q * 1.25, hjust = 1, size = 3.1, color = "grey35",
label = "uncorrected 0.05") +
annotate("text", x = m, y = q / m * 1.9, hjust = 1, size = 3.1,
color = "firebrick", label = "Bonferroni") +
labs(title = "Benjamini-Hochberg as a step-up rule",
subtitle = sprintf("Blue line: threshold q*i/m. The largest crossing index is k = %d; everything up to it is rejected", k_hat),
x = "Rank i", y = "p-value") +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(x = seq_len(m), y = pvals, type = "scatter", mode = "markers",
marker = list(size = 12), name = "p-values") |>
add_lines(x = c(0, m), y = c(0, q), name = "BH threshold q*i/m") |>
add_lines(x = c(0, m), y = rep(q, 2), name = "Uncorrected 0.05",
line = list(dash = "dash")) |>
add_lines(x = c(0, m), y = rep(q / m, 2), name = "Bonferroni",
line = list(dash = "dot")) |>
layout(title = "Benjamini-Hochberg FDR procedure",
xaxis = list(title = "Rank i"), yaxis = list(title = "p-value"))# The adjusted p-value is the step-up statistic, made monotone and capped at 1
adj_manual <- rev(cummin(rev(m / seq_len(m) * pvals)))
adj_manual <- pmin(adj_manual, 1)
data.frame(p = pvals,
manual = round(adj_manual, 5),
p.adjust_BH = round(p.adjust(pvals, "BH"), 5),
identical = abs(adj_manual - p.adjust(pvals, "BH")) < 1e-12)The rev(cummin(rev(.))) enforces monotonicity, adjusted
p-values must be non-decreasing in rank, and pmin(., 1)
caps them. Omitting either produces values that disagree with
p.adjust() and can exceed 1.
19 Stability selection
BH needs p-values. Most selectors, LASSO, trees, RFE, do not produce them. Stability selection (Meinshausen & Bühlmann, 2010) wraps error control around any selector by asking a different question: how often is this feature selected across perturbations of the data?
Algorithm. For \(b=1,\dots,B\): draw a subsample of size \(\lfloor n/2\rfloor\), run the selector, record the selected set \(\hat S_b\). The selection probability is
\[\hat\Pi_j=\frac1B\sum_{b=1}^{B}\mathbb 1\{j\in\hat S_b\},\]
and the stable set is \(\hat S^{\text{stable}}=\{j:\hat\Pi_j\ge\pi_{\text{thr}}\}\).
The bound. If the selector chooses at most \(q\) variables per subsample and the selection procedure is not worse than random exchangeable, then the expected number of falsely selected variables satisfies \[\boxed{\;\mathbb E[V]\ \le\ \frac{1}{2\pi_{\text{thr}}-1}\cdot\frac{q^2}{p}\;}\] for \(\pi_{\text{thr}}\in(0.5,1)\).
Two things follow. The threshold must exceed 0.5, or the denominator is non-positive and there is no bound. And the guarantee is on the expected count \(\mathbb E[V]\), not on the FDR, to convert, divide by the number selected.
stability_selection <- function(X, y, B = 100, q_sel = 10, seed = 91) {
set.seed(seed)
n <- nrow(X); p <- ncol(X)
counts <- numeric(p); names(counts) <- colnames(X)
for (b in seq_len(B)) {
idx <- sample(n, floor(n / 2))
f <- glmnet(X[idx, , drop = FALSE], y[idx], alpha = 1,
dfmax = q_sel, standardize = TRUE)
# Take the LAST lambda on the path with at most q_sel active variables
ok <- which(f$df > 0 & f$df <= q_sel)
if (!length(ok)) next
sel <- which(as.vector(f$beta[, max(ok)]) != 0)
counts[sel] <- counts[sel] + 1
}
counts / B
}
set.seed(93)
n_st <- 200; p_st <- 100
X_st <- matrix(rnorm(n_st * p_st), n_st, p_st,
dimnames = list(NULL, paste0("V", 1:p_st)))
true_st <- 1:6
b_st <- numeric(p_st); b_st[true_st] <- 2
y_st <- as.vector(X_st %*% b_st) + rnorm(n_st)
q_sel <- 10
pi_hat <- stability_selection(X_st, y_st, B = 100, q_sel = q_sel)
pi_thr <- 0.7
stable <- names(which(pi_hat >= pi_thr))
c(threshold = pi_thr,
q_per_subsample = q_sel,
expected_false_selections_bound = round(q_sel^2 / (p_st * (2 * pi_thr - 1)), 3),
n_selected = length(stable),
n_true_selected = sum(as.integer(sub("V", "", stable)) %in% true_st),
n_false_selected = sum(!as.integer(sub("V", "", stable)) %in% true_st))#> threshold q_per_subsample
#> 0.7 10.0
#> expected_false_selections_bound n_selected
#> 2.5 6.0
#> n_true_selected n_false_selected
#> 6.0 0.0#> V1 V2 V3 V4 V5 V6 V95 V70 V16 V31
#> 1.00 1.00 1.00 1.00 1.00 1.00 0.22 0.20 0.17 0.17data.frame(feature = names(pi_hat), pi = as.numeric(pi_hat),
truth = ifelse(seq_len(p_st) %in% true_st, "true signal", "null")) |>
arrange(desc(pi)) |> mutate(rank = row_number()) |>
ggplot(aes(rank, pi, color = truth)) +
geom_hline(yintercept = pi_thr, linetype = "dashed", color = "grey35") +
geom_hline(yintercept = 0.5, linetype = "dotted", color = "firebrick") +
geom_point(size = 2, alpha = 0.85) +
scale_color_manual(values = c(null = "grey65", `true signal` = "#D8433B")) +
labs(title = "Selection probabilities under subsampling",
subtitle = "Dashed: the chosen threshold. Dotted red: 0.5, below which the bound does not apply",
x = "Feature, ranked by selection probability",
y = expression(hat(Pi)[j]), color = NULL) +
theme_dspa()
The true signals separate cleanly from the nulls, and the bound tells you how many false selections to expect. Choosing \(\pi_{\text{thr}}\) is choosing a point on that tradeoff, higher threshold, fewer false selections, less power.
Selection probability is a function of both the feature and the regularization strength, which is the surface stability selection thresholds:
lam_grid_st <- exp(seq(log(0.5), log(0.03), length.out = 25))
B_st <- 60
set.seed(95)
Pi_mat <- matrix(0, p_st, length(lam_grid_st))
for (b in seq_len(B_st)) {
idx <- sample(n_st, floor(n_st / 2))
f <- glmnet(X_st[idx, ], y_st[idx], alpha = 1, lambda = lam_grid_st,
standardize = TRUE)
Pi_mat <- Pi_mat + (as.matrix(f$beta) != 0)
}
Pi_mat <- Pi_mat / B_st
ord_st <- order(rowMeans(Pi_mat), decreasing = TRUE)[1:30]
plot_ly(x = log(lam_grid_st), y = seq_along(ord_st),
z = Pi_mat[ord_st, ], type = "surface",
colorscale = "Viridis",
colorbar = list(title = "Selection\nprobability")) |>
layout(title = "Stability paths: selection probability over lambda and feature",
scene = list(xaxis = list(title = "log(lambda)"),
yaxis = list(title = "Feature (ranked)"),
zaxis = list(title = "Selection probability")))The six true signals form a plateau at probability 1 across a wide range of \(\lambda\); the nulls form a low, ragged foreground that rises only as \(\lambda\) becomes small enough to admit everything. Stability selection takes the maximum over \(\lambda\), which is why it is far less sensitive to the choice of \(\lambda\) than a single LASSO fit.
20 Knockoffs
Stability selection bounds the expected count of false selections. Knockoffs (Barber & Candès, 2015; Candès et al., 2018) give exact finite-sample FDR control, for any sample size, with no asymptotics.
20.1 Construction
Build a matrix \(\tilde X\) of decoy variables satisfying two properties:
\[ \begin{aligned} \textbf{(1) Exchangeability: }&\quad \big(X,\tilde X\big)_{\mathrm{swap}(S)}\ \overset{d}{=}\ \big(X,\tilde X\big)\quad\text{for every }S\subseteq\{1,\dots,p\}\\ \textbf{(2) Conditional independence: }&\quad \tilde X\ \perp\ Y\ \big|\ X \end{aligned} \]
Property (2) makes the decoys guaranteed nulls. Property (1) makes them indistinguishable from the real variables as far as the null features are concerned, so any importance statistic treats a null \(X_j\) and its decoy symmetrically, which is what converts a ranking into an error guarantee.
| Construction | Requires | Assumes about \(Y\mid X\) | Function |
|---|---|---|---|
| Fixed-X | \(n\ge2p\); Gaussian noise | Homoscedastic linear model | create.fixed() |
| Model-X, Gaussian | The law of \(X\) known | Nothing | create.gaussian() |
| Model-X, second-order | \(\mu,\Sigma\) estimated from data | Nothing | create.second_order() |
Model-X knockoffs move the assumption from \(Y\mid X\) to \(X\). That is usually the better trade in biomedical work: the design is often observational with many samples available to estimate \(\Sigma\), while the outcome model is exactly what is unknown.
Common misconception: “knockoffs are assumption-free.” They place no assumption on \(Y\mid X\), which is the headline result, and a genuine advance. But model-X knockoffs assume the joint law of \(X\) is known, and the FDR guarantee follows from exchangeability of \((X,\tilde X)\) under swaps, which holds only if the decoys are generated from the right distribution.
Supplying a covariance that does not describe the data, for instance \(\Sigma=I\) when the features are strongly correlated, produces decoys that are not exchangeable with the real variables, and the guarantee silently evaporates (§11.17.5 measures exactly this).
create.second_order()estimates \(\mu\) and \(\Sigma\) from the data and is the right default when they are unknown.
20.2 The \(W_j\) statistic
For each feature compute an importance for the real variable and its decoy, then combine them antisymmetrically.
For stat.glmnet_lambdasmax, let
\[Z_j=\sup\big\{\lambda:\ \hat\beta_j(\lambda)\ne0\big\},\qquad \tilde Z_j=\sup\big\{\lambda:\ \hat{\tilde\beta}_j(\lambda)\ne0\big\}\]
be the largest \(\lambda\) at which each first enters the LASSO path, a variable that enters early (at large \(\lambda\)) is important. Then
\[\boxed{\;W_j=\max\big(Z_j,\tilde Z_j\big)\cdot\operatorname{sign}\big(Z_j-\tilde Z_j\big)\;}\]
For stat.glmnet_coefdiff, \(W_j=|\hat\beta_j|-|\hat{\tilde\beta}_j|\)
at a cross-validated \(\lambda\).
The antisymmetry is what matters. Swapping \(X_j\) with \(\tilde X_j\) must flip the sign of \(W_j\). Combined with exchangeability, this makes the sign of \(W_j\) for a null feature a fair coin flip, so negative \(W_j\)’s estimate the number of false positives among the positive ones.
20.3 The threshold
\[\boxed{\;\tau=\min\left\{t>0:\ \frac{\#\{j:W_j\le-t\}}{\max\big(1,\#\{j:W_j\ge t\}\big)}\ \le\ q\right\}\;}\]
and select \(\hat S=\{j:W_j\ge\tau\}\). The numerator estimates the number of false discoveries at threshold \(t\) using the negative \(W_j\)’s as a proxy — that is the whole idea.
Knockoff+ adds 1 to the numerator:
\[\tau_+=\min\left\{t>0:\ \frac{1+\#\{j:W_j\le-t\}}{\max\big(1,\#\{j:W_j\ge t\}\big)}\le q\right\},\]
which upgrades the guarantee from modified FDR to exact FDR control at any finite \(n\), at a small cost in power.
library(knockoff)
set.seed(101)
n_k <- 1000; p_k <- 300; k_k <- 30; amplitude <- 3.5
X_k <- matrix(rnorm(n_k * p_k), n_k, p_k)
nonzero_k <- sample(p_k, k_k)
beta_k <- amplitude * (seq_len(p_k) %in% nonzero_k) / sqrt(n_k)
y_k <- as.vector(X_k %*% beta_k) + rnorm(n_k)
fdp <- function(selected) sum(!selected %in% nonzero_k) / max(1, length(selected))
power <- function(selected) sum(selected %in% nonzero_k) / k_k
res_k <- knockoff.filter(X_k, y_k, fdr = TARGET_FDR,
knockoffs = create.second_order,
statistic = stat.glmnet_coefdiff)
c(target_fdr = TARGET_FDR,
n_selected = length(res_k$selected),
observed_fdp = round(fdp(res_k$selected), 4),
power = round(power(res_k$selected), 4))#> target_fdr n_selected observed_fdp power
#> 0.1000 25.0000 0.0800 0.7667Wj <- res_k$statistic
data.frame(W = Wj, truth = ifelse(seq_len(p_k) %in% nonzero_k,
"true signal", "null")) |>
ggplot(aes(W, fill = truth)) +
geom_histogram(bins = 60, alpha = 0.8, position = "identity") +
geom_vline(xintercept = res_k$threshold, linetype = "dashed",
color = "black", linewidth = 0.9) +
geom_vline(xintercept = -res_k$threshold, linetype = "dotted",
color = "grey40") +
scale_fill_manual(values = c(null = "grey70", `true signal` = "#D8433B")) +
labs(title = "Knockoff statistics and the selection threshold",
subtitle = "Nulls are symmetric about zero by construction; the negatives estimate the false positives among the positives",
x = expression(W[j]), y = "Count", fill = NULL) +
theme_dspa()
c(nulls_positive = sum(Wj > 0 & !seq_len(p_k) %in% nonzero_k),
nulls_negative = sum(Wj < 0 & !seq_len(p_k) %in% nonzero_k),
symmetry_check = "these two counts should be similar")#> nulls_positive nulls_negative
#> "40" "32"
#> symmetry_check
#> "these two counts should be similar"The null \(W_j\)’s are symmetric about zero, the coin-flip property, while the true signals pile up on the positive side. Counting the negatives is therefore a valid estimate of how many nulls have crossed to the positive side.
# A custom antisymmetric statistic: the marginal-correlation difference
stat_marginal_diff <- function(X, X_k, y) {
abs(as.vector(crossprod(X, y))) - abs(as.vector(crossprod(X_k, y)))
}
set.seed(103)
res_custom <- knockoff.filter(X_k, y_k, fdr = TARGET_FDR,
knockoffs = create.second_order,
statistic = stat_marginal_diff)
data.frame(
statistic = c("glmnet coefficient difference", "LASSO signed max",
"marginal correlation difference"),
n_selected = c(length(res_k$selected),
length(knockoff.filter(X_k, y_k, fdr = TARGET_FDR,
knockoffs = create.second_order,
statistic = stat.glmnet_lambdasmax)$selected),
length(res_custom$selected)),
fdp = round(c(fdp(res_k$selected),
fdp(knockoff.filter(X_k, y_k, fdr = TARGET_FDR,
knockoffs = create.second_order,
statistic = stat.glmnet_lambdasmax)$selected),
fdp(res_custom$selected)), 4))Any antisymmetric statistic is valid; they differ only in power. The LASSO-based statistics use the multivariate fit and are stronger than a marginal one, which is blind to the same interactions as a univariate filter (§11.4).
20.4 Verifying the guarantee
FDR control is a statement about an expectation, so a single run cannot demonstrate it.
set.seed(105)
replicate_ko <- function(q, amp, reps = KO_REPS, n = 400, p = 150, k = 15) {
out <- vapply(seq_len(reps), function(r) {
X <- matrix(rnorm(n * p), n, p)
nz <- sample(p, k)
b <- amp * (seq_len(p) %in% nz) / sqrt(n)
y <- as.vector(X %*% b) + rnorm(n)
s <- tryCatch(knockoff.filter(X, y, fdr = q,
knockoffs = create.second_order,
statistic = stat.glmnet_coefdiff)$selected,
error = function(e) integer(0))
c(fdp = sum(!s %in% nz) / max(1, length(s)),
pow = sum(s %in% nz) / k)
}, numeric(2))
c(q = q, amplitude = amp,
mean_FDP = mean(out["fdp", ]), mean_power = mean(out["pow", ]))
}
ko_tab <- as.data.frame(do.call(rbind, lapply(
c(0.05, 0.10, 0.20), replicate_ko, amp = 4)))
ko_tab |> mutate(across(everything(), \(z) round(z, 4)))The mean FDP stays at or below the target \(q\) in every row, which is the guarantee, and note that it is typically below rather than at \(q\), since knockoffs are conservative.
q_grid <- c(0.05, 0.10, 0.15, 0.20, 0.30)
amp_grid <- c(2, 3, 4, 5, 6)
ko_grid <- expand.grid(q = q_grid, amp = amp_grid)
ko_res <- t(vapply(seq_len(nrow(ko_grid)), function(i)
replicate_ko(ko_grid$q[i], ko_grid$amp[i], reps = max(8, KO_REPS %/% 3)),
numeric(4)))
Z_fdp <- matrix(ko_res[, 3], length(q_grid), length(amp_grid))
Z_pow <- matrix(ko_res[, 4], length(q_grid), length(amp_grid))
plot_ly() |>
add_surface(x = amp_grid, y = q_grid, z = Z_pow, opacity = 0.9,
showscale = FALSE, colorscale = "Viridis", name = "Power") |>
add_surface(x = amp_grid, y = q_grid, z = Z_fdp, opacity = 0.9,
showscale = FALSE, colorscale = "Reds", name = "FDP") |>
layout(title = "Knockoff power (upper, viridis) and realized FDP (lower, red)",
scene = list(xaxis = list(title = "Signal amplitude"),
yaxis = list(title = "Target FDR q"),
zaxis = list(title = "Proportion")))Two things to read. Power rises steeply with signal amplitude and gently with \(q\), accepting more false discoveries buys some, but not much, additional detection. Realized FDP tracks the target \(q\) along the \(q\) axis and is essentially flat in amplitude, which is exactly the guarantee: the error rate is controlled regardless of how strong the signal happens to be.
# What breaks when the assumed covariance is wrong
set.seed(107)
n_b <- 400; p_b <- 100
Sig_b <- 0.8^abs(outer(1:p_b, 1:p_b, "-")) # strongly correlated design
X_b <- MASS::mvrnorm(n_b, rep(0, p_b), Sig_b)
nz_b <- sample(p_b, 10)
b_b <- 4 * (seq_len(p_b) %in% nz_b) / sqrt(n_b)
y_b <- as.vector(X_b %*% b_b) + rnorm(n_b)
correct <- knockoff.filter(X_b, y_b, fdr = TARGET_FDR,
knockoffs = create.second_order,
statistic = stat.glmnet_coefdiff)
wrong <- knockoff.filter(X_b, y_b, fdr = TARGET_FDR,
knockoffs = function(X)
create.gaussian(X, rep(0, ncol(X)), diag(ncol(X))),
statistic = stat.glmnet_coefdiff)
data.frame(
construction = c("second-order (Sigma estimated)", "Gaussian with Sigma = I"),
n_selected = c(length(correct$selected), length(wrong$selected)),
realized_fdp = round(c(sum(!correct$selected %in% nz_b) / max(1, length(correct$selected)),
sum(!wrong$selected %in% nz_b) / max(1, length(wrong$selected))), 4),
target = TARGET_FDR)The features here are correlated at \(\rho=0.8\) between neighbours. Asserting \(\Sigma=I\) constructs decoys that are not exchangeable with the real variables, and the realized FDP is no longer controlled. The guarantee is only as good as the model for \(X\).
21 PART V: WRAPPERS AND EMBEDDED METHODS
22 Case study: amyotrophic lateral sclerosis
The ALS
case study comes from a large clinical trial in amyotrophic
lateral sclerosis, a rare neurodegenerative disorder. The training
set has 2,223 observations and 131 numeric variables. The outcome,
ALSFRS_slope, measures clinical decline
over a year, and many predictors are max, min,
and median summaries of the same underlying
clinical measurement, so the design is heavily correlated by
construction.
als_raw <- if (ONLINE) {
dspa_read("https://umich.instructure.com/files/1789624/download?download_frd=1",
"ALS_train.csv")
} else NULL
if (is.null(als_raw)) {
set.seed(111); n_a <- 600; p_a <- 60
als_raw <- as.data.frame(matrix(rnorm(n_a * p_a), n_a, p_a))
names(als_raw) <- paste0("f", seq_len(p_a))
als_raw$ID <- seq_len(n_a)
als_raw$ALSFRS_slope <- rowSums(als_raw[, 1:6]) * 0.3 + rnorm(n_a)
message("Note: ALS file unavailable; using a synthetic substitute.")
}
c(rows = nrow(als_raw), columns = ncol(als_raw))#> rows columns
#> 2223 101stopifnot("ALSFRS_slope" %in% names(als_raw))
# The feature set is defined ONCE, by name, and reused everywhere below
OUTCOME <- "ALSFRS_slope"
DROP_IDS <- intersect(c("ID", "SubjectID"), names(als_raw))
FEATURES <- setdiff(names(als_raw), c(OUTCOME, DROP_IDS))
c(outcome = OUTCOME, dropped = paste(DROP_IDS, collapse = ", "),
n_features = length(FEATURES))#> outcome dropped n_features
#> "ALSFRS_slope" "ID, SubjectID" "98"als <- als_raw[, c(OUTCOME, FEATURES)]
als <- als[complete.cases(als), ]
als <- als[, c(TRUE, vapply(als[, -1, drop = FALSE], \(v)
is.numeric(v) && var(v) > 0, logical(1)))]
FEATURES <- setdiff(names(als), OUTCOME)
c(rows_complete = nrow(als), features_retained = length(FEATURES))#> rows_complete features_retained
#> 2223 98set.seed(1234)
sp_als <- rsample::initial_split(als, prop = 0.75)
als_train <- rsample::training(sp_als)
als_test <- rsample::testing(sp_als)
c(train = nrow(als_train), test = nrow(als_test))#> train test
#> 1667 556cm_als <- cor(als_train[, FEATURES])
offdiag <- cm_als[upper.tri(cm_als)]
c(features = length(FEATURES),
pairs = length(offdiag),
median_abs_correlation = round(median(abs(offdiag)), 4),
pairs_above_0.9 = sum(abs(offdiag) > 0.9),
pairs_above_0.99 = sum(abs(offdiag) > 0.99))#> features pairs median_abs_correlation
#> 98.0000 4753.0000 0.0549
#> pairs_above_0.9 pairs_above_0.99
#> 6.0000 0.0000ggplot(data.frame(r = offdiag), aes(abs(r))) +
geom_histogram(bins = 60, fill = "steelblue", color = "white") +
geom_vline(xintercept = 0.9, linetype = "dashed", color = "firebrick") +
labs(title = "Pairwise absolute correlations among ALS predictors",
subtitle = "The mass above 0.9 is what makes plain LASSO selection unstable here",
x = "|correlation|", y = "Pairs") +
theme_dspa()
That correlation structure is the whole difficulty. Where predictors are near-duplicates, LASSO picks one arbitrarily (§11.8.1), the irrepresentable condition (§11.12) is likely violated, and any single selection run is unstable.
23 Recursive feature elimination
RFE is a backward wrapper: fit a model on all features, rank them by importance, drop the weakest, refit, repeat. The subset size is chosen by resampling.
RFE must re-rank inside every resampling fold. Ranking once on the full data and then cross-validating only the subset size reintroduces the selection bias of §11.3.
caret::rfe()does the right thing by default, the importance ranking is recomputed within each fold — which is why its accuracy estimates are honest and why any selector compared against it must be run the same way.
library(caret)
library(ranger)
set.seed(121)
ctrl_rfe <- rfeControl(
functions = rfFuncs,
method = "cv",
number = 5,
verbose = TRUE,
allowParallel = FALSE
)
# ctrl_rfe <- rfeControl(functions = rfFuncs, method = "cv", number = 5,
# verbose = FALSE)
sizes_rfe <- c(5, 10, 20, 40)
# rfe_fit <- rfe(x = als_train[, FEATURES], y = als_train[[OUTCOME]],
# sizes = sizes_rfe, rfeControl = ctrl_rfe)
rfe_fit <- rfe(x = als_train[, FEATURES], y = als_train[[OUTCOME]],
sizes = sizes_rfe, rfeControl = ctrl_rfe)#> +(rfe) fit Fold1 size: 98
#> -(rfe) fit Fold1 size: 98
#> +(rfe) imp Fold1
#> -(rfe) imp Fold1
#> +(rfe) fit Fold1 size: 40
#> -(rfe) fit Fold1 size: 40
#> +(rfe) fit Fold1 size: 20
#> -(rfe) fit Fold1 size: 20
#> +(rfe) fit Fold1 size: 10
#> -(rfe) fit Fold1 size: 10
#> +(rfe) fit Fold1 size: 5
#> -(rfe) fit Fold1 size: 5
#> +(rfe) fit Fold2 size: 98
#> -(rfe) fit Fold2 size: 98
#> +(rfe) imp Fold2
#> -(rfe) imp Fold2
#> +(rfe) fit Fold2 size: 40
#> -(rfe) fit Fold2 size: 40
#> +(rfe) fit Fold2 size: 20
#> -(rfe) fit Fold2 size: 20
#> +(rfe) fit Fold2 size: 10
#> -(rfe) fit Fold2 size: 10
#> +(rfe) fit Fold2 size: 5
#> -(rfe) fit Fold2 size: 5
#> +(rfe) fit Fold3 size: 98
#> -(rfe) fit Fold3 size: 98
#> +(rfe) imp Fold3
#> -(rfe) imp Fold3
#> +(rfe) fit Fold3 size: 40
#> -(rfe) fit Fold3 size: 40
#> +(rfe) fit Fold3 size: 20
#> -(rfe) fit Fold3 size: 20
#> +(rfe) fit Fold3 size: 10
#> -(rfe) fit Fold3 size: 10
#> +(rfe) fit Fold3 size: 5
#> -(rfe) fit Fold3 size: 5
#> +(rfe) fit Fold4 size: 98
#> -(rfe) fit Fold4 size: 98
#> +(rfe) imp Fold4
#> -(rfe) imp Fold4
#> +(rfe) fit Fold4 size: 40
#> -(rfe) fit Fold4 size: 40
#> +(rfe) fit Fold4 size: 20
#> -(rfe) fit Fold4 size: 20
#> +(rfe) fit Fold4 size: 10
#> -(rfe) fit Fold4 size: 10
#> +(rfe) fit Fold4 size: 5
#> -(rfe) fit Fold4 size: 5
#> +(rfe) fit Fold5 size: 98
#> -(rfe) fit Fold5 size: 98
#> +(rfe) imp Fold5
#> -(rfe) imp Fold5
#> +(rfe) fit Fold5 size: 40
#> -(rfe) fit Fold5 size: 40
#> +(rfe) fit Fold5 size: 20
#> -(rfe) fit Fold5 size: 20
#> +(rfe) fit Fold5 size: 10
#> -(rfe) fit Fold5 size: 10
#> +(rfe) fit Fold5 size: 5
#> -(rfe) fit Fold5 size: 5#>
#> Recursive feature selection
#>
#> Outer resampling method: Cross-Validated (5 fold)
#>
#> Resampling performance over subset size:
#>
#> Variables RMSE Rsquared MAE RMSESD RsquaredSD MAESD Selected
#> 5 0.346 0.689 0.248 0.0183 0.0326 0.00608
#> 10 0.340 0.702 0.243 0.0172 0.0325 0.00546
#> 20 0.337 0.706 0.243 0.0191 0.0339 0.00805 *
#> 40 0.339 0.704 0.244 0.0209 0.0355 0.00926
#> 98 0.341 0.701 0.245 0.0225 0.0357 0.00948
#>
#> The top 5 variables (out of 20):
#> ALSFRS_Total_range, trunk_range, hands_range, mouth_range, ALSFRS_Total_minrfe_selected <- predictors(rfe_fit)
c(optimal_size = rfe_fit$bestSubset, n_selected = length(rfe_selected))#> optimal_size n_selected
#> 20 20ggplot(rfe_fit$results, aes(Variables, RMSE)) +
geom_errorbar(aes(ymin = RMSE - RMSESD / sqrt(5),
ymax = RMSE + RMSESD / sqrt(5)),
width = 0.6, color = "grey60") +
geom_line(linewidth = 0.9, color = "steelblue") +
geom_point(size = 2.4) +
geom_point(data = rfe_fit$results[which.min(rfe_fit$results$RMSE), ],
color = "firebrick", size = 4) +
scale_x_log10(breaks = rfe_fit$results$Variables) +
labs(title = "Recursive feature elimination: cross-validated RMSE by subset size",
subtitle = "Ranking is recomputed inside every fold, so these estimates are not selection-biased",
x = "Number of features (log scale)", y = "CV RMSE") +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(rfe_fit$results, x = ~Variables, y = ~RMSE, type = "scatter",
mode = "lines+markers",
error_y = ~list(array = RMSESD / sqrt(5))) |>
layout(title = "RFE: CV RMSE by subset size",
xaxis = list(title = "Number of features", type = "log"),
yaxis = list(title = "RMSE"))24 Boruta
Boruta is an all-relevant selector built around random forests. Its logic is the same as the knockoff idea, arrived at independently.
Algorithm. For each run: duplicate every feature and permute each copy to produce shadow features that are guaranteed nulls; fit a random forest on the extended data; record each real feature’s importance and the maximum shadow importance; count a real feature as a “hit” if it beats that maximum. After \(M\) runs, test each feature’s hit count against \(\mathrm{Binomial}(M,0.5)\) with a multiplicity correction, and label it Confirmed, Rejected, or Tentative.
The decision rule is a binomial test against the maximum shadow importance, not against the mean or against a fixed cutoff. Using the maximum is what makes the comparison conservative and is the direct analogue of the knockoff threshold.
Common misconception: “Boruta and LASSO disagree, so one of them is wrong.” They are solving different problems. Boruta is all-relevant: it aims to find every feature carrying information about the outcome, including redundant ones. The LASSO is minimal-optimal: it seeks the smallest set sufficient for prediction, and will drop a perfect duplicate because keeping it adds nothing.
On data with near-duplicate predictors the two must diverge, and both are behaving correctly. Which you want depends on the question: all-relevant for “what is involved in this biology,” minimal-optimal for “what do I need to predict.” Overlap between them measures redundancy in the design, not the quality of either method.
library(Boruta)
set.seed(123)
bor <- Boruta(x = als_train[, FEATURES], y = als_train[[OUTCOME]],
maxRuns = BORUTA_RUNS, doTrace = 0)
bor#> Boruta performed 59 iterations in 29.4317 secs.
#> 26 attributes confirmed important: ALSFRS_Total_max,
#> ALSFRS_Total_median, ALSFRS_Total_min, ALSFRS_Total_range,
#> Creatinine_min and 21 more;
#> 58 attributes confirmed unimportant: Age_mean, Albumin_max,
#> Albumin_median, Albumin_min, Albumin_range and 53 more;
#> 14 tentative attributes left: ALT.SGPT._min, bp_systolic_range,
#> Chloride_range, Creatinine_max, Creatinine_median and 9 more;dec <- bor$finalDecision
c(confirmed = sum(dec == "Confirmed"),
tentative = sum(dec == "Tentative"),
rejected = sum(dec == "Rejected"),
total = length(dec))#> confirmed tentative rejected total
#> 26 14 58 98bor_fixed <- TentativeRoughFix(bor)
boruta_selected <- getSelectedAttributes(bor_fixed, withTentative = FALSE)
c(after_rough_fix = length(boruta_selected))#> after_rough_fix
#> 34imp_hist <- as.data.frame(bor$ImpHistory)
imp_long <- imp_hist |>
pivot_longer(everything(), names_to = "feature", values_to = "importance") |>
filter(is.finite(importance))
med <- imp_long |> summarise(m = median(importance), .by = feature) |>
arrange(desc(m))
top_show <- head(med$feature, 28)
imp_long |> filter(feature %in% top_show) |>
mutate(feature = factor(feature, levels = rev(top_show)),
kind = case_when(grepl("^shadow", feature) ~ "shadow",
feature %in% boruta_selected ~ "confirmed",
TRUE ~ "other")) |>
ggplot(aes(feature, importance, fill = kind)) +
geom_boxplot(outlier.size = 0.4, linewidth = 0.3) +
coord_flip() +
scale_fill_manual(values = c(shadow = "#D8433B", confirmed = "#3B7DD8",
other = "grey75")) +
labs(title = "Boruta importance history: real features against shadows",
subtitle = "Red boxes are permuted shadow features. A real feature is confirmed by beating the shadow MAXIMUM",
x = NULL, y = "Importance (Z-score of mean decrease in accuracy)",
fill = NULL) +
theme_dspa(8)
# --- Interactive equivalent ------------------------------------------------
plot_ly(imp_long |> filter(feature %in% top_show),
x = ~feature, y = ~importance, color = ~feature, type = "box") |>
layout(title = sprintf("Boruta importance across %d runs", BORUTA_RUNS),
xaxis = list(title = "Feature", categoryorder = "median descending"),
yaxis = list(title = "Importance"), showlegend = FALSE)The shadow boxes sit low and tight; confirmed features clear them consistently. That separation is the evidence, not the absolute importance values, which have no scale of their own.
25 Minimum redundancy, maximum relevance
mRMR is an information-theoretic multivariate filter. It selects features that are individually informative about the outcome and mutually non-redundant.
Let \(I(f;y)\) be the mutual information between feature \(f\) and the outcome,
\[I(f;y)=\iint p(f,y)\log\frac{p(f,y)}{p(f)p(y)}\,df\,dy,\]
and \(I(f_i;f_j)\) the mutual information between two features. The greedy criterion adds, at each step,
\[\boxed{\;f_{\text{next}}=\arg\max_{f\notin S}\left[\underbrace{I(f;y)}_{\text{relevance}}-\underbrace{\frac{1}{|S|}\sum_{f_i\in S}I(f;f_i)}_{\text{redundancy}}\right]\;}\]
Because it explicitly penalizes redundancy, mRMR is well matched to a design like the ALS data where many predictors are near-duplicates, and it costs \(O(p|S|)\) mutual-information evaluations rather than the \(O(p)\) model refits an RFE-style wrapper needs.
MRMR_OK <- requireNamespace("mRMRe", quietly = TRUE)
if (MRMR_OK) {
suppressMessages(library(mRMRe))
mr_df <- als_train[, c(OUTCOME, FEATURES)]
mr_df[] <- lapply(mr_df, as.numeric)
dd <- mRMR.data(data = mr_df)
n_keep <- 20
mr <- mRMR.classic(data = dd, target_indices = 1, feature_count = n_keep)
mrmr_selected <- FEATURES[unlist(solutions(mr)) - 1]
c(n_selected = length(mrmr_selected))
head(mrmr_selected, 10)
} else {
message("mRMRe unavailable; the praznik package provides mRMR, JMI, CMIM and DISR.")
mrmr_selected <- character(0)
}#> [1] "ALSFRS_Total_range" "Potassium_range" "hands_range"
#> [4] "mouth_range" "trunk_range" "leg_range"
#> [7] "ALSFRS_Total_min" "onset_delta_mean" "respiratory_range"
#> [10] "pulse_range"26 LASSO and elastic net on the ALS data
Xa <- as.matrix(als_train[, FEATURES])
ya <- als_train[[OUTCOME]]
Xa_test <- as.matrix(als_test[, FEATURES])
ya_test <- als_test[[OUTCOME]]
set.seed(131)
cv_lasso_a <- cv.glmnet(Xa, ya, alpha = 1, nfolds = 10)
cv_enet_a <- cv.glmnet(Xa, ya, alpha = 0.5, nfolds = 10)
plot_cv_glmnet(cv_lasso_a, "ALS, LASSO")
lasso_selected <- FEATURES[which(as.vector(coef(cv_lasso_a, s = "lambda.1se"))[-1] != 0)]
enet_selected <- FEATURES[which(as.vector(coef(cv_enet_a, s = "lambda.1se"))[-1] != 0)]
c(lasso_selected = length(lasso_selected),
elastic_net_selected = length(enet_selected))#> lasso_selected elastic_net_selected
#> 1 6set.seed(133)
pi_als <- stability_selection(Xa, ya, B = 60, q_sel = 25, seed = 133)
stable_als <- names(sort(pi_als[pi_als >= 0.7], decreasing = TRUE))
c(stable_features = length(stable_als),
expected_false_bound = round(25^2 / (length(FEATURES) * (2 * 0.7 - 1)), 2))#> stable_features expected_false_bound
#> 5.00 15.94#> ALSFRS_Total_median ALSFRS_Total_range Potassium_range respiratory_max
#> 1.000000 1.000000 0.833333 0.733333
#> trunk_min ALSFRS_Total_min onset_delta_mean hands_range
#> 0.733333 0.666667 0.650000 0.633333
#> Calcium_range mouth_range
#> 0.566667 0.53333327 A like-for-like comparison
These selectors answer different questions, so agreement is not the goal. LASSO and elastic net are minimal-optimal: smallest set sufficient for prediction. Boruta is all-relevant: every feature carrying information. RFE optimizes a predictive criterion at a chosen subset size. Stability selection controls the expected number of false selections.
On a design with correlations above 0.99, minimal-optimal methods will drop members of a duplicated group that all-relevant methods keep. That is correct behavior from both, and their overlap measures redundancy in the data rather than the quality of either method.
sel_sets <- list(LASSO = lasso_selected, ElasticNet = enet_selected,
RFE = rfe_selected, Boruta = boruta_selected,
Stability = stable_als)
if (length(mrmr_selected)) sel_sets$mRMR <- mrmr_selected
data.frame(method = names(sel_sets),
n_selected = vapply(sel_sets, length, integer(1)),
row.names = NULL)ov <- outer(seq_along(sel_sets), seq_along(sel_sets),
Vectorize(function(i, j) length(intersect(sel_sets[[i]], sel_sets[[j]]))))
dimnames(ov) <- list(names(sel_sets), names(sel_sets))
ov#> LASSO ElasticNet RFE Boruta Stability mRMR
#> LASSO 1 1 1 1 1 1
#> ElasticNet 1 6 6 6 1 6
#> RFE 1 6 20 20 3 12
#> Boruta 1 6 20 34 3 13
#> Stability 1 1 3 3 5 3
#> mRMR 1 6 12 13 3 20as.data.frame(as.table(ov)) |>
setNames(c("m1", "m2", "n")) |>
ggplot(aes(m1, m2, fill = n)) +
geom_tile(color = "white") +
geom_text(aes(label = n), size = 3.4) +
scale_fill_viridis_c(option = "mako", direction = -1, name = "Shared") +
coord_fixed() +
labs(title = "Features selected in common by each pair of methods",
subtitle = "Diagonal entries are set sizes. Low off-diagonal counts reflect different objectives, not disagreement about the data",
x = NULL, y = NULL) +
theme_dspa(10)
# The comparison that IS commensurable: held-out predictive performance
# of an OLS model refitted on each selected set.
eval_set <- function(feats, label) {
if (!length(feats)) return(NULL)
f <- as.formula(paste(OUTCOME, "~", paste(sprintf("`%s`", feats), collapse = " + ")))
m <- lm(f, data = als_train)
p <- predict(m, newdata = als_test) # newdata, NOT newx
data.frame(method = label, n_features = length(feats),
test_rmse = sqrt(mean((p - ya_test)^2)),
test_mae = mean(abs(p - ya_test)),
test_r2 = 1 - sum((p - ya_test)^2) / sum((ya_test - mean(ya_test))^2))
}
pred_lasso <- as.vector(predict(cv_lasso_a, newx = Xa_test, s = "lambda.1se"))
baseline <- data.frame(
method = c("Intercept only", "LASSO (shrunk coefficients)"),
n_features = c(0, length(lasso_selected)),
test_rmse = c(sqrt(mean((mean(ya) - ya_test)^2)),
sqrt(mean((pred_lasso - ya_test)^2))),
test_mae = c(mean(abs(mean(ya) - ya_test)), mean(abs(pred_lasso - ya_test))),
test_r2 = c(1 - sum((mean(ya) - ya_test)^2) / sum((ya_test - mean(ya_test))^2),
1 - sum((pred_lasso - ya_test)^2) / sum((ya_test - mean(ya_test))^2)))
comparison <- bind_rows(baseline,
bind_rows(lapply(names(sel_sets), \(nm) eval_set(sel_sets[[nm]], paste(nm, "(refit OLS)")))))
comparison |> mutate(across(where(is.numeric), \(z) round(z, 4))) |> arrange(test_rmse)Every row is evaluated on the same held-out set with
predict(..., newdata =), so the numbers are comparable.
Note the two LASSO rows: refitting OLS on the selected variables removes
the shrinkage, which usually increases variance and can
worsen held-out error, the “relaxed LASSO” tradeoff.
comparison |>
filter(n_features > 0) |>
ggplot(aes(n_features, test_rmse)) +
geom_hline(yintercept = comparison$test_rmse[1], linetype = "dashed",
color = "grey45") +
annotate("text", x = Inf, y = comparison$test_rmse[1], hjust = 1.05, vjust = -0.6,
size = 3.1, color = "grey35", label = "intercept-only baseline") +
geom_point(aes(color = method), size = 4) +
scale_color_viridis_d(option = "turbo", end = 0.9) +
labs(title = "Held-out RMSE against model size",
subtitle = "Lower and further left is better: comparable accuracy from fewer variables",
x = "Number of features selected", y = "Test RMSE", color = NULL) +
theme_dspa(10)
28 Complexity
\(n\) = observations, \(p\) = features, \(s\) = selected size, \(K\) = path length, \(B\) = resampling replicates, \(M\) = Boruta runs, \(T\) = trees per forest.
| Method | Cost | Error control | Sees interactions? |
|---|---|---|---|
| Univariate filter | \(O(np)\) | Per-test only | No |
| SIS | \(O(np)\) | Sure screening (asymptotic) | No |
| mRMR (greedy) | \(O(ps)\) MI evaluations | None | Pairwise |
| LASSO path | \(O(npK)\) with warm starts | None | Yes |
| LARS | \(O(p^3+np^2)\) | None | Yes |
| Elastic net | as LASSO | None | Yes |
| Best subset | \(O(2^p)\) | None | Yes |
| OMP | \(O(nps)\) | None | Yes |
| RFE (\(k\)-fold) | \(k\cdot O(p)\) model refits | None (CV estimate only) | Yes |
| Boruta | \(O(M\cdot T\cdot np\log n)\) | Binomial test, corrected | Yes |
| Stability selection | \(B\times\) selector cost | \(\mathbb E[V]\le\frac{q^2}{p(2\pi_{\mathrm{thr}}-1)}\) | Inherited |
| Knockoffs (SDP) | \(O(p^3)\) + one path | Exact finite-sample FDR | Yes |
| Knockoffs (ASDP) | \(O(p\,b^2)\) for block size \(b\) | Same | Yes |
Three consequences. The \(O(p^3)\) SDP is the knockoff
bottleneck, use method = "asdp" beyond a few
hundred features. Stability selection multiplies whatever it
wraps by \(B\), which is why
it is usually wrapped around the LASSO rather than around an RFE. And
only knockoffs give a finite-sample guarantee;
everything else is either asymptotic, resampling-based, or
uncontrolled.
29 PART VI: SYNTHESIS
30 Common pitfalls
| # | Pitfall | Consequence | Fix |
|---|---|---|---|
| 1 | Selecting features on all the data, then cross-validating the model | Near-perfect accuracy on pure noise | Re-run selection inside every training fold |
| 2 | Comparing a selector run on full data against one run inside folds | The two estimates are not commensurable | Run both under the same resampling scheme |
| 3 | predict(lm_fit, newx = ...) |
Argument silently ignored; returns training fits | newdata = |
| 4 | Comparing a training error against test errors | The regularized methods appear better than they are | Every arm evaluated on the same held-out set |
| 5 | Naming a mean squared error “RMS” | Off by a square root wherever it is quoted | Report RMSE explicitly, or name it MSE |
| 6 | Positional train/test splits on ordered data | Systematic, not random, held-out set | rsample::initial_split(strata = ) |
| 7 | Quoting a soft threshold without its fidelity convention | Off by a factor of 2 | \(\frac{1}{2n}\|\cdot\|^2\Rightarrow\) threshold \(n\lambda\) |
| 8 | Identifying best subset with OMP | One is NP-hard and exact, the other greedy | \(\ell_0\) exact vs. forward greedy |
| 9 | Writing the elastic net without \(\lambda\) | \(\alpha\) sets the mix; \(\lambda\) sets the weight | Both are required |
| 10 | One soft-threshold pass called “LASSO” | Correct only under orthonormal design | Iterate coordinate descent to convergence |
| 11 | Expecting LASSO to find the truth given enough data | Selection consistency needs the irrepresentable condition | Check it; use adaptive LASSO or stability selection |
| 12 | Reading LASSO output on correlated groups as “the important variable” | One of the group is picked arbitrarily | Elastic net for the grouping effect |
| 13 | Reporting naive p-values after selection | Anti-conservative; rejects on pure noise | Selective inference, or data splitting |
| 14 | Applying BH as a step-down rule | Misses rejections below the largest crossing index | Find \(\hat k\), reject all \(i\le\hat k\) |
| 15 | Adjusted p-values not made monotone or capped | Disagrees with p.adjust; values above 1 |
rev(cummin(rev(.))) then pmin(., 1) |
| 16 | Stability threshold at or below 0.5 | The bound’s denominator is non-positive | \(\pi_{\mathrm{thr}}\in(0.5,1)\) |
| 17 | Reading stability’s bound as an FDR | It bounds the expected count \(\mathbb E[V]\) | Divide by the number selected |
| 18 | Model-X knockoffs with an assumed \(\Sigma = I\) | Exchangeability fails; FDR control is void | create.second_order() |
| 19 | Judging FDR control from one knockoff run | The guarantee is on an expectation | Replicate and average the FDP |
| 20 | A \(W_j\) statistic that is not antisymmetric | Null signs are no longer coin flips | Swapping \(X_j,\tilde X_j\) must flip \(W_j\) |
| 21 | Expecting Boruta and LASSO to agree | All-relevant vs. minimal-optimal objectives | Compare on held-out prediction instead |
| 22 | Univariate screening for interaction effects | Both true predictors rank with the noise | Multivariate or model-based selection |
| 23 | attach() or rm(list = ls()) in an analysis
script |
Shadow copies; destroyed environments | with(), $, explicit
data = |
| 24 | Transcribing selection counts into prose | Desynchronizes on any change to seed or data | Compute inline with sprintf() |
31 Practice problems
31.1 Problem 1: How big is selection bias?
Measure the inflation as a function of how many features are screened in.
Solution
set.seed(201)
keep_grid <- c(2, 5, 10, 25, 50)
p1 <- do.call(rbind, lapply(keep_grid, function(k) {
reps <- replicate(12, { d <- simulate_null(n = 80, p = 2000)
c(biased_cv(d, k_keep = k), honest_cv(d, k_keep = k)) })
data.frame(k_kept = k, biased = mean(reps[1, ]), honest = mean(reps[2, ]))
}))
p1 |> mutate(inflation = biased - honest,
across(where(is.numeric), \(z) round(z, 4)))p1 |> pivot_longer(c(biased, honest), names_to = "protocol", values_to = "acc") |>
ggplot(aes(k_kept, acc, color = protocol)) +
geom_hline(yintercept = 0.5, linetype = "dashed", color = "grey40") +
geom_line(linewidth = 1) + geom_point(size = 2.4) +
scale_x_log10(breaks = keep_grid) +
scale_color_manual(values = c(biased = "#D8433B", honest = "#3B7DD8")) +
labs(title = "Selection bias against the number of screened-in features",
subtitle = "Data contain no signal. The honest protocol tracks 0.5 throughout",
x = "Features retained by the screen (log scale)", y = "CV accuracy",
color = NULL) +
theme_dspa()
31.2 Problem 2: Verify both closed forms
Confirm the ridge shrinkage factor and the LASSO soft threshold numerically under orthonormal design, for several \(\lambda\).
Solution
set.seed(203)
n2 <- 300; p2 <- 6
Q2 <- qr.Q(qr(matrix(rnorm(n2 * p2), n2, p2)))
y2 <- as.vector(Q2 %*% c(4, -3, 2, 1, 0, 0)) + rnorm(n2, sd = 0.4)
b2_ols <- as.vector(crossprod(Q2, y2))
do.call(rbind, lapply(c(0.002, 0.01, 0.05), function(l) {
r <- as.vector(coef(glmnet(Q2, y2, alpha = 0, lambda = l,
standardize = FALSE, intercept = FALSE)))[-1]
s <- as.vector(coef(glmnet(Q2, y2, alpha = 1, lambda = l,
standardize = FALSE, intercept = FALSE)))[-1]
data.frame(lambda = l,
ridge_max_err = max(abs(r - b2_ols / (1 + 2 * n2 * l))),
lasso_max_err = max(abs(s - soft_threshold(b2_ols, n2 * l))))
})) |> mutate(across(-lambda, \(z) signif(z, 3)))glmnet to numerical precision,
using the same \(\frac{1}{2n}\)
fidelity convention for each. Substituting a threshold derived
under a different scaling would fail here immediately.
31.3 Problem 3: Grouping effect across \(\alpha\)
Trace how the elastic net distributes weight within a correlated group as \(\alpha\) varies.
Solution
set.seed(205)
n3 <- 150; z3 <- rnorm(n3)
X3 <- cbind(sapply(1:4, \(i) z3 + rnorm(n3, sd = 0.08)),
matrix(rnorm(n3 * 8), n3, 8))
colnames(X3) <- c(paste0("g", 1:4), paste0("n", 1:8))
y3 <- 4 * z3 + rnorm(n3)
a_seq <- c(0.05, 0.2, 0.4, 0.6, 0.8, 1)
p3 <- do.call(rbind, lapply(a_seq, function(a) {
set.seed(207)
cv <- cv.glmnet(X3, y3, alpha = a, nfolds = 10)
cf <- as.vector(coef(cv, s = "lambda.min"))[-1]
data.frame(alpha = a, n_group_selected = sum(cf[1:4] != 0),
group_coef_sd = sd(cf[1:4]),
total_group_weight = sum(abs(cf[1:4])),
n_noise_selected = sum(cf[5:12] != 0))
}))
p3 |> mutate(across(where(is.numeric), \(z) round(z, 4)))ggplot(p3, aes(alpha, n_group_selected)) +
geom_line(linewidth = 1, color = "steelblue") + geom_point(size = 2.4) +
scale_y_continuous(breaks = 0:4, limits = c(0, 4)) +
labs(title = "Members of a correlated group retained, against alpha",
subtitle = "Four predictors correlated above 0.99, all genuinely associated with y",
x = expression(alpha), y = "Group members selected") +
theme_dspa()
31.4 Problem 4: Construct an irrepresentable failure
Build designs that violate the condition by increasing amounts and measure recovery.
Solution
set.seed(209)
n4 <- 300; p4 <- 15; S4 <- 1:3
make_design <- function(rho) {
X <- matrix(rnorm(n4 * p4), n4, p4)
X[, 4] <- rho * rowSums(X[, S4]) / 3 + sqrt(1 - rho^2) * rnorm(n4)
X
}
p4tab <- do.call(rbind, lapply(seq(0, 0.98, length.out = 7), function(rho) {
X <- make_design(rho)
ic <- irrep(X, S4, rep(1, 3))
rec <- mean(replicate(15, {
b <- numeric(p4); b[S4] <- 2
y <- as.vector(X %*% b) + rnorm(n4)
sel <- which(as.vector(coef(cv.glmnet(X, y, alpha = 1, nfolds = 5),
s = "lambda.min"))[-1] != 0)
identical(sort(sel), S4)
}))
data.frame(rho = rho, irrepresentable = ic, holds = ic < 1, recovery = rec)
}))
p4tab |> mutate(across(where(is.numeric), \(z) round(z, 4)))ggplot(p4tab, aes(irrepresentable, recovery)) +
geom_vline(xintercept = 1, linetype = "dashed", color = "firebrick") +
geom_line(linewidth = 1, color = "steelblue") + geom_point(size = 2.6) +
annotate("text", x = 1.03, y = 0.9, hjust = 0, size = 3.2, color = "firebrick",
label = "condition fails") +
labs(title = "Exact support recovery against the irrepresentable quantity",
subtitle = "Recovery collapses as the quantity crosses 1, at fixed n and fixed signal strength",
x = "Irrepresentable quantity", y = "Exact recovery rate") +
theme_dspa()
31.5 Problem 5: Stability selection’s bound
Check the bound empirically across thresholds.
Solution
set.seed(211)
n5 <- 150; p5 <- 200; true5 <- 1:5
X5 <- matrix(rnorm(n5 * p5), n5, p5, dimnames = list(NULL, paste0("V", 1:p5)))
b5 <- numeric(p5); b5[true5] <- 2.5
q5 <- 15
p5tab <- do.call(rbind, lapply(c(0.55, 0.6, 0.7, 0.8, 0.9), function(pt) {
V <- replicate(8, {
y <- as.vector(X5 %*% b5) + rnorm(n5)
pi_h <- stability_selection(X5, y, B = 50, q_sel = q5,
seed = sample(1e6, 1))
sel <- as.integer(sub("V", "", names(which(pi_h >= pt))))
sum(!sel %in% true5)
})
data.frame(pi_thr = pt, observed_EV = mean(V),
bound = q5^2 / (p5 * (2 * pt - 1)))
}))
p5tab |> mutate(bound_respected = observed_EV <= bound,
across(where(is.numeric), \(z) round(z, 3)))p5tab |> pivot_longer(c(observed_EV, bound), names_to = "series", values_to = "v") |>
ggplot(aes(pi_thr, v, color = series)) +
geom_line(linewidth = 1) + geom_point(size = 2.4) +
scale_color_manual(values = c(observed_EV = "#3B7DD8", bound = "#D8433B")) +
labs(title = "Expected false selections against the Meinshausen-Buhlmann bound",
subtitle = "The bound is conservative; observed counts sit well below it",
x = expression(pi[thr]), y = "False selections", color = NULL) +
theme_dspa()
31.6 Problem 6: Antisymmetry is required
Show that a non-antisymmetric \(W_j\) breaks FDR control.
Solution
set.seed(213)
# VALID: swapping X_j and Xk_j flips the sign
stat_valid <- function(X, X_k, y)
abs(as.vector(crossprod(X, y))) - abs(as.vector(crossprod(X_k, y)))
# INVALID: always non-negative, so swapping cannot flip the sign
stat_broken <- function(X, X_k, y)
abs(as.vector(crossprod(X, y))) + abs(as.vector(crossprod(X_k, y)))
run_stat <- function(stat_fn, reps = 15, n = 300, p = 100, k = 10) {
vapply(seq_len(reps), function(r) {
X <- matrix(rnorm(n * p), n, p)
nz <- sample(p, k); b <- 4 * (seq_len(p) %in% nz) / sqrt(n)
y <- as.vector(X %*% b) + rnorm(n)
s <- tryCatch(knockoff.filter(X, y, fdr = 0.1,
knockoffs = create.second_order,
statistic = stat_fn)$selected,
error = function(e) integer(0))
sum(!s %in% nz) / max(1, length(s))
}, numeric(1))
}
c(target_fdr = 0.10,
mean_FDP_antisymmetric = round(mean(run_stat(stat_valid)), 4),
mean_FDP_not_antisymmetric = round(mean(run_stat(stat_broken)), 4))#> target_fdr mean_FDP_antisymmetric
#> 0.1 0.0
#> mean_FDP_not_antisymmetric
#> 0.931.7 Problem 7: Refit or shrink?
Compare LASSO predictions against OLS refitted on the selected support.
Solution
set.seed(215)
snr_grid <- c(0.5, 1, 2, 4, 8)
p7 <- do.call(rbind, lapply(snr_grid, function(snr) {
res <- replicate(20, {
n <- 120; p <- 60
X <- matrix(rnorm(n * p), n, p)
b <- c(rep(1, 5), rep(0, p - 5)) * snr / sqrt(5)
y <- as.vector(X %*% b) + rnorm(n)
Xt <- matrix(rnorm(2000 * p), 2000, p)
yt <- as.vector(Xt %*% b) + rnorm(2000)
cv <- cv.glmnet(X, y, alpha = 1, nfolds = 10)
sel <- which(as.vector(coef(cv, s = "lambda.min"))[-1] != 0)
pl <- as.vector(predict(cv, newx = Xt, s = "lambda.min"))
pr <- if (length(sel) && length(sel) < n - 2)
as.vector(cbind(1, Xt[, sel, drop = FALSE]) %*%
coef(lm(y ~ X[, sel, drop = FALSE]))) else rep(mean(y), 2000)
c(sqrt(mean((pl - yt)^2)), sqrt(mean((pr - yt)^2)))
})
data.frame(snr = snr, lasso = mean(res[1, ]), refit_ols = mean(res[2, ]))
}))
p7 |> mutate(refit_better = refit_ols < lasso,
across(where(is.numeric), \(z) round(z, 4)))p7 |> pivot_longer(c(lasso, refit_ols), names_to = "estimator", values_to = "rmse") |>
ggplot(aes(snr, rmse, color = estimator)) +
geom_line(linewidth = 1) + geom_point(size = 2.4) +
scale_x_log10(breaks = snr_grid) +
scale_color_manual(values = c(lasso = "#3B7DD8", refit_ols = "#D8433B")) +
labs(title = "Shrunk LASSO coefficients against OLS refitted on the selected support",
subtitle = "Refitting removes the shrinkage bias and adds variance; which wins depends on the signal-to-noise ratio",
x = "Signal-to-noise ratio (log scale)", y = "Test RMSE", color = NULL) +
theme_dspa()
31.8 Problem 8: Screening then selecting
Combine SIS with the LASSO and measure what each stage contributes.
Solution
set.seed(217)
n8 <- 150; p8 <- 10000; k8 <- 10
X8 <- matrix(rnorm(n8 * p8), n8, p8)
nz8 <- sample(p8, k8); b8 <- numeric(p8); b8[nz8] <- 3
y8 <- as.vector(X8 %*% b8) + rnorm(n8)
t_direct <- system.time({
cv_d <- cv.glmnet(X8, y8, alpha = 1, nfolds = 5)
sel_d <- which(as.vector(coef(cv_d, s = "lambda.1se"))[-1] != 0)
})[["elapsed"]]
t_sis <- system.time({
d8 <- floor(n8 / log(n8))
scr <- order(abs(crossprod(X8, y8)), decreasing = TRUE)[1:d8]
cv_s <- cv.glmnet(X8[, scr], y8, alpha = 1, nfolds = 5)
sel_s <- scr[which(as.vector(coef(cv_s, s = "lambda.1se"))[-1] != 0)]
})[["elapsed"]]
data.frame(
approach = c("LASSO on all p", "SIS then LASSO"),
features_entering_lasso = c(p8, floor(n8 / log(n8))),
n_selected = c(length(sel_d), length(sel_s)),
true_recovered = c(sum(sel_d %in% nz8), sum(sel_s %in% nz8)),
false_selected = c(sum(!sel_d %in% nz8), sum(!sel_s %in% nz8)),
seconds = round(c(t_direct, t_sis), 3))#> sure_screening_held
#> FALSE32 Checkpoint
- A colleague screened 20,000 genes down to 50 by \(t\)-test, then reported 10-fold CV accuracy of 0.92 for a classifier on those 50. What do you ask?
- Your LASSO selects
gene_Abut notgene_B, which correlates with it at 0.98. What can you conclude aboutgene_B? - You have \(n = 200\), \(p = 5000\), and believe about 10 predictors matter. Which method gives a guarantee, and what does it guarantee?
- Stability selection with \(\pi_{\mathrm{thr}} = 0.5\) returns 40 features. What is wrong?
- Boruta confirms 45 features; LASSO selects 12. Which is right?
- You want confidence intervals for the coefficients LASSO selected. What are your options?
Answers
- Whether the screening was repeated inside every fold. If the 50 genes were chosen using all 20,000 columns and all the labels, then every fold’s training set was selected with knowledge of that fold’s held-out labels, and the CV estimate is contaminated. On data with no signal at all this protocol routinely returns accuracy well above chance, and the inflation grows as \(p\) grows and \(n\) shrinks. The fix is nested: re-run the entire screening procedure within each training fold. Ask also for the sample size, the class balance, and the no-information rate.
- Very little. With two predictors correlated at 0.98, the LASSO’s \(\ell_1\) penalty makes selecting one nearly as good as selecting either, so it picks one essentially arbitrarily, and the choice flips under resampling or a small perturbation of the data. The irrepresentable condition is likely violated here, which means more data will not resolve it. Use the elastic net for the grouping effect, or stability selection to see how often each is chosen; if both are near 0.5, that is the finding.
- Knockoffs, which give exact finite-sample
FDR control: the expected proportion of false discoveries among
the selected set is at most \(q\), at
any \(n\), with no asymptotics,
provided the knockoffs are constructed correctly. With \(p \gg n\) use model-X knockoffs with
create.second_order(), andmethod = "asdp"since \(p^3\) at \(p=5000\) is prohibitive. Stability selection is the alternative: it bounds \(\mathbb{E}[V]\), the expected number of false selections, and wraps around any selector. - The bound does not apply. The Meinshausen–Bühlmann guarantee is \(\mathbb{E}[V]\le\frac{1}{2\pi_{\mathrm{thr}}-1}\cdot\frac{q^2}{p}\), whose denominator \(2\pi_{\mathrm{thr}}-1\) is zero at \(\pi_{\mathrm{thr}}=0.5\) and negative below it. The threshold must lie strictly in \((0.5, 1)\). Forty features is also a lot for a threshold that low, raise \(\pi_{\mathrm{thr}}\) and report the bound alongside the selection.
- Both, for different questions. Boruta is all-relevant: it aims to find every feature carrying information about the outcome, including redundant ones, so on correlated data it keeps whole groups. LASSO is minimal-optimal: it seeks the smallest set sufficient for prediction and drops duplicates. For biological interpretation, which pathways are involved, the all-relevant answer is what you want. For a parsimonious predictive model, the minimal-optimal one. The comparison that is meaningful is held-out predictive performance at each set size.
- Three, and the naive one is not among them.
Refitting OLS on the selected variables and reporting the usual
intervals is invalid, the model was chosen because it fit these data.
Instead: selective inference, which conditions on the
polyhedral selection event and gives truncated-Gaussian intervals
(
selectiveInference::fixedLassoInf); data splitting, selecting on one half and inferring on the other, which is always valid and costs efficiency; or abandon coefficient-level inference and use knockoffs or stability selection to control an error rate over the selection itself.
33 Summary
The problem
- Selection is discrete and dimension reduction continuous; the discrete choice adds estimator variance that everything downstream inherits.
- Selection bias is the cardinal error: screening on all the data and then cross-validating the model produces excellent apparent accuracy on pure noise. The entire selection procedure must run inside every training fold.
- Univariate filters are \(O(np)\) and blind to interactions. SIS retains all true predictors with probability tending to 1, but only those that are marginally associated.
Regularization
- Under orthonormal design, ridge scales every coefficient by \((1+2n\lambda)^{-1}\) and LASSO soft-thresholds at \(n\lambda\), for the \(\frac{1}{2n}\) fidelity convention. The threshold constant follows the convention, and the two must match.
- Best subset (\(\ell_0\), exact, NP-hard) and OMP (greedy, polynomial) are different algorithms.
- The elastic net needs both \(\alpha\) (mix) and \(\lambda\) (weight). Its \(\ell_2\) component produces the grouping effect, selecting correlated predictors together.
- A single thresholding pass is not the LASSO unless the design is orthonormal. Coordinate descent iterates to convergence.
When selection succeeds
- LASSO selection consistency requires the irrepresentable condition on the design. When it fails, more data does not help.
- Support recovery has a phase transition at roughly \(n\propto s\log p\).
- \(\mathrm{df}(\hat\beta^{L})=\mathbb E[\#\{j:\hat\beta_j\ne0\}]\) exactly, which makes \(C_p\)/AIC computable along the path.
- Naive p-values after selection are invalid. Use selective inference, data splitting, or control an error rate over the selection instead.
Error control
- BH is a step-up rule: find \(\hat k=\max\{i:p_{(i)}\le qi/m\}\) and reject everything up to it.
- Stability selection bounds \(\mathbb E[V]\le\frac{q^2}{p(2\pi_{\mathrm{thr}}-1)}\) and requires \(\pi_{\mathrm{thr}}>0.5\).
- Knockoffs give exact finite-sample FDR control. They need exchangeable decoys, an antisymmetric \(W_j\), and a correctly specified model for \(X\), asserting \(\Sigma=I\) on correlated data voids the guarantee.
- Verify FDR control by replication; the guarantee is about an expectation.
Choosing a method
- All-relevant (Boruta) and minimal-optimal (LASSO) answer different questions. Compare them on held-out prediction, not on overlap.
- Only knockoffs give a finite-sample guarantee; everything else is asymptotic, resampling-based, or uncontrolled.
Where these threads continue
| Thread | Continues in |
|---|---|
| Selection for repeated-measures and longitudinal designs | Longitudinal analysis |
| The optimization behind coordinate descent and ADMM | Function optimization |
| Sparsity and pruning in learned representations | Deep learning |
34 Chapter roadmap
- Chapter 1: Foundations. R toolchain, reproducibility conventions, simulation.
- Chapter 2: Data quality and exploratory visual analytics. Missingness, multiplicity, FDR.
- Chapter 3: Linear algebra and regression. SVD, conditioning, the hat matrix.
- Chapter 4: Dimensionality reduction. PCA, spurious correlation when \(p\gg n\).
- Chapter 5: Supervised classification. Evaluation metrics, the leakage taxonomy.
- Chapter 6: Black-box methods. Kernels, ensembles, regularized networks.
- Chapter 7: Text mining and association rules. Multiplicity control over an exponential search.
- Chapter 8: Unsupervised clustering. Internal and external validation, stability.
- Chapter 9: Model assessment and validation. Optimism, calibration, nested resampling.
- Chapter 10: Formats, streams, performance. Columnar storage, sketches, Amdahl’s law.
- Longitudinal and time-series analysis. Mixed models, ARIMA, forecast evaluation.
- Function optimization. Coordinate descent, proximal methods, duality.
- Deep learning. Sparsity, pruning, representation learning.
35 Session information
#> R version 4.3.3 (2024-02-29 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 11 x64 (build 26200)
#>
#> Matrix products: default
#>
#>
#> locale:
#> [1] LC_COLLATE=English_United States.utf8
#> [2] LC_CTYPE=English_United States.utf8
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United States.utf8
#>
#> time zone: America/New_York
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] mRMRe_2.1.2.2 igraph_2.0.3 survival_3.7-0 Boruta_8.0.0
#> [5] ranger_0.16.0 caret_6.0-94 lattice_0.22-6 knockoff_0.3.6
#> [9] rsample_1.2.1 glmnet_4.1-8 Matrix_1.6-5 plotly_4.12.1
#> [13] patchwork_1.3.0 tidyr_1.3.1 dplyr_1.1.4 ggplot2_4.0.1
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.1 viridisLite_0.4.2 timeDate_4032.109
#> [4] farver_2.1.2 S7_0.2.1 fastmap_1.2.0
#> [7] pROC_1.18.5 digest_0.6.37 rpart_4.1.23
#> [10] timechange_0.3.0 lifecycle_1.0.5 magrittr_2.0.3
#> [13] compiler_4.3.3 rlang_1.1.5 sass_0.4.9
#> [16] tools_4.3.3 yaml_2.3.10 data.table_1.16.4
#> [19] knitr_1.51 labeling_0.4.3 htmlwidgets_1.6.4
#> [22] plyr_1.8.9 RColorBrewer_1.1-3 selectiveInference_1.2.5
#> [25] withr_3.0.2 purrr_1.0.2 stats4_4.3.3
#> [28] nnet_7.3-19 grid_4.3.3 Rdsdp_1.0.5.2.1
#> [31] future_1.33.2 globals_0.16.3 scales_1.4.0
#> [34] iterators_1.0.14 MASS_7.3-60.0.1 isoband_0.2.7
#> [37] cli_3.6.3 rmarkdown_2.31 intervals_0.15.5
#> [40] generics_0.1.3 otel_0.2.0 rstudioapi_0.18.0
#> [43] future.apply_1.11.2 RSpectra_0.16-1 reshape2_1.4.4
#> [46] httr_1.4.7 cachem_1.1.0 stringr_1.5.1
#> [49] splines_4.3.3 parallel_4.3.3 vctrs_0.6.5
#> [52] hardhat_1.4.3 jsonlite_1.8.9 listenv_0.9.1
#> [55] crosstalk_1.2.1 foreach_1.5.2 gower_1.0.1
#> [58] jquerylib_0.1.4 recipes_1.4.0 glue_1.8.0
#> [61] parallelly_1.37.1 codetools_0.2-20 stringi_1.8.4
#> [64] lubridate_1.9.3 shape_1.4.6.1 gtable_0.3.6
#> [67] tibble_3.2.1 furrr_0.3.1 pillar_1.10.1
#> [70] adaptMCMC_1.5 htmltools_0.5.8.1 randomForest_4.7-1.1
#> [73] ipred_0.9-14 lava_1.8.0 R6_2.6.1
#> [76] doParallel_1.0.17 evaluate_1.0.3 bslib_0.9.0
#> [79] class_7.3-22 Rcpp_1.0.14 coda_0.19-4.1
#> [82] nlme_3.1-165 prodlim_2024.06.25 xfun_0.52
#> [85] ModelMetrics_1.2.2.2 pkgconfig_2.0.3