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(100) 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(100) 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(100) 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.
cv.glmnet()to choose the penalty
lambdathat minimizes cross-validated error. $\endgroup$
glmnetfunction by default chooses a set of 100 lambda values to evaluate, while the others compute the entire sequence of fits at once if
p<n(as you seem to have in your data). If your pre-ordained choice of
s=.1is between two of the
glmnetchoices and is close to the threshold for inclusion/exclusion of a predictor, you could find
glmnetincluding or excluding predictors relative to the others. $\endgroup$