SOCR ≫ DSPA ≫ DSPA3 Topics ≫ 

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 ggplot2 and rendered statically. Immediately after each one, the equivalent plot_ly() code appears in a chunk marked eval=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:

  1. Distinguish feature selection from dimension reduction, and filter, wrapper, and embedded strategies from one another.
  2. Recognize and avoid selection bias, and demonstrate that selecting on all the data can produce excellent cross-validated accuracy on pure noise.
  3. Derive the ridge and LASSO closed forms under orthonormal design and state the correct soft-threshold level.
  4. Explain the elastic net’s grouping effect and when it is preferable to the LASSO.
  5. Implement cyclic coordinate descent and validate it against glmnet.
  6. State the irrepresentable condition and explain what it does and does not guarantee about support recovery.
  7. Use the exact degrees-of-freedom result for the LASSO to connect the path to information criteria.
  8. Apply Benjamini–Hochberg correctly as a step-up procedure, and explain why the rejection set is contiguous by construction.
  9. Run stability selection and state its bound on the expected number of false selections.
  10. 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.3650

There 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.857797
c(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.7
data.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                        1

All 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")))