Randomized LASSO for a Gaussian linear model

This vignette shows the basic workflow of lassoinf. We fit a LASSO to a randomized response, then use LassoInference for confidence intervals and p-values for the selected coefficients that are valid after selection.

Randomization and the proportion used for selection

In data splitting, a proportion $\pi$ of the data is used to select the model and the rest for inference. For the score $Z = X'y$, with variance $\Sigma = \sigma^2 X'X$, the score of a $\pi$ fraction of the data has variance $\Sigma / \pi$.

Randomization instead selects with $Z_{noisy} = Z + \omega$, $\omega \sim N(0, \Omega)$. Matching $\text{Var}(Z_{noisy}) = \Sigma + \Omega = \Sigma / \pi$ gives

\[ \Omega = \frac{1 - \pi}{\pi} \Sigma. \]

For least squares, this is the same as fitting the LASSO to $y_{noisy} = y + \epsilon$ with $\epsilon \sim N(0, \frac{1-\pi}{\pi}\sigma^2 I)$. As $\pi \to 1$ almost all the information is used for selection and little is left for inference. Smaller $\pi$ leaves more.

Data

We use $n = 100$ and $p = 20$, with one strong effect, two weaker ones and $\sigma^2 = 2$. Selection uses $\pi = 0.9$.

library(lassoinf)
set.seed(0)

n <- 100; p <- 20
signal_strength <- 4
X <- matrix(rnorm(n * p), n, p)
beta <- numeric(p)
beta[1:3] <- c(2, 1, -1) * signal_strength / sqrt(n)
sigma2 <- 2
y <- drop(X %*% beta) + rnorm(n, sd = sqrt(sigma2))

With $n > p$ we estimate $\sigma^2$ from the full least squares fit, and build $\Sigma$ and $\Omega$:

pi_selection <- 0.9
Q <- crossprod(X)
Z_full <- drop(crossprod(X, y))
sigma2_hat <- sum(lm.fit(X, y)$residuals^2) / (n - p)
Sigma <- sigma2_hat * Q
sigma2_rand <- (1 - pi_selection) / pi_selection * sigma2_hat
Sigma_noise <- sigma2_rand * Q
y_noisy <- y + rnorm(n, sd = sqrt(sigma2_rand))
c(sigma2_hat = sigma2_hat, sigma2_rand = sigma2_rand)
#>  sigma2_hat sigma2_rand 
#>   2.1757329   0.2417481

Fit the LASSO on the randomized response

We fit the LASSO $\frac{1}{2}\|y_{noisy} - Xb\|^2 + \lambda \|b\|_1$. glmnet scales the loss by $1/n$, so we pass it $\lambda / n$. LassoInference needs the solution, the gradient of the smooth part at it, $G = X'(X\hat{\beta} - y_{noisy})$, and its Hessian, $Q = X'X$.

lam <- 2 * sqrt(n)
fit <- glmnet::glmnet(X, y_noisy, lambda = lam / n, intercept = FALSE,
                      standardize = FALSE, control = list(thresh = 1e-14))
beta_hat <- as.numeric(coef(fit))[-1]
G_hat <- drop(crossprod(X, X %*% beta_hat - y_noisy))
which(beta_hat != 0)
#> [1]  1  3 10 13 14

Inference

The arguments are the solution, gradient and Hessian, the penalty weights $D$, the bounds $L$, $U$ (none here), the full-data score and its covariance, and the randomization covariance.

li <- LassoInference$new(beta_hat, G_hat, Q, rep(lam, p), NULL, NULL,
                         Z_full, Sigma, Sigma_noise, level = 0.9)
summ <- li$summary()
summ
#>    index    beta_hat  lower_conf  upper_conf      p_value
#> 0      0  0.41432718  0.50799234  1.08675205 7.329526e-05
#> 2      2 -0.03783531 -0.57080194  0.04672353 1.456693e-01
#> 9      9 -0.16256389 -0.65087642 -0.01061865 9.066102e-02
#> 12    12 -0.02702932 -0.02090107  1.03813586 1.129303e-01
#> 13    13  0.13122199 -0.49668608  0.46893379 8.089657e-01

index is 0-based. beta_hat is the LASSO solution, which is shrunk toward zero and can lie outside its interval. The intervals are for the coefficients of the selected model, $\beta_E = (X_E'X_E)^{-1} X_E' \mu$ with $\mu = E[y]$. They come from the contrasts in li$contrasts, so we can compute the targets and check coverage:

Z_mu <- drop(Q %*% beta)
summ$truth <- sapply(as.character(summ$index), function(j) sum(li$contrasts[[j]] * Z_mu))
summ$cover <- summ$lower_conf <= summ$truth & summ$truth <= summ$upper_conf
summ
#>    index    beta_hat  lower_conf  upper_conf      p_value        truth cover
#> 0      0  0.41432718  0.50799234  1.08675205 7.329526e-05  0.855889010  TRUE
#> 2      2 -0.03783531 -0.57080194  0.04672353 1.456693e-01 -0.406730210  TRUE
#> 9      9 -0.16256389 -0.65087642 -0.01061865 9.066102e-02  0.009152833 FALSE
#> 12    12 -0.02702932 -0.02090107  1.03813586 1.129303e-01  0.006070596  TRUE
#> 13    13  0.13122199 -0.49668608  0.46893379 8.089657e-01 -0.040488544  TRUE

summary() also stores naive intervals, which ignore selection, and data-splitting intervals in li$"_naive" and li$"_splitting":

li$`_naive`
#>    index   beta_hat  lower_conf  upper_conf      p_value
#> 0      0  0.8038264  0.52088364  1.08676909 2.969066e-06
#> 2      2 -0.3414098 -0.57683903 -0.10598055 1.706525e-02
#> 9      9 -0.3485781 -0.57251641 -0.12463971 1.045677e-02
#> 12    12 -0.1524278 -0.39815843  0.09330274 3.075812e-01
#> 13    13  0.2684742  0.01106352  0.52588484 8.624453e-02
li$`_splitting`
#>    index   beta_hat lower_conf upper_conf     p_value
#> 0      0  1.4619576  0.5672141  2.3567010 0.007196907
#> 2      2 -1.4585552 -2.2030478 -0.7140625 0.001270877
#> 9      9 -0.4672224 -1.1753776  0.2409328 0.277819090
#> 12    12  0.6308092 -0.1462591  1.4078776 0.181791012
#> 13    13 -0.3747552 -1.1887592  0.4392488 0.448890562

Coverage over repeated experiments

The randomized intervals should cover their targets 90% of the time. The naive intervals should not.

one_replication <- function() {
  X <- matrix(rnorm(n * p), n, p)
  y <- drop(X %*% beta) + rnorm(n, sd = sqrt(sigma2))
  Q <- crossprod(X)
  sigma2_hat <- sum(lm.fit(X, y)$residuals^2) / (n - p)
  sigma2_rand <- (1 - pi_selection) / pi_selection * sigma2_hat
  y_noisy <- y + rnorm(n, sd = sqrt(sigma2_rand))
  fit <- glmnet::glmnet(X, y_noisy, lambda = lam / n, intercept = FALSE,
                        standardize = FALSE, control = list(thresh = 1e-14))
  beta_hat <- as.numeric(coef(fit))[-1]
  if (all(beta_hat == 0)) return(NULL)
  G_hat <- drop(crossprod(X, X %*% beta_hat - y_noisy))
  li <- LassoInference$new(beta_hat, G_hat, Q, rep(lam, p), NULL, NULL,
                           drop(crossprod(X, y)), sigma2_hat * Q, sigma2_rand * Q, level = 0.9)
  summ <- li$summary()
  truth <- sapply(as.character(summ$index), function(j) sum(li$contrasts[[j]] * drop(Q %*% beta)))
  covers <- function(df) df$lower_conf <= truth & truth <= df$upper_conf
  data.frame(method = rep(c("randomized", "naive", "splitting"), each = nrow(summ)),
             cover = c(covers(summ), covers(li$`_naive`), covers(li$`_splitting`)),
             length = c(summ$upper_conf - summ$lower_conf,
                        li$`_naive`$upper_conf - li$`_naive`$lower_conf,
                        li$`_splitting`$upper_conf - li$`_splitting`$lower_conf))
}

set.seed(1)
results <- do.call(rbind, replicate(200, one_replication(), simplify = FALSE))
data.frame(coverage = tapply(results$cover, results$method, mean),
           median_length = tapply(results$length, results$method, median))
#>             coverage median_length
#> naive      0.7381166     0.4782847
#> randomized 0.8914798     0.8712314
#> splitting  0.9013453     1.5124690

The randomized intervals use the information that selection left over, so they are typically shorter than data splitting while keeping the nominal coverage.