LASSO 回归 - p 值和系数

机器算法验证 p 值 套索
2022-03-12 18:06:54

我已经在 R 中使用cv.glmnet. 我想为选择的系数生成 p 值。

我发现 boot.lass.proj 可以生成引导 p 值 https://rdrr.io/rforge/hdi/man/boot.lasso.proj.html

虽然该boot.lasso.proj程序产生了 p 值,但我认为它正在做自己的套索 - 但我没有看到获得系数的方法。

将 p 值hdi用于由 产生的系数是否安全cv.glmnet

2个回答

为了扩展 Ben Bolker 在对另一个答案的评论中指出的内容,频率论者p值对 LASSO 中的回归系数意味着什么的问题一点也不简单。您测试系数值的实际零假设是什么?您如何考虑 LASSO 对来自同一群体的多个样本执行的这一事实可能会返回完全不同的预测变量集,尤其是在实践中经常看到的相关预测变量类型?您如何考虑您已将结果值用作模型构建过程的一部分,例如在交叉验证或您用于选择惩罚水平以及保留预测变量数量的其他方法中?

这些问题在本网站上进行了讨论。此页面是一个很好的起点,其中包含指向hdi您提到的 R 包的链接以及指向该selectiveInference包的链接,这也在此页面上进行了讨论。稀疏的统计学习在第 6 章中介绍了 LASSO 的推理,并参考了几年前的文献。

请不要简单地将这些或任何其他 LASSO 方法返回的p值用作简单的即插即用结果。重要的是要考虑为什么/是否需要p值以及它们在 LASSO 中的真正含义。如果你的主要兴趣是预测而不是推理,那么预测性能的衡量标准对你和你的听众来说会更有用。

回想一下 LASSO 的作用是消除过程。换句话说,它使用 CV 保持“最佳”特征空间。一种可能的补救措施是选择最终的特征空间并将其反馈到lm命令中。这样,您将能够计算最终选择的 X 变量的统计显着性。例如,请参见以下代码:

library(ISLR)
library(glmnet)
ds <- na.omit(Hitters)
X <- as.matrix(ds[,1:10])
lM_LASSO <- cv.glmnet(X,y = log(ds$Salary),
                intercept=TRUE, alpha=1, nfolds=nrow(ds),
                parallel = T)
opt_lam <- lM_LASSO$lambda.min
lM_LASSO <- glmnet(X,y = log(ds$Salary),
        intercept=TRUE, alpha=1, lambda = opt_lam)
W <- as.matrix(coef(lM_LASSO))
W

                        1
(Intercept)  4.5630727825
AtBat       -0.0021567122
Hits         0.0115095746
HmRun        0.0055676901
Runs         0.0003147141
RBI          0.0001307846
Walks        0.0069978218
Years        0.0485039070
CHits        0.0003636287

keep_X <- rownames(W)[W!=0]
keep_X <- keep_X[!keep_X == "(Intercept)"]
X <- X[,keep_X]
summary(lm(log(ds$Salary)~X))

Call:
lm(formula = log(ds$Salary) ~ X)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.23409 -0.45747  0.06435  0.40762  3.02005 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.5801734  0.1559086  29.377  < 2e-16 ***
XAtBat      -0.0025470  0.0010447  -2.438  0.01546 *  
XHits        0.0126216  0.0039645   3.184  0.00164 ** 
XHmRun       0.0057538  0.0103619   0.555  0.57919    
XRuns        0.0003510  0.0048428   0.072  0.94228    
XRBI         0.0002455  0.0045771   0.054  0.95727    
XWalks       0.0072372  0.0026936   2.687  0.00769 ** 
XYears       0.0487293  0.0206030   2.365  0.01877 *  
XCHits       0.0003622  0.0001564   2.316  0.02138 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6251 on 254 degrees of freedom
Multiple R-squared:  0.5209,    Adjusted R-squared:  0.5058 
F-statistic: 34.52 on 8 and 254 DF,  p-value: < 2.2e-16

请注意,这些系数与从模型中得出的系数几乎没有什么不同glmnet最后,您可以使用该stargazer包输出到格式良好的表格中。在这种情况下,我们有

stargazer::stargazer(lm(log(ds$Salary)~X),type = "text")
===============================================
                        Dependent variable:    
                    ---------------------------
                              Salary)          
-----------------------------------------------
XAtBat                       -0.003**          
                              (0.001)          

XHits                        0.013***          
                              (0.004)          

XHmRun                         0.006           
                              (0.010)          

XRuns                         0.0004           
                              (0.005)          

XRBI                          0.0002           
                              (0.005)          

XWalks                       0.007***          
                              (0.003)          

XYears                        0.049**          
                              (0.021)          

XCHits                       0.0004**          
                             (0.0002)          

Constant                     4.580***          
                              (0.156)          

-----------------------------------------------
Observations                    263            
R2                             0.521           
Adjusted R2                    0.506           
Residual Std. Error      0.625 (df = 254)      
F Statistic           34.521*** (df = 8; 254)  
===============================================
Note:               *p<0.1; **p<0.05; ***p<0.01

引导程序

使用引导方法,我将上述标准错误与引导错误进行比较,作为稳健性检查:

library(boot)

W_boot <- function(ds, indices) {
  ds_boot <- ds[indices,] 
  X <- as.matrix(ds_boot[,1:10])
  y <- log(ds$Salary)
  lM_LASSO <- glmnet(X,y = log(ds$Salary),
                     intercept=TRUE, alpha=1, lambda = opt_lam)
  W <- as.matrix(coef(lM_LASSO))
  return(W)
} 

results <- boot(data=ds, statistic=W_boot, 
                R=10000)

se1 <- summary(lm(log(ds$Salary)~X))$coef[,2]
se2 <- apply(results$t,2,sd)
se2 <- se2[W!=0]
plot(se2~se1)
abline(a=0,b=1)

在此处输入图像描述 截距似乎有一个小的偏差。否则,临时方法似乎是合理的。无论如何,您可能想查看此线程以进一步讨论此问题。