6
$\begingroup$

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?

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

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Digio May 27 at 11:50
  • $\begingroup$ 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. $\endgroup$ – jpryan28 May 27 at 12:16
  • 1
    $\begingroup$ @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. $\endgroup$ – EdM May 27 at 16:18
  • $\begingroup$ @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. $\endgroup$ – EdM May 27 at 16:27
  • $\begingroup$ @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? $\endgroup$ – Digio May 27 at 18:47
1
$\begingroup$

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)

enter image description here 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.

$\endgroup$
  • 6
    $\begingroup$ 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 ...) $\endgroup$ – Ben Bolker May 26 at 20:08
  • 1
    $\begingroup$ 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. $\endgroup$ – majeed simaan May 26 at 22:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.