After the LASSO selects an active set \(E\), it is natural to ask about the variables it left out. A common heuristic looks at the KKT “near misses”: the inactive variables whose gradient at the LASSO solution is largest in absolute value. Here we keep the top 10 inactive variables by \(|\nabla \ell(\hat{\beta})_j|\) and, for each, test its coefficient in the model \(E \cup \{j\}\).
Both the LASSO and the screening step look at the data, so naive
tests of the screened variables are invalid.
ScreenedSelection adds the screening step to the selection
event, and the polyhedral lemma then gives valid intervals and
p-values.
With \(Q = X'X\) and \(W = Q_{E,E}^{-1}\), the LASSO selection event depends on the data only through
On the selection event the inactive gradient at the LASSO solution is \(\nabla \ell(\hat{\beta})_{-E} = -U_{-E}(Z + \omega) + \text{const}\), an affine function of the same coordinates, where \(Z + \omega\) is the randomized score used for selection. Keeping the top \(K\) variables by \(|\nabla \ell(\hat{\beta})_j|\), with the signs of their gradients, is then also a polyhedral event:
conditioning = "first_dropped"
(default): also condition on the largest dropped variable \(l^*\) and its sign \(s^*\). This gives the rows \(s_k G_k \geq s^* G_{l^*}\) for kept \(k\) and \(|G_l|
\leq s^* G_{l^*}\) for the other dropped \(l\), which is \(O(p)\) rows.conditioning = "exact": the pairwise
rows \(s_k G_k \geq |G_l|\) for every
kept \(k\) and dropped \(l\), which is \(2K(p - |E| - K)\) rows. It conditions on
less, at the price of more constraints.The coefficient of inactive variable \(j\) in the model \(E \cup \{j\}\) is \(U_j / S_{jj}\) with \(S_{jj} = Q_{jj} - Q_{j,E} W Q_{E,j}\): a contrast \(\eta'Z\) in the same coordinates. All of this is computed with matrix-vector products by \(Q\), so \(X'X\) is never formed.
We take \(n = 100\) observations and \(p = 300\) features, with 15 nonzero coefficients of moderate size, and treat \(\sigma^2 = 1\) as known. Selection uses a randomized response \(y_{noisy} = y + \epsilon\) with \(\epsilon \sim N(0, \frac{1-\pi}{\pi}\sigma^2 I)\) and \(\pi = 0.8\).
library(lassoinf)
n <- 100; p <- 300; k <- 15
pi_selection <- 0.8
scalar_noise <- (1 - pi_selection) / pi_selection # Var(omega) = scalar_noise * Sigma
make_data <- function() {
X <- matrix(rnorm(n * p), n, p)
beta <- numeric(p)
beta[1:k] <- 3 / sqrt(n) * sample(c(-1, 1), k, replace = TRUE)
mu <- drop(X %*% beta)
y <- mu + rnorm(n)
y_noisy <- y + rnorm(n, sd = sqrt(scalar_noise))
list(X = X, beta = beta, mu = mu, y = y, y_noisy = y_noisy)
}
set.seed(0)
d <- make_data()We fit \(\frac{1}{2}\|y_{noisy} - Xb\|^2 +
\lambda \|b\|_1\) with glmnet, which scales the loss
by \(1/n\), so we pass it \(\lambda / n\). \(Q = X'X\) and \(\Sigma = \sigma^2 X'X\) are matrix-free
operators.
lam <- 0.8 * sqrt(2 * n * log(p))
lasso_inference <- function(d) {
fit <- glmnet::glmnet(d$X, d$y_noisy, lambda = lam / n, intercept = FALSE,
standardize = FALSE, control = list(thresh = 1e-14))
beta_hat <- as.numeric(coef(fit))[-1]
beta_hat[abs(beta_hat) < 1e-6] <- 0
G_hat <- drop(crossprod(d$X, d$X %*% beta_hat - d$y_noisy))
Q <- xtvx_operator(d$X)
LassoInference$new(beta_hat, G_hat, Q, rep(lam, p), NULL, NULL,
drop(crossprod(d$X, d$y)), Q, NULL, scalar_noise, level = 0.9)
}
li <- lasso_inference(d)
E <- li$E + 1 # the package reports 0-based indices
cat(length(E), "active variables,", sum(d$beta[E] != 0), "of them true signals;",
sum(d$beta[-E] != 0), "true signals were left out\n")
#> 19 active variables, 9 of them true signals; 6 true signals were left outscreen <- ScreenedSelection$new(li, top_k = 10)
screen$screened + 1
#> [1] 189 162 44 294 61 2 106 283 244 109
screen$screened_signs
#> [1] -1 -1 -1 -1 1 -1 1 1 1 -1
c(first_dropped = screen$first_dropped + 1, sign = screen$first_dropped_sign)
#> first_dropped sign
#> 17 1
c(lasso_rows = li$A$rows(), screening_rows = screen$A_screen$rows())
#> lasso_rows screening_rows
#> 581 550inactive_summary gives the coefficient of each screened
variable in the model \(E \cup \{j\}\),
conditional on both selection steps. Its target is the population
coefficient \(\eta_j' X'\mu\),
the coefficient of \(x_j\) in the
regression of \(\mu = E[y]\) on \(X_{E \cup \{j\}}\). Its index
column is 0-based.
coords <- SelectionCoordinates$new(li)
Z_mu <- drop(crossprod(d$X, d$mu))
screened_df <- inactive_summary(screen)
screened_df$truth <- sapply(screened_df$index, function(j) estimand_value(coords$inactive_coef(j), Z_mu))
screened_df$signal <- d$beta[screened_df$index + 1] != 0
screened_df
#> index estimate lower_conf upper_conf p_value truth signal
#> 1 188 0.05345127 -0.6260857 0.2062085 0.407933313 0.01669426 FALSE
#> 2 161 0.07852045 -0.6665216 0.1821812 0.347894500 -0.02136290 FALSE
#> 3 43 0.09475950 -0.1453832 0.6777984 0.287533164 0.03568281 FALSE
#> 4 293 0.08365067 -0.6433786 0.1338483 0.280942388 0.06358084 FALSE
#> 5 60 -0.32050108 -1.0656255 -0.2459266 0.008475311 0.02334345 FALSE
#> 6 1 0.13222270 -0.6541625 0.2259877 0.423620180 0.10483530 TRUE
#> 7 105 -0.20172439 -0.5126651 0.2701743 0.610137441 -0.00571950 FALSE
#> 8 282 -0.07898855 -0.3933067 0.3385964 0.902141289 -0.13044266 FALSE
#> 9 243 -0.19587896 -0.5974074 0.1855058 0.386711164 -0.01032484 FALSE
#> 10 108 0.07080740 -0.4755348 0.3123230 0.733555690 0.04547556 FALSEUsing the LASSO constraints alone treats the 10 variables as if they had been chosen in advance, and naive intervals ignore selection altogether. Both are shown for the same variables:
naive_p_value <- function(li, est) {
sd <- sqrt(li$si$compute_contrast(est$eta)$naive_variance)
2 * pnorm(-abs(estimand_value(est, li$Z_full)) / sd)
}
data.frame(variable = screen$screened + 1,
screened = inactive_summary(screen)$p_value,
lasso_only = inactive_summary(li, variables = screen$screened)$p_value,
naive = sapply(screen$screened, function(j) naive_p_value(li, coords$inactive_coef(j))),
signal = d$beta[screen$screened + 1] != 0)
#> variable screened lasso_only naive signal
#> 1 189 0.407933313 0.4841260049 0.648586668 FALSE
#> 2 162 0.347894500 0.4561303544 0.497085348 FALSE
#> 3 44 0.287533164 0.3097101103 0.406756644 FALSE
#> 4 294 0.280942388 0.3471326764 0.432193496 FALSE
#> 5 61 0.008475311 0.0004677782 0.004563512 FALSE
#> 6 2 0.423620180 0.2642190567 0.270856858 TRUE
#> 7 106 0.610137441 0.0663404152 0.062894733 FALSE
#> 8 283 0.902141289 0.6421688279 0.428212923 FALSE
#> 9 244 0.386711164 0.2244077368 0.074154618 FALSE
#> 10 109 0.733555690 0.8226034959 0.515091813 FALSEIgnoring the screening step tends to give smaller p-values, since variables are kept because their gradients are large. The repeated experiments below show that those smaller p-values are not valid.
The same estimands can be centered at the LASSO solution instead of
the refit: with basis = "lasso" the coefficient is the
one-step estimate \(-\nabla
\ell(\hat{\beta})_j / S_{jj}\). It uses the same contrast, so the
intervals are shifted copies of the refit intervals, and the estimand is
shifted by the same constant. The p-values change because \(\theta = 0\) is now a null hypothesis about
the shifted estimand.
head(inactive_summary(screen, basis = "lasso"))
#> index estimate lower_conf upper_conf p_value
#> 1 188 0.2668280 -0.4127089 0.4195853 0.986359660
#> 2 161 0.2692984 -0.4757436 0.3729591 0.842114625
#> 3 43 0.3626536 0.1225110 0.9456925 0.032789959
#> 4 293 0.1813887 -0.5456406 0.2315863 0.506402483
#> 5 60 -0.3944575 -1.1395819 -0.3198830 0.003394698
#> 6 1 0.2465080 -0.5398772 0.3402731 0.709159031exact <- ScreenedSelection$new(li, top_k = 10, conditioning = "exact")
c(first_dropped = screen$A_screen$rows(), exact = exact$A_screen$rows())
#> first_dropped exact
#> 550 5420
cols <- c("lower_conf", "upper_conf", "p_value")
cbind(index = screen$screened,
first_dropped = inactive_summary(screen)[, cols],
exact = inactive_summary(exact)[, cols])
#> index first_dropped.lower_conf first_dropped.upper_conf
#> 1 188 -0.6260857 0.2062085
#> 2 161 -0.6665216 0.1821812
#> 3 43 -0.1453832 0.6777984
#> 4 293 -0.6433786 0.1338483
#> 5 60 -1.0656255 -0.2459266
#> 6 1 -0.6541625 0.2259877
#> 7 105 -0.5126651 0.2701743
#> 8 282 -0.3933067 0.3385964
#> 9 243 -0.5974074 0.1855058
#> 10 108 -0.4755348 0.3123230
#> first_dropped.p_value exact.lower_conf exact.upper_conf exact.p_value
#> 1 0.407933313 -0.6260857 0.2062085 0.407933313
#> 2 0.347894500 -0.6665216 0.1821812 0.347894500
#> 3 0.287533164 -0.1453832 0.6777984 0.287533164
#> 4 0.280942388 -0.6433786 0.1338483 0.280942388
#> 5 0.008475311 -1.0656255 -0.2459266 0.008475311
#> 6 0.423620180 -0.6541625 0.2259877 0.423620180
#> 7 0.610137441 -0.5061937 0.2744241 0.624927146
#> 8 0.902141289 -0.3933067 0.3385964 0.902141289
#> 9 0.386711164 -0.5974074 0.1855058 0.386711164
#> 10 0.733555690 -0.4755348 0.3123230 0.733555690Conditioning on the first-dropped variable restricts the selection event further, so its truncation intervals are subsets of the exact ones, and often the same.
We repeat the experiment and record, for each screened variable, whether the 90% interval covers its target. With 50 replications that is 500 intervals per method.
one_replication <- function() {
d <- make_data()
li <- lasso_inference(d)
screen <- ScreenedSelection$new(li, top_k = 10)
coords <- SelectionCoordinates$new(li)
Z_mu <- drop(crossprod(d$X, d$mu))
q <- qnorm(0.95)
rows <- list()
for (j in screen$screened) {
est <- coords$inactive_coef(j)
target <- estimand_value(est, Z_mu)
sd <- sqrt(li$si$compute_contrast(est$eta)$naive_variance)
theta <- estimand_value(est, li$Z_full)
signal <- d$beta[j + 1] != 0
for (method in c("screened", "lasso only")) {
res <- contrast_inference(if (method == "screened") screen else li, est)
rows[[length(rows) + 1]] <- data.frame(method = method,
cover = res$lower_conf <= target && target <= res$upper_conf,
length = res$upper_conf - res$lower_conf, reject = res$p_value < 0.1, signal = signal)
}
rows[[length(rows) + 1]] <- data.frame(method = "naive", cover = abs(theta - target) <= q * sd,
length = 2 * q * sd, reject = abs(theta) > q * sd, signal = signal)
}
do.call(rbind, rows)
}
set.seed(1)
results <- do.call(rbind, replicate(50, one_replication(), simplify = FALSE))
data.frame(coverage = tapply(results$cover, results$method, mean),
median_length = tapply(results$length, results$method, median))
#> coverage median_length
#> lasso only 0.824 0.6987183
#> naive 0.858 0.3595961
#> screened 0.912 0.7932155Only the intervals that condition on the screening step come close to the nominal 90% coverage. Rejection rates at level 0.1, split by whether the screened variable is a true signal:
with(results, tapply(reject, list(method, signal = signal), mean))
#> signal
#> FALSE TRUE
#> lasso only 0.2622951 0.5753425
#> naive 0.2786885 0.5479452
#> screened 0.1077283 0.2602740Among screened null variables, the screened tests reject at about the nominal 10%, while the naive and LASSO-only tests reject far more often. The screened tests still find a share of the true signals the LASSO missed.