Inference after a glmnet fit: carving for logistic regression

This vignette shows carving in logistic regression. We select variables with glmnet on a random half of the data, then do inference with all of it, accounting for selection. glmnet_inference does this from the fitted glmnet object. It uses glmnet_problem_from_fit, which returns the LASSO problem glmnet solved in the original coordinates: the solution, the gradient and Hessian of the smooth part, the penalty weights and the bounds. For inference without held-out data see vignette("no_randomization", package = "lassoinf").

Carving as randomization

glmnet divides the loss by the number of observations, and we use this per-observation scale throughout. On the subsample of size $m = \pi n$, glmnet’s problem gives the LASSO solution $\hat{\beta}$, the gradient $G$ at $\hat{\beta}$, and the information $\hat{Q}$ at the one-step relaxed fit on the selected model $E$ (see ?glmnet_problem). LassoInference linearizes the selection with $\hat{Q}$ and selects with

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

The full-data statistic uses the same matrix and the full-data gradient $G_{full}$, at one Newton step for the full-data loss on $E$ from $\hat{\beta}$, $\tilde{\beta} = \hat{\beta} - H_{E,E}^{-1} G_{full,E}$:

\[ Z = -G_{full}(\tilde{\beta}) + \hat{Q}\tilde{\beta}. \]

Both are $\hat{Q}$ times approximations to the unpenalized MLE, on the subsample and on all the data. With $I$ the per-observation Fisher information, the subsample MLE minus the full-data MLE is asymptotically independent of the full-data MLE, with variance $I^{-1}(1/m - 1/n)$, and $\hat{Q} \approx I$. So

\[ \Sigma = \text{Var}(Z) \approx \hat{Q} / n, \qquad \text{Var}(\omega) \approx \frac{1 - \pi}{\pi}\, \Sigma, \]

and selecting on a subsample is randomized selection with scalar_noise $= (1 - \pi) / \pi$, which is 1 for $\pi = 1/2$. glmnet_inference with selection_rows does all of this.

Data

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)

Select on half of the data with glmnet

We fit a path that ends at a fixed $\lambda$, with the default standardization and an unpenalized intercept. The coordinates are $(\text{intercept}, \beta)$, so index 0 is the intercept. glmnet_inference warns if the fit does not satisfy the KKT conditions closely enough, so a small thresh matters.

pi_selection <- 0.5
lam <- 0.06
lambda_path <- lam * exp(seq(log(10), 0, length.out = 30))

select <- function(d) {
  rows <- sample(nrow(d$X), round(pi_selection * nrow(d$X)))
  fit <- glmnet::glmnet(d$X[rows, ], d$y[rows], family = "binomial",
                        lambda = lambda_path, control = list(thresh = 1e-14))
  list(fit = fit, rows = rows)
}
sel <- select(d)
which(coef(sel$fit, s = lam)[-1] != 0)
#> [1] 1 2 3

Inference

li <- glmnet_inference(sel$fit, d$X, d$y, s = lam, selection_rows = sel$rows, level = 0.9)
li$scalar_noise
#> [1] 1
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.09794322 -0.1604769  0.2091305 8.675623e-01 (Intercept)
#> 1     1  0.37009397  0.6443921  0.9557616 2.042810e-14          x1
#> 2     2 -0.25187169 -0.7054899 -0.3597215 3.304804e-06          x2
#> 3     3  0.09067492  0.1822282  0.5512053 1.654996e-03          x3

The estimate for the selected model is $W Z_E$, with $W = \hat{Q}_{E,E}^{-1}$, and its target is $W (\hat{Q}\beta^*)_E$, with $\beta^* = (0, \beta)$ the true coefficients including the intercept. When $E$ contains all the true signals, the target is $\beta^*_E$. Otherwise it is the projection of $\beta^*$ onto the coordinates in $E$, in the metric given by the information matrix $\hat{Q}$.

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.1604769  0.2091305 8.675623e-01 -1.796405e-20  TRUE
#> 1          x1  0.6443921  0.9557616 2.042810e-14  6.000000e-01 FALSE
#> 2          x2 -0.7054899 -0.3597215 3.304804e-06 -6.000000e-01  TRUE
#> 3          x3  0.1822282  0.5512053 1.654996e-03  4.000000e-01  TRUE

For comparison, li$"_naive" holds Wald intervals that use all the data and ignore selection. Data splitting refits the selected model on the half of the data not used for selection:

splitting <- function(d, rows, li, level = 0.9) {
  held_out <- setdiff(seq_len(nrow(d$X)), rows)
  X1 <- cbind(1, d$X)[held_out, li$E + 1, drop = FALSE]
  fit <- glm.fit(X1, d$y[held_out], family = binomial())
  est <- unname(coef(fit))
  mu <- fit$fitted.values
  sd <- sqrt(diag(solve(crossprod(X1, X1 * (mu * (1 - mu))))))
  q <- qnorm(1 - (1 - level) / 2)
  data.frame(index = li$E, estimate = est, lower_conf = est - q * sd, upper_conf = est + q * sd)
}
li$`_naive`
#>   index    beta_hat lower_conf  upper_conf      p_value
#> 0     0 -0.05541122 -0.2022343  0.09141183 5.347510e-01
#> 1     1  0.80118217  0.6465780  0.95578637 1.542972e-17
#> 2     2 -0.55326082 -0.7105093 -0.39601230 7.155738e-09
#> 3     3  0.42623038  0.2744639  0.57799682 3.846473e-06
splitting(d, sel$rows, li)
#>   index    estimate lower_conf upper_conf
#> 1     0  0.04189285 -0.1713128  0.2550985
#> 2     1  0.91872578  0.6617891  1.1756625
#> 3     2 -0.53284884 -0.7485171 -0.3171805
#> 4     3  0.44722400  0.2334200  0.6610280

The held-out refit targets the coefficients of the best-fitting logistic model on $X_E$. That is the same as the target above when $E$ contains all the true signals.

Coverage over repeated experiments

one_replication <- function() {
  d <- make_data(n)
  sel <- select(d)
  li <- glmnet_inference(sel$fit, d$X, d$y, s = lam, selection_rows = sel$rows, level = 0.9)
  summ <- li$summary()
  keep <- summ$index != 0   # report the selected variables, not the intercept
  truth <- target(li)[keep]
  covers <- function(df) (df$lower_conf <= truth & truth <= df$upper_conf)
  lengths <- function(df) df$upper_conf - df$lower_conf
  tabs <- list(carving = summ[keep, ], naive = li$`_naive`[keep, ],
               splitting = splitting(d, sel$rows, li)[keep, ])
  data.frame(method = rep(names(tabs), each = sum(keep)),
             cover = unlist(lapply(tabs, covers)), length = unlist(lapply(tabs, lengths)),
             signal = rep(beta[summ$index[keep]] != 0, length(tabs)))
}

set.seed(1)
results <- do.call(rbind, replicate(100, one_replication(), simplify = FALSE))
data.frame(coverage = tapply(results$cover, results$method, mean),
           median_length = tapply(results$length, results$method, median))
#>            coverage median_length
#> carving   0.8969359     0.3641605
#> naive     0.8328691     0.3093225
#> splitting 0.9247911     0.4400042

Carving and data splitting both have close to the nominal 90% coverage. Carving also uses the half of the data that selected the model, so its intervals are shorter. The naive intervals ignore selection and undercover, almost entirely on the selected null variables:

with(results, tapply(cover, list(method, signal = signal), mean))
#>            signal
#>                 FALSE      TRUE
#>   carving   0.8815789 0.9010601
#>   naive     0.5000000 0.9222615
#>   splitting 0.9210526 0.9257951