Inference after a glmnet fit without randomization

The simplest use of lassoinf takes a standard glmnet fit on all the data, at a given \(\lambda\), and gives intervals and p-values for the selected coefficients that account for the selection. Nothing is held out and no noise is added. This is the polyhedral approach of Lee et al. (2016), here with the linearization of a GLM used throughout the package. glmnet_inference does it in one call; we also show what it computes.

Setup

glmnet divides the loss by the number of observations \(n\). glmnet_problem_from_fit returns the LASSO solution \(\hat{\beta}\), the gradient \(G\) of the smooth part at \(\hat{\beta}\), and the Hessian \(\hat{Q}\), all in the original coordinates \((\text{intercept}, \beta)\). Linearizing the KKT conditions, the selection is a function of the score

\[ Z = -G + \hat{Q}\hat{\beta}, \]

a one-step estimator, \(\hat{Q}\) times an approximation to the unpenalized MLE. With \(I\) the per-observation Fisher information, \(\text{Var}(Z) \approx I / n\). We take \(\Sigma = \hat{Q} / n\). By default glmnet_problem_from_fit evaluates \(\hat{Q}\) at the relaxed fit: one Newton step for the unpenalized loss on the intercept and the selected variables \(E\) from the LASSO solution, \(\tilde{\beta}_E = \hat{\beta}_E - H_{E,E}^{-1} \nabla \ell(\hat{\beta})_E\), with \(H\) the Hessian at \(\hat{\beta}\). At the LASSO solution \(\nabla \ell(\hat{\beta})_E = -D_E \, \text{sign}(\hat{\beta}_E)\), so the step undoes the shrinkage of \(\hat{\beta}\) toward zero, which makes the information at \(\hat{\beta}\) overstate \(I\). Pass information = "lasso" to use \(\hat{Q} = H\).

There is no randomization: the score used for selection is \(Z\) itself, which is scalar_noise = 0.

Data and fit

library(lassoinf)
set.seed(0)

n <- 600; p <- 30
beta <- c(0.6, -0.6, 0.4, rep(0, p - 3))
make_data <- function(n) {
  X <- matrix(rnorm(n * p), n, p)
  y <- rbinom(n, 1, plogis(drop(X %*% beta)))
  list(X = X, y = y)
}
d <- make_data(n)

lam <- 0.03
lambda_path <- lam * exp(seq(log(10), 0, length.out = 30))
fit <- glmnet::glmnet(d$X, d$y, family = "binomial", lambda = lambda_path, control = list(thresh = 1e-14))

We fit a path ending at \(\lambda\) with a small thresh: glmnet_problem_from_fit warns if the fit does not satisfy the KKT conditions closely enough.

Inference

li <- glmnet_inference(fit, d$X, d$y, s = lam, level = 0.9)
summ <- li$summary()
summ$variable <- ifelse(summ$index == 0, "(Intercept)", paste0("x", summ$index))
summ
#>    index    beta_hat lower_conf upper_conf      p_value    variable
#> 0      0 -0.04390486 -0.1926541  0.5088327 8.227524e-01 (Intercept)
#> 1      1  0.60206784  0.7002772  1.0398212 7.327472e-15          x1
#> 2      2 -0.37446602 -0.7584731 -0.4454181 2.797866e-08          x2
#> 3      3  0.23466365  0.2731540  0.5905153 1.765597e-04          x3
#> 7      7  0.01072090 -0.5636646  0.2200132 4.570290e-01          x7
#> 21    21 -0.16211621 -0.5225519 -0.1506638 1.047598e-02         x21

Index 0 is the intercept. The estimate for the selected model \(E\) is \(W Z_E\) with \(W = \hat{Q}_{E,E}^{-1}\), and its target is \(W (\hat{Q}\beta^*)_E\), with \(\beta^* = (0, \beta)\). When \(E\) contains all the true signals, this is \(\beta^*_E\).

target <- function(li) drop(li$W %*% (li$Q_hat %*% c(0, beta))[li$E + 1])
summ$truth <- target(li)
summ$cover <- summ$lower_conf <= summ$truth & summ$truth <= summ$upper_conf
summ[, c("variable", "lower_conf", "upper_conf", "p_value", "truth", "cover")]
#>       variable lower_conf upper_conf      p_value         truth cover
#> 0  (Intercept) -0.1926541  0.5088327 8.227524e-01 -3.870803e-18  TRUE
#> 1           x1  0.7002772  1.0398212 7.327472e-15  6.000000e-01 FALSE
#> 2           x2 -0.7584731 -0.4454181 2.797866e-08 -6.000000e-01  TRUE
#> 3           x3  0.2731540  0.5905153 1.765597e-04  4.000000e-01  TRUE
#> 7           x7 -0.5636646  0.2200132 4.570290e-01  3.330745e-18  TRUE
#> 21         x21 -0.5225519 -0.1506638 1.047598e-02  7.648649e-18 FALSE

For comparison, li$"_naive" holds Wald intervals that ignore selection:

li$`_naive`
#>    index    beta_hat  lower_conf  upper_conf      p_value
#> 0      0 -0.06214006 -0.21203798  0.08775786 4.953194e-01
#> 1      1  0.87004951  0.70027786  1.03982116 3.469010e-17
#> 2      2 -0.60219853 -0.75846975 -0.44592731 2.319811e-10
#> 3      3  0.43744706  0.28436545  0.59052866 2.597222e-06
#> 7      7  0.16635893  0.01732356  0.31539431 6.635057e-02
#> 21    21 -0.36205334 -0.52299283 -0.20111385 2.153428e-04

glmnet_inference is the following, with prob the LASSO problem glmnet solved and scalar_noise = 0 for no randomization:

prob <- glmnet_problem_from_fit(fit, d$X, d$y, s = lam)
Z <- -prob$G_hat + drop(prob$Q_hat %*% prob$beta_hat)
li_manual <- LassoInference$new(prob$beta_hat, prob$G_hat, prob$Q_hat, prob$D, prob$L, prob$U,
                                Z_full = Z, Sigma = prob$Q_hat / n, scalar_noise = 0,
                                level = 0.9)
all.equal(li_manual$summary(), li$summary())
#> [1] TRUE

For a gaussian response

For family = "gaussian" the loss is \(\frac{1}{2n}\|y - X\beta\|^2\), so \(\hat{Q} = X'X / n\) and \(\text{Var}(Z) = \sigma^2 \hat{Q} / n\). glmnet_inference estimates \(\sigma^2\) from the full least squares fit when \(n > p\); otherwise pass dispersion. The relaxed and LASSO information agree.

Coverage over repeated experiments

one_replication <- function() {
  d <- make_data(n)
  fit <- glmnet::glmnet(d$X, d$y, family = "binomial", lambda = lambda_path, control = list(thresh = 1e-14))
  li <- glmnet_inference(fit, d$X, d$y, s = lam, level = 0.9)
  summ <- li$summary()
  keep <- summ$index != 0   # report the selected variables, not the intercept
  truth <- target(li)[keep]
  tabs <- list(polyhedral = summ[keep, ], naive = li$`_naive`[keep, ])
  data.frame(method = rep(names(tabs), each = sum(keep)),
             cover = unlist(lapply(tabs, function(df) df$lower_conf <= truth & truth <= df$upper_conf)),
             length = unlist(lapply(tabs, function(df) df$upper_conf - df$lower_conf)),
             signal = rep(beta[summ$index[keep]] != 0, length(tabs)))
}

set.seed(1)
results <- do.call(rbind, replicate(100, one_replication(), simplify = FALSE))
with(results, tapply(cover, list(method, signal = signal), mean))
#>             signal
#>                  FALSE      TRUE
#>   naive      0.0433213 0.8266667
#>   polyhedral 0.8880866 0.8700000
with(results, tapply(length, list(method, signal = signal), median))
#>             signal
#>                  FALSE      TRUE
#>   naive      0.2995776 0.3128654
#>   polyhedral 0.7467743 0.3401971

The naive intervals almost never cover for the selected null variables: they were selected because their estimates are large. The polyhedral intervals cover close to the nominal 90% for both, a little below on the true signals, which reflects the large-sample approximations behind the Gaussian model for \(Z\). For strong signals, far from the selection boundary, they are about as long as the naive intervals. For variables selected near the boundary they are much longer. Randomization, for example carving as in vignette("glmnet_logistic", package = "lassoinf"), gives shorter intervals there.

Reference

Lee, J. D., Sun, D. L., Sun, Y. and Taylor, J. E. (2016). Exact post-selection inference, with application to the lasso. The Annals of Statistics, 44(3), 907–927.