Different sets of features selected by three different functions in R for running LASSO Regressions despite the same random seed for each

Question

The GitHub Repository for this research project has all of the code included in this question. Brief background context: I am just finishing up the work on my part as a coauthor on a research project exploring the properties of a newly proposed Automated Optimal Variable Selection algorithm. My role is to decide which existing automated feature selection techniques to use as the benchmarks, then run them and evaluate their performance using standard classification performance metrics in R. I ended up choosing 3 benchmarks: Backward Elimination Stepwise Regression, Forward Selection Stepwise Regression, and LASSO Regression.

I have already done this, but in the process I ran into something most unexpected which I have not been able to find any journal articles, R documentation, or textbook sub sections to explain. My collaborator randomly generated 260,000 synthetic datasets for me to run the benchmarks on, so I ran 260k LASSO Regressions (initially using the enet() function from R's elastic net package), one for each dataset, got my results, and calculated my performance metrics. And then, I did so again using the glmnet() function from the glmnet package and found that the set of variables it selected in each dataset was not always identical to the set of variables selected by enet on the same dataset even when using the same random seed value for both.

So, at that point, being exasperated, I threw my hands up and I re-did all again for a third time using the lars() function from the package of the same name and once again, the variables selected were slightly different. In the aforementioned Repository on GitHub, in the Stage 2 Results folder, the fact that each of these selected different variables can be verified by inspecting LASSO's Selections via glmnet.xlsx, LASSO's Selections via lars.xlsx, Variables Selected by LASSO ran via enet.xlsx, or Overall LASSO Performance Metrics.xlsx.

To take an example at random per a helpful suggestion below in the comments, for dataset 0.25-3-1-1, these are each of their sets of selected factors respectively:

• enet: X5, X21, X22
• lars: X21, X22
• glment: X5, X21, X22, X30

And for completeness with respect to this example, the correct set of regressors, i.e. the structural variables for the 0.25-3-1-1 dataset is: X5, X21, X22

Which means only enet got the exact right answer (a True Positive Rate of 1 and a True Negative Rate of 1).

Here is the code I used to run them via the enet function:

set.seed(11)     
enet_LASSO_fits <- lapply(datasets, function(i) 
               elasticnet::enet(x = as.matrix(dplyr::select(i, 
                                         starts_with("X"))), 
               y = i$Y, lambda = 0, normalize = FALSE))

# This separates out and stores just the coefficient estimates from each LASSO.
LASSO_Coeffs <- lapply(enet_LASSO_fits, 
                          function(i) predict(i, 
                                           x = as.matrix(dplyr::select(i, starts_with("X"))), 
                                           s = 0.1, mode = "fraction", 
                                           type = "coefficients")[["coefficients"]])

# Write my own custom lapply which will separate out and return a 
# new list containing just the Independent Variables
# which are 'selected' or chosen for each individual dataset.
IVs_Selected <- lapply(LASSO_Coeffs, function(i) names(i[i > 0]))

Here is the code for running them via the glmnet function:

set.seed(11)     
glmnet_lasso.fits <- lapply(datasets, function(i) 
               glmnet(x = as.matrix(select(i, starts_with("X"))), 
                      y = i$Y, alpha = 1))

# This stores and prints out all of the regression 
# equation specifications selected by LASSO when called
lasso.coeffs = glmnet_lasso.fits |> 
  Map(f = \(model) coef(model, s = .1))

Variables.Selected <- lasso.coeffs |>
  Map(f = \(matr) matr |> as.matrix() |> 
       as.data.frame() |> filter(s1 != 0) |> rownames())
Variables.Selected = lapply(seq_along(datasets), \(j)
                            j <- (Variables.Selected[[j]][-1]))

Here is the code for running them via the lars function:

set.seed(11)     
lars_LASSO_fits <- lapply(datasets, function(i) 
  lars(x = as.matrix(select(i, starts_with("X"))), 
         y = i$Y, type = "lasso"))

# This stores and prints out all of the regression 
# equation specifications selected by LASSO when called
Lars.Coeffs <- lapply(lars_LASSO_fits, 
                       function(i) predict(i, 
                                           x = as.matrix(dplyr::select(i, starts_with("X"))), 
                                           s = 0.1, mode = "fraction", 
                                           type = "coefficients")[["coefficients"]])

IVs.Selected.by.Lars <- lapply(Lars.Coeffs, function(i) names(i[i > 0]))

What is going on here? Am I doing something wrong or do each of these fitting functions use different underlying stopping conditions or starting conditions or something like that?

p.s. 1 - The script I used to run my 260k LASSOs for the first time (by way of enet()) is the one in the GitHub Repo called "LASSO Regressions.R", the script in which I estimated them using the glmnet function is fittingly called "LASSO using the 'glmnet' package.R", and the one which in used lars to fit them is called "LASSO using Lars.R".

p.s. 2 - By the way, I had re-ran all of my 260k Backward and Forward Stepwise Regressions using the stepAIC() function from R's MASS package instead what I used the first time around, namely, the step() function from the stats library and all of the variables it selected were identical in every case, and as a result of that, I had no doubts that this would be the same for LASSO.

p.s. As was suggested in a comment thread with the person who answered this question, I asked a follow up question to this one afterwards.

Answer 1

First, as a comment suggested, none of these functions uses random sampling. The data sets I examined all had the number of candidate parameters less than the number of observations. If that's the case, you can get full linear least-squares (LS) regression fits as a reference.

With lars() and the enet() function that's based on it, the functions can return the entire path of solutions from fully penalized (all coefficients 0) to that LS solution. That's explained in Chapter 3 of Elements of Statistical Learning. When you make predictions from such models with a penalty-factor choice s, you must specify what that penalty means.

You chose s = 0.1 and mode = "fraction" for those models. According to the lars manual, with that mode:

If mode="fraction", then s should be a number between 0 and 1, and it refers to the ratio of the L1 norm of the coefficient vector, relative to the norm at the full LS solution.

So in those cases you are (arbitrarily) choosing the set of coefficients along the path that has a sum of absolute values equal to one-tenth of the sum of coefficient absolute values at the LS solution. I think that you should get the same results from both types of models. You specified normalize=FALSE in the call to enet() but didn't to lars(), which might be responsible (perhaps due to numerical precision issues) for the difference you found even if the predictor variables were nominally standardized to start with.

The penalty factor in glmnet() is different. With its more general application, that function chooses by default 100 penalty parameter values. The maximum is the smallest penalty that sets all coefficients to 0, and the minimum is a small fraction of it; "the default depends on the sample size nobs relative to the number of variables nvars" according to the manual.

Furthermore, the value of s means something completely different from what you specified for the other functions. Again, quoting from the manual:

Value(s) of the penalty parameter lambda at which predictions are required.

That's the actual penalty applied to the L1 norm of the coefficients, not a fraction of the L1 norm at the LS solution as you specified to the other functions. The value used by predict.glmnet() is equivalent to mode = "penalty" for predict.enet() and mode = "lambda" for predict.lars().

Finally, I respectfully suggest that this approach with a fixed choice of s (even when you eventually specify the same type of s to all models) is not consistent with your intent to use LASSO as a "benchmark" against which to compare another predictor-selection method. Choosing an arbitrary penalty factor is not how these functions are used in practice. Typically, one would use cv.glmnet() to find an optimal penalty factor, and return the coefficients that were retained in the model at that optimal penalty.* I think you will find it very hard to convince reviewers that your approach to feature selection with LASSO, with a single fixed value of s, constitutes valid "benchmarking."

---

*Unlike your code, that would require random sampling and a corresponding set.seed() for reproducibility.