23
$\begingroup$

I'm reading about best subset selection in the Elements of statistical learning book. If I have 3 predictors $x_1,x_2,x_3$, I create $2^3=8$ subsets:

  1. Subset with no predictors
  2. subset with predictor $x_1$
  3. subset with predictor $x_2$
  4. subset with predictor $x_3$
  5. subset with predictors $x_1,x_2$
  6. subset with predictors $x_1,x_3$
  7. subset with predictors $x_2,x_3$
  8. subset with predictors $x_1,x_2,x_3$

Then I test all these models on the test data to choose the best one.

Now my question is why is best subset selection not favored in comparison to e.g. lasso?

If I compare the thresholding functions of best subset and lasso, I see that the best subset sets some of the coefficients to zero, like lasso. But, the other coefficient (non-zero ones) will still have the ols values, they will be unbiasd. Whereas, in lasso some of the coefficients will be zero and the others (non-zero ones) will have some bias. The figure below shows it better: enter image description here

From the picture the part of the red line in the best subset case is laying onto the gray one. The other part is laying in the x-axis where some of the coefficients are zero. The gray line defines the unbiased solutions. In lasso, some bias is introduced by $\lambda$. From this figure I see that best subset is better than lasso! What are the disadvantages of using best subset?

$\endgroup$
3
  • 1
    $\begingroup$ .. and what do the curves look like when the randomness in the data causes you to select one of the many wrong subsets and the associated coefficient estimates are far from zero relative to their standard errors? $\endgroup$
    – jbowman
    Jun 9, 2018 at 13:57
  • 3
    $\begingroup$ @jbowman I don't understand it very clearly, why would the randomness in the data cause me to select the wrong one? If I would use cross validation to select the best subset, I would then have smaller chances to select the wrong subset. $\endgroup$
    – Ville
    Jun 9, 2018 at 14:11
  • 1
    $\begingroup$ You seem to be equating "less bias" with "better". What brings you to place such a high value on unbiasedness? $\endgroup$ Jun 19, 2018 at 16:47

2 Answers 2

23
$\begingroup$

In subset selection, the nonzero parameters will only be unbiased if you have chosen a superset of the correct model, i.e., if you have removed only predictors whose true coefficient values are zero. If your selection procedure led you to exclude a predictor with a true nonzero coefficient, all coefficient estimates will be biased. This defeats your argument if you will agree that selection is typically not perfect.

Thus, to make "sure" of an unbiased model estimate, you should err on the side of including more, or even all potentially relevant predictors. That is, you should not select at all.

Why is this a bad idea? Because of the bias-variance tradeoff. Yes, your large model will be unbiased, but it will have a large variance, and the variance will dominate the prediction (or other) error.

Therefore, it is better to accept that parameter estimates will be biased but have lower variance (regularization), rather than hope that our subset selection has only removed true zero parameters so we have an unbiased model with larger variance.

Since you write that you assess both approaches using cross-validation, this mitigates some of the concerns above. One remaining issue for Best Subset remains: it constrains some parameters to be exactly zero and lets the others float freely. So there is a discontinuity in the estimate, which isn't there if we tweak the lasso $\lambda$ beyond a point $\lambda_0$ where a predictor $p$ is included or excluded. Suppose that cross-validation outputs an "optimal" $\lambda$ that is close to $\lambda_0$, so we are essentially unsure whether p should be included or not. In this case, I would argue that it makes more sense to constrain the parameter estimate $\hat{\beta}_p$ via the lasso to a small (absolute) value, rather than either completely exclude it, $\hat{\beta}_p=0$, or let it float freely, $\hat{\beta}_p=\hat{\beta}_p^{\text{OLS}}$, as Best Subset does.

This may be helpful: Why does shrinkage work?

$\endgroup$
13
  • $\begingroup$ Hmm. I don't think this answers why best subset is worse than lasso (which is the main question here). $\endgroup$
    – amoeba
    Jun 12, 2018 at 8:04
  • $\begingroup$ @amoeba: would you like to elaborate? $\endgroup$ Jun 12, 2018 at 20:39
  • $\begingroup$ Well, I understood the question as asking why lasso is preferred to best subset. Imagine we put both in a cross-validation loop, and then either tune the lasso parameter or find the best subset. The lasso is usually recommended. I understood the question as asking Why? (see e.g. the title of the Q) and I am not sure your answer actually answers that. Or did I misunderstand your answer? $\endgroup$
    – amoeba
    Jun 12, 2018 at 21:00
  • 1
    $\begingroup$ One remaining issue for Best Subset is that it constrains some parameters to be exactly zero and lets the others float freely, so there is a discontinuity in the estimate, which isn't there if we tweak the lasso $\lambda$ beyond a point $\lambda_0$ where a predictor $p$ is included or excluded. I'd argue that if we are essentially unsure whether $p$ should be included or not, because $\lambda\approx\lambda_0$, then it makes more sense to constrain the parameter estimate $\hat{\beta}_p$ via the lasso, rather than let it float freely. $\endgroup$ Jun 15, 2018 at 7:23
  • 1
    $\begingroup$ Agree that this answer doesn't really answer the question - I've added my take on this below... $\endgroup$ Jun 19, 2018 at 16:41
24
$\begingroup$

In principle, if the best subset can be found (i.e. for large enough datasets or if the solution is extremely sparse or if there are nonnegativity constraints which can help to identify the best subset), it is indeed better than the LASSO, in terms of (1) selecting the variables that actually contribute to the fit, (2) not selecting the variables that do not contribute to the fit, (3) prediction accuracy and (4) producing essentially unbiased estimates for the selected variables. One recent paper that argued for the superior quality of best subset over LASSO is that by Bertsimas et al (2016) "Best subset selection via a modern optimization lens". Another older one giving a concrete example (on the deconvolution of spike trains) where best subset was better than LASSO or ridge is that by de Rooi & Eilers (2011). (Well that one actually approximates best subset using an iterative adaptive ridge regression approach as in Frommlet & Nuel (2016) and as used in my L0glm package that is in development)

The reason that the LASSO is still preferred in practice is mostly due to it being computationally much easier to calculate and the coefficient paths being quite stable. Best subset selection, i.e. using an $L_0$ pseudonorm penalty, is often seen as a combinatorial problem, and is then NP hard, whereas the LASSO solution is convex and easy to calculate over a regularization path using pathwise coordinate descent. In fact, the LASSO ($L_1$ norm penalized regression) is the tightest convex relaxation of $L_0$ pseudonorm penalized regression / best subset selection and so is uniquely defined (bridge regression, i.e. $L_q$ norm penalized regression with q close to 0 would in principle be closer to best subset selection than LASSO, but this is no longer a convex optimization problem, and so is quite tricky to fit).

To reduce the bias of the LASSO one can use derived multistep approaches, such as the adaptive LASSO (where coefficients are differentially penalized based on a prior estimate from a least squares or ridge regression fit, e.g. just using the LASSO selected variables and then using penalty weights = 1/abs(LASSO coefficients of 1st LASSO fit)) or relaxed LASSO (a simple solution being to do a least squares fit of the variables selected by the LASSO). In comparison to best subset, LASSO tends to select slightly too many variables though. Best subset selection is better, but harder to fit.

That being said, there are also efficient computational methods now to approximate best subset selection / $L_0$ penalized regression in a very stable way. One is using the adaptive ridge approach described in the paper "An Adaptive Ridge Procedure for L0 Regularization" by Frommlet & Nuel (2016) (this approach has a nice empirical Bayesian interpretation), as implemented in the l0ara package (there referred to as broken adaptive ridge) and my own L0glm package, still in development. Here one is approximating the non-convex best subset objective by a series of convex optimization problems, which then still has a unique solution at convergence. Another one is the approach for fast best subset selection as described in the abess package, which executes in polynomial time. And there is also the L0Learn package for sparse learning & approximating the L0 penalty, also potentially in combination with L2 penalty regularisation to deal better with collinearity. All are well worth checking out. Traditional combinatorial approaches to approximate best subset, such as used in the glmulti, leaps or bestglm packages I would not recommend - they do not scale to large problem sizes & tend to be unstable.

Note that also for best subset selection you'll still have to use either cross validation or some information criterion to determine what number of predictors gives you either the best prediction performance / explanatory power for the number of variables in your model (e.g. based on the AIC or for high dimensional settings the GIC criterion) or approximates best the true support (e.g. based on BIC or for high dimensional settings eBIC or MBIC), which is essential to avoid overfitting. The paper "Extended Comparisons of Best Subset Selection, Forward Stepwise Selection, and the Lasso" by Hastie et al (2017) provides an extensive comparison of best subset, LASSO and some LASSO variants like the relaxed LASSO, and they claim that the relaxed LASSO was the one that produced the highest model prediction accuracy under the widest range of circumstances, i.e. they came to a different conclusion than Bertsimas, but that was also because they mainly looked at predictive performance and not at consistency of variable selection. So the conclusion about which is best depends a lot on what your aim is, e.g. achieving highest prediction accuracy, or achieving variable selection consistency and being best at picking out relevant variables and not including irrelevant ones or recovering the true support size most accurately. Ridge regression e.g. typically selects way too many variables but the prediction accuracy for cases with highly collinear variables can nevertheless be really good. If the true solution is really sparse there are cases where the optimal level of regularisation that results in optimal predictive performance and optimal support recovery coincide though - that is of course ideal. But of course it also depends how you would like to choose the tradeoff between Type I & Type II errors. Ideally you would like to have both a low number of false positives and a low number of false negatives, but that may only be possible with relatively large datasets & when the true solution is very sparse.

For a problem with many more cases than variables like you describe it is plain clear best subset selection is the preferred option though.

In my tests I found that an iterative adaptive ridge approach / broken adaptive ridge / L0glm approach has a series of advantages: (1) it nicely approximates best subset, (2) does so in a stable way (i.e. it has stable coefficient paths over the level of regularisation), (3) has good variable selection consistency (e.g. the union of selected variables over bootstrapped datasets typically converges quickly to a stable set, whereas LASSO would keep on adding spurious variables, which could only be mitigated by selecting the variables occurring in more than 80% of the bootstrap replicates, as in the bolasso method), (4) deals well with collinearity and high dimensional problems (e.g. with several perfectly collinear variables in the model their effect would just be spread equally - it does grouped variable selection like elastic net), (5) resolves common problems one sometimes runs into with unregularised approaches (e.g. complete separation in logistic regression), (6) has the advantage that only variables with small or zero effect are removed and that variables with large effect will stay in the model and will be unbiased, implying that omitted variable bias will be small, (7) is fast and scales well to large problems with millions of variables (computational speed is similar to glmnet), (8) allows one to choose the optimal level of regularisation a priori based on what information criterion one would like to optimise, (9) allows one to do inference via bootstrapping for both low & high dimensional problems and (10) allows one to specify prior beliefs based on the starting values given (which correspond to the square root of the SD of the Gaussian prior implied in the first iteration of the algorithm) and the prior mean (normally zero, but can also be specified to be nonzero for particular parameters if nonzero centered ridge regression updates are used). Only in high dimensional settings with an extremely sparse true solution I found MCP & SCAD penalized regression to perform better (these can be fit using ncvreg or ordinis), but even that could be covered by using MCP or SCAD penalized regression estimates as warm start starting values in an iterative adaptive ridge L0glm approach, in which case L0glm can find an even better solution still & remove some false positives.

But all this is talking about particular approximations of best subset that are uniquely determined, and are calculated either using a nonconvex penalty (in MCP and SCAD) or as a sequence of convex optimization problems (in iterative / broken adaptive ridge / L0glm). Many best subset algorithms approach the problem as a combinatorial optimization problem, and they suffer from poor stability. That, combined with the high computational cost, is probably the reason many stay away from them. But given that in my tests iterative / broken adaptive ridge / L0glm models were superior to LASSO in terms of predictive performance, biasedness of the coefficients and support recovery, stable best subset approximations should imho receive a lot more attention. See e.g. here for some examples and benchmarks that I made relating to the single or multichannel deconvolution of a blurred sparse spike train, where I compared some different approaches (L0glm, abess, L0Learn, glmnet & ordinis) (this is a problem where for single channel n=p, so in that sense a relatively small dataset, though with nonnegativity constraints on the coefficients), or a high dimensional example (with n=500 and p=1000000 variables and 50 contributing variables, so with a very sparse solution and also with nonnegativity constraints) where nonnegative MCP penalized regression fit using ordinis or L0glm with MCP as initialisation performed very well & where in both cases these approximations of best subset performed better than LASSO, both in terms of predictive performance and in terms of support recovery / variable selection consistency (ending up with just a couple of false negatives & 1 false positive), whilst also scaling well to large problems. Both of those examples are ones with nonnegativity constraints on the coefficients, which is a situation that makes recovery of the best subset a lot easier, as it acts as an effective regularisation method in itself - in fact, under that circumstance a hard thresholded version of the nonnegative LASSO outperforms plain nonnegative LASSO in support recovery and unconstrained LASSO in predictive performance, and even a hard thresholded version of nonnegative least squares would then outperform the plain LASSO.

$\endgroup$
34
  • 3
    $\begingroup$ Why is best subset the same as using L0 penalty? Best subset selects the best subset (with the lowest validation error) among subsets of any size; at least that's what OP suggested in their question. L0 penalty requires the subset to be of size $k$ (that is defined by the regularization parameter $\lambda$); one can search for best $k$ using a validation set, and then it's the best subset of size of $k$ across all possible $k$... okay, now I see that it's the same :-) $\endgroup$
    – amoeba
    Jun 19, 2018 at 18:43
  • 3
    $\begingroup$ I don't think any of the answers are addressing the problem of stability. Like stepwise and all possible subsets regression, lasso is notoriously unstable. In other words if you were to bootstrap the entire process you'll find too much arbitrariness in the list of features selected. $\endgroup$ Jun 29, 2018 at 11:16
  • 2
    $\begingroup$ Yes the variables selected by LASSO can be unstable, and this is even more so the case for best subset regression - elastic net regression is a little bit better in this respect - that tends to include far too many variables then, but selected in a more stable way, and can give better prediction accuracy under high collinearity. But a lot depends on what is the most important criterion for your application - prediction accuracy, the false positive rate of including irrelevant variables or the false negative rate of not including highly relevant variables... $\endgroup$ Jun 29, 2018 at 15:14
  • 2
    $\begingroup$ @FrankHarrell Thanks - will take a look & have a try! In the meantime here are a couple of examples for a specific deconvolution problem (single or multichannel) where the specific interest was to recover the correct support size & not just get optimal predictive performance: github.com/tomwenseleers/L0glm/blob/master/benchmarks/…. And high dimensional example here: github.com/tomwenseleers/L0glm/blob/master/benchmarks/…. Sorry code not tidied up yet... $\endgroup$ Jul 13, 2023 at 21:07
  • 3
    $\begingroup$ @FrankHarrell If you were to do a similar plot for that Gaussian high dimensional problem that I linked with n=500 and p=1000000 and 50 contributing variables but with nonnegativity constraints, github.com/tomwenseleers/L0glm/blob/master/benchmarks/…, you'd see that the accuracy & F1 score to recover the true support set with MCP or SCAD is very good, but not when using LASSO. So it's all very problem specific. Looking at one particular setting is not so informative... $\endgroup$ Jul 14, 2023 at 11:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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