LASSO Regression - p-values and coefficients

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?

• What is hdi here? – Richard Hardy May 26 '19 at 18:16
• HDI - the hdi package that contains boot.lass.proj rdrr.io/rforge/hdi – jpryan28 May 26 '19 at 22:29

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.

• My question is, if you somehow derive empirical p-values (via Bootstrap or random perturbations or otherwise), won't it be OK to trust them then? – Digio May 27 '19 at 11:50
• Thanks, EdM. My main concern is that when I give the results to my collaborators (or reviewers) the first question is always "but is it statistically significant!?" I think I can give them a bit more explanation now of why that isn't really the point. – jpryan28 May 27 '19 at 12:16
• @jpryan28 I don't mean to discourage you from trying these methods for p-value calculations in LASSO. Just wanted to make sure that you think carefully about the issues. Your practical concerns with your audiences are certainly valid. It might be a good strategy to have such values in hand, but then present them (if required) in a way that demonstrates your deep understanding of the statistical issues. This page, also linked from another answer, is a further helpful resource. – EdM May 27 '19 at 16:18
• @Digio LASSO models based on different bootstrap samples will typically differ in terms of the predictors that are chosen, particularly with the correlated predictors typical in practice. If $X_1$ and $X_2$ are correlated predictors, each of which shows up exclusive of the other in half of the bootstrap models, what do you mean by a confidence interval for the coefficient of either of them? Note that this is not such a problem for predictions made from the models, as in this situation $X_1$ and $X_2$ serve as proxies for each other. – EdM May 27 '19 at 16:27
• @EdM What I'm asking is not specific to bootstrap nor Lasso. The way I see it, the problem has to do with the traditional way of calculating p-values and the fact that any model selection method makes them biased. Therefore, if you choose to use an empirical method instead of Wald test to derive p-values, having applied a model selection technique (lasso or otherwise) should not be posing a problem anymore? – Digio May 27 '19 at 18:47

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
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.

• I have very serious concerns about this approach (I know it's widely done), because it doesn't take the selection process (nor the process of selecting the penalty by cross-validation) into account. Can you point to references that assert that this approach actually has correct frequentist properties? (I can imagine that it's OK in particular cases where you have lots of data and the parameters have either zero (or very small) effects or very large effects ...) – Ben Bolker May 26 '19 at 20:08
• Good point. I don't have a specific paper to justify this methodology, but my thinking that this should not be that different from the "true" standard errors. As you mentioned above, running a bootstrap would serve as a more robust solution. I reran the code with the addition of bootstrap and compared the results with the lm standard errors. Ignoring the intercept, the results seem to be consistent. – majeed simaan May 26 '19 at 22:00
• I thought I would reference "An Introduction to Statistical Learning" (Gareth James, Daniela Witten, Trevor Hastie, Robert Tibshirani) to expand on the way "best" is used above for high dimensions: "In the high-dimensional setting, ... we can never know exactly which variables (if any) truly are predictive of the outcome, and we can never identify the best coefficients for use in the regression. At most, we can hope to assign large regression coefficients to variables that are correlated with the variables that truly are predictive of the outcome." – user271536 Aug 27 '20 at 19:53
• @majeedsimaan Should one not re-estimate the optimal value of lambda in each bootstrap replication? Furthermore, the size of standard errors should be judged in comparison to the size of the estimated coefficient. I.e., is the bias still small if we would zoom in on the bottom-left part of the plot? Also, the df value of the lm need adjusting, the open question obviously being by how much. – Marjolein Fokkema Feb 26 at 2:01