R：如何将 cv.glmnet() 函数中的 lasso lambda 值转换为 selectionInference 包？

Posted 2023-02-14

【中文标题】R：如何将 cv.glmnet() 函数中的 lasso lambda 值转换为 selectionInference 包？

【英文标题】：R: How to translate lasso lambda values from the cv.glmnet() function into the selectiveInference package?

【发布时间】：2021-09-09 19:11:04

【问题描述】：

我正在使用selectiveInference 包使用套索（“l1 norm”）执行选择后推理。这个包假定 lambda 是固定的——也就是说，我们事先确定了它。但是，我需要使用交叉验证。

Taylor & Tibshirani (2018) 使用模拟表明使用交叉验证来确定 lambda 会产生有效的推理统计，使用 selectiveInference::fixedLassoInf() 方法。 （Another paper 提出了一种处理由交叉验证确定的 lambda 的方法，但它似乎还没有包含在包中，并且 2018 年论文中的模拟对我来说已经足够好。）

我在文档中看到它说glmnet 使用 1/n 套索参数化，而 selectiveInference 使用通用参数化。文档展示了如何从普通 lambda 转换为 glmnet 可以使用的东西。

我需要做相反的事情：从 cv.glmnet() 给我的东西出发，然后将它变成 fixedLassoInf() 想要的通用尺度的 lambda。

具体来说，glmnet 文档内容如下：

另请注意，对于“高斯”，glmnet 在计算其 lambda 序列之前将 y 标准化为具有单位方差（使用 1/n 而不是 1/(n-1) 公式）（然后对结果系数进行非标准化）；如果您想用其他软件重现/比较结果，最好提供一个标准化的 y

虽然selectiveInference 说：

估计的套索系数（例如，来自 glmnet）。这是长度 p （因此截距不包括在第一个组件中）。当心！此函数使用“标准”套索目标......相比之下，glmnet 将第一项乘以 1/n 的因子。所以运行glmnet后，要提取一个值lambda对应的beta，需要使用beta = coef(obj,s=lambda/n)[-1]...

有关可重现的示例，请参见下面的代码。

我的问题特别关注如何调整这条线：si_lambda <- glmnet_lambda。也就是说，我要做什么转换从 lambda cv.glmnet() 给我（我将其分配给 glmnet_lambda）到 selectiveInference 将使用的 lambda（我称之为 si_lambda )?

我最初的想法是，由于文档说要除以 n，我的想法是将 cv.glmnet() 给我的值乘以我的样本量。运行时不会发出警告或错误，但它给了我一个 188.5121 的 lambda，感觉不对。抱歉，如果这是答案并且我只是很密集 - 但我想确保我以适当的方式从一个软件转移到另一个软件。

library(glmnet)
library(selectiveInference)
library(tidyverse)
set.seed(1839)

n <- 1000       # sample size
B <- c(0, 1, 0) # intercept 0, beta1 = 1, beta2 = 0
eps_sd <- 1     # sd of the error

# make data
X <- cbind(1, replicate(length(B) - 1, rnorm(n, 0, 1)))
y <- X %*% B + rnorm(n, 0, eps_sd)
dat <- as.data.frame(X[, -1])
dat <- as_tibble(cbind(dat, y))

# get lambda by way of cross-validation
glmnet_lambda <- cv.glmnet(
  x = as.matrix(select(dat, -y)),
  y = dat$y
) %>% 
  getElement("lambda.1se")

# run glmnet with that lambda
m1 <- glmnet(
  x = as.matrix(select(dat, -y)),
  y = dat$y,
  lambda = glmnet_lambda
)

# get coefs from that model, dropping intercept, per the docs
m1_coefs <- coef(m1)[-1]

# what reparameterization do I do here?
si_lambda <- glmnet_lambda

# do post-selection inference with m1
# runs with warning, so I assume parameterized incorrectly -- how to fix?
m2 <- fixedLassoInf(
  x = as.matrix(select(dat, -y)),
  y = dat$y,
  beta = m1_coefs,
  lambda = si_lambda
)

以及会话信息：

> sessionInfo()
R version 4.1.0 (2021-05-18)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 11.4

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] forcats_0.5.1            stringr_1.4.0            dplyr_1.0.6             
 [4] purrr_0.3.4              readr_1.4.0              tidyr_1.1.3             
 [7] tibble_3.1.2             ggplot2_3.3.3            tidyverse_1.3.1         
[10] selectiveInference_1.2.5 MASS_7.3-54              adaptMCMC_1.4           
[13] coda_0.19-4              survival_3.2-11          intervals_0.15.2        
[16] glmnet_4.1-1             Matrix_1.3-3            

【参考方案1】：

需要翻一下fixedLassoInf文档中的例子；将其调整为您的示例将授予以下代码：

library(glmnet)
library(selectiveInference)

# Make dataset
set.seed(1839)
n <- 1000       # sample size
B <- c(0, 1, 0) # intercept 0, beta1 = 1, beta2 = 0
eps_sd <- 1     # sd of the error
X <- cbind(1, replicate(length(B) - 1, rnorm(n, 0, 1)))
y <- X %*% B + rnorm(n, 0, eps_sd)

# Cross-validation to find lambda
gfit = cv.glmnet(X[,-1], y) # we need to remove the intercept variable (glmnet will add another one)
lambda = gfit$lambda.min

# Obtain coefficients (properly scaling lambda and removing the intercept coefficient)
(beta = coef(gfit, x=X[,-1], y=y, s=lambda, exact=TRUE)[-1])
# [1]  0.99297607 -0.04300646

# Compute fixed lambda p-values and selection intervals
(out = fixedLassoInf(X[,-1], y, beta, lambda*n))
# Call:fixedLassoInf(x = X[, -1], y = y, beta = beta, lambda = lambda * n)
# 
# Standard deviation of noise (specified or estimated) sigma = 1.012
# 
# Testing results at lambda = 4.562, with alpha = 0.100
# 
# Var   Coef Z-score P-value LowConfPt UpConfPt LowTailArea UpTailArea
# 1  0.998  31.475   0.000     0.945    1.050       0.049      0.049
# 2 -0.048  -1.496   0.152    -0.100    0.032       0.050      0.049
# 
# Note: coefficients shown are partial regression coefficients

【讨论】：

文档涉及在选择推理公式中获取给定的 lambda，并将其转换为套索的 glmnet 公式。我对相反的事情感兴趣：所以我对coef() 的调用将涉及到s = lambda，然后必须以某种方式将lambda 转换为selectiveInference 公式。文档假设我们有一个 lambda，需要将其转换为 glmnet——我有一个来自 glmnet 的 lambda 值，需要将其转换为选择性推理。

你是对的。我现在已经更正了答案。

此外，有用户在选择性推理方面存在问题的历史记录 (stats.stackexchange.com/questions/tagged/selectiveinference)

字，所以你读到的，和我的一样，是乘以n而不是除？它有时会给我奇怪的 lambda 值。是的，我确保我正确执行此操作的原因之一是我正在做一些模拟研究，以查看包是否以有效的误报率等执行，我遇到了奇怪的问题做乘以 n 的事情。

是的，确实，未缩放的 lambda 看起来相当大，但这只是因为您已经习惯了 glmnet 处理的比例（它使用自己的标准来定义 lambda，因此会产生很多混乱）。请务必阅读有关 fixedLassoInf 的详细信息；您的模拟可能无法测量此处实际计算的内容（例如，它使用部分回归系数）。对 lasso 回归进行推断并不简单！