The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.
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.
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.
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.2417481We 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 14The 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-01index 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 TRUEsummary() 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.448890562The 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.5124690The randomized intervals use the information that selection left over, so they are typically shorter than data splitting while keeping the nominal coverage.
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.