I’m not really sure if a package or separate function that does what you want exists. I happen to have tried to do something like this myself some time ago and ran into the same problem. I then wrote a very simple bootstrap function that performs the original fit and the bootstrap for it (relies on the `glmnet`

package). I won’t post all the code here as it’s quite ugly and written very specifically for my own dataset and needs at that moment, but the part that you would need looked somewhat like this:

```
# Requires as input a prediction and output matrix (pred.mat and out.mat) as well as a lamba grid over which to use the glmnet functions
lasso_bootstrapper <- function(pred.mat, out.mat, lambda.grid) {
# Calculates lamba.min (which I used in my analysis)
lasso_lamda_min <- cv.glmnet(x = pred.mat, y = out.mat, lambda = lambda.grid, alpha = 1)$lambda.min
# Performs the original fit on the complete data.
# Note: I was working with a continous outcome, hence family = "gaussian". Change this if necessary
orig.fit <- glmnet(x = pred.mat, y = out.mat, family = "gaussian", lambda = lambda.grid, alpha = 1)
# Extracts the coefficients from the fit on the full dataset
coef.orig <- predict(orig.fit, s = lasso_lamda_min, newx = pred.mat, type = "coefficients")
# Dataframe for storing results from the bootstrap. Change ncol according to the number of predictors/variables in your analysis
boot.res <- data.frame(matrix(nrow = B, ncol = nrow(coef.orig))
for(i in 1:B) {
# Create indicators for sub-sampling the original dataframe
boot.train <- sample(n.orig, replace = T)
# Calculate lamba.min for the bootstrap sample
lasso_lamda_boot_min <- cv.glmnet(x = pred.mat[boot.train], y = out.mat[boot.train], lambda = lambda.grid, alpha = 1)$lambda.min
# Fits the lasso model in the bootstrap sample and extracts coefficients
boot.fit <- glmnet(x = pred.mat[boot.train,], y = out.mat[boot.train], family = "gaussian", lambda = lambda.grid, alpha = 1)
coef.boot <- predict(boot.fit, s = lasso_lamda_boot_min , type = "coefficients")
# Stores coefficients from the fit in the bootstrap sample in boot.res
boot.res[i,] <- coef.boot
}
# Returns the boot.res as output for the function
return(boot.res)
}
```

This will result in a dataframe of all the coefficients in your bootstrap sample. If a variable was not selected it results in a zero if I recall correctly, so using `apply`

you can quickly summarize the selection frequency.

In reply to @f2harrell below: we used it in a context where we first had a well specified basic model and then attempted to test if we could find any improvement when adding in (a large number of) predictors from a -omics platform. The eventual decision about their usefullness was based on actual increments in predictive performance, in large part inspired by your blogpost. The added performance was (almost) nothing and this was just one small analysis to get a better idea of what was going on in the background.