LASSO Regression - p-values and coefficients

Question

I've run a LASSO in R using cv.glmnet. I would like to generate p-values for the coefficients that are selected.

I found the boot.lass.proj to produce bootstrapped p-values https://rdrr.io/rforge/hdi/man/boot.lasso.proj.html

While the boot.lasso.proj program produced p-values, I assume it is doing its own lasso - but I'm not seeing a way to get the coefficients.

Would it be safe to use the p-values from hdi for the coefficients produced by cv.glmnet?

Answer 1

To expand on what Ben Bolker notes in a comment on another answer, the issue of what a frequentist p-value means for a regression coefficient in LASSO is not at all easy. What's the actual null hypothesis against which you are testing the coefficient values? How do you take into account the fact that LASSO performed on multiple samples from the same population may return wholly different sets of predictors, particularly with the types of correlated predictors that often are seen in practice? How do you take into account that you have used the outcome values as part of the model-building process, for example in the cross-validation or other method you used to select the level of penalty and thus the number of retained predictors?

These issues are discussed on this site. This page is one good place to start, with links to the R hdi package that you mention and also to the selectiveInference package, which is also discussed on this page. Statistical Learning with Sparsity covers inference for LASSO in Chapter 6, with references to the literature as of a few years ago.

Please don't simply use the p-values returned by those or any other methods for LASSO as simple plug-and-play results. It's important to think why/whether you need p-values and what they really mean in LASSO. If your main interest is in prediction rather than inference, measures of predictive performance would be much more useful to you and to your audience.

Answer 2

Recall that LASSO functions as an elimination process. In other words, it keeps the "best" feature space using CV. One possible remedy is to select the final feature space and feed it back into an lm command. This way, you would be able to compute the statistical significance of the final selected X variables. For instance, see the following code:

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

Note that the coefficients are little different from the ones derived from the glmnet model. Finally, you can use the stargazer package to output into a well-formatted table. In this case, we have

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

Bootstrap

Using a bootstrap approach, I compare the above standard errors with the bootstrapped one as a robustness check:

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)

There seems to be a small bias for the intercept. Otherwise, the ad-hoc approach seems to be justified. In any case, you may wanna check this thread for further discussion on this.