Discussion: Deterioration of performance of the lasso with many predictors

1 Deterioration of the lasso

There is no question that regularization methods have become very popular in recent years. This popularity has come from three related factors: the availability of very high-dimensional data, the fact that (unlike previously available methods) regularization approaches can handle such data and the belief that the consequence of using such methods (compared to classically ‘optimal’ ones) is small.

This belief comes from so-called ‘oracle inequalities’ probability bounds on the loss

L ( λ ) = | | X ( β 0 − β ̂ λ ) | | 2 n .

L ( λ ) ≤ k σ 2 log ( p ) p 0 n

Unfortunately, these inequalities do not correspond to practical implementation of the lasso. In practice, λ is selected using a data-dependent method, such as using an information measure, generalized cross-validation or k-fold cross-validation (see, for example, Fan and Li, 2001; Wang et al., 2007; Zou et al., 2007; Zhang et al., 2010; and Flynn et al., 2013, for applications of these selectors to penalized regression estimators). Chatterjee (2014) showed that the existing theoretical results do not apply to data-dependent choice of λ, so it is unclear what oracle inequalities say about practical performance of the lasso.

We (Flynn et al., 2016) examined the question of the effect of the presence of spurious predictors on the lasso when using a data-dependent choice of λ . Specifically, in the simple case of orthonormal predictors and a single predictor with non-zero coefficient β 1 , the probability that the loss of the lasso using the optimal data-dependent choice of λ ( λ p * ) when including p predictors is larger than that when using the optimal choice based on one predictor with non-zero slope ( λ 1 * ) is arbitrarily close to 1 for an appropriately high signal to noise ratio and large p ; specifically,

Theorem: For all 1 < p ≤ n ,

Pr L p ( λ p * ) L 1 ( λ 1 * ) > 1 = Φ | β 1 | σ − 1 2 p ,

That is, the lasso deteriorates when additional superfluous predictors are added to the model being fit even when the amount of regularization is chosen optimally, something that is not true for estimation based on classical variable selection. Further, the expected amount of deterioration is infinite:

Theorem: For all 1 < p ≤ n ,

E L p ( λ p * ) L 1 ( λ 1 * ) = ∞ .

Computing the probability of deterioration in the more general case where β 0 is s -sparse (i.e., exactly s predictors have non-zero coefficients, for some 1 < = s < n ) is possible with an exact formula for λ s * . In Flynn et al., (2016), we demonstrate this calculation for s = 2 and investigate more realistic situations using simulations. The simulations show that the amount of deterioration can be far greater than oracle inequalities would imply, highlighting the danger in adding predictors with zero coefficients to the model, and the potentially misleading nature of thinking of log ( p ) as the ‘price’ of using the lasso to weed out large numbers of (unnecessary) potential predictors.

2 Categorical predictors

It is clear that this problem is a potentially serious one for data with categorical predictors, as each such predictor X j with K j categories adds K j − 1 variables to the design. For example, the data set analyzed by TG in Section 2 has only 22 categorical predictors, but this corresponds to 91 predictors from the point of view of the underlying regression model. The inclusion of interaction effects would obviously further magnify these numbers greatly. On the other hand, the group-related penalties discussed by TG effectively reduce the available degrees of freedom by imposing structure on the allowed slope estimates, potentially reducing the seriousness of this problem. In this section, we explore this question using simulations.

The underlying structure here is fairly simple. There are p potential categorical predictors, p = 2 , 3 or 4 , each having five levels (note that here p is the number of categorical predictors, not the number of underlying variables used to code any effects in fitted regressions). Only the first predictor X 1 has any predictive power for the response y , with E ( y | X 1 ) = { 12 , 9 , 6 , 3 , 0 } for the five levels. The underlying model is y i = E ( y i | X 1 i ) + ε i for the i th observation, with ε i ∼ iid N ( 0 , σ 2 ) and σ = 2 . The design is balanced, with k = 25 observations at each of the 5 p design points, making the total sample size n = 5 p + 2 . On each simulation run, for a given estimator, first the correct true model is fit (i.e., the one based only on the variable X 1 ) and the minimum average squared error of the fits (over a grid of values of the regularization parameter ranging from zero to the value that shrinks all the estimated slope coefficients to zero) is recorded, and then the model using all p predictors is fit and the minimum average squared error of the fits is determined. The ratio of these two values is the deterioration in best possible performance that comes from adding unnecessary predictors, and those values are averaged over the simulation replications. Categorical variable effects are fit based on indicator variable coding, and models for p ≥ 2 can also include interaction effects. All runs are based on 1 000 simulation replications.

Table 1 gives the results for different estimators. The first set of numbers refers to the ordinary lasso which does not take into account any of the categorical (group) structure. It is apparent that even in the case where p = 2 , there is serious deterioration in the optimal performance of the lasso, particularly when the two-way interaction is included (note that in that case, the number of potential predictors is 24, compared to the need for only four in the true model). As p increases, deterioration even with only main effects goes up, and when including interactions (corresponding to 124 predictors for the highest-order interaction model when p = 3 and 624 predictors when p = 4 , respectively), the optimal average squared error can be 15 to 30 times that when using only X 1 . These patterns (and those for the other estimators) are relatively insensitive to sample size and the value of σ , based on simulation results not presented here.

Of course, the ordinary lasso does not take into account the group structure implied by the categorical predictors. The picture is very different for the group lasso, the next entries in the table. While there is still substantial deterioration, with the optimal loss when using unnecessary predictors as much as three times that when using only X 1 , it is much smaller than that for the ordinary lasso. It is also not very sensitive to the complexity of the model being fit, even for the model including the four-way interaction when p = 4 . The ordinal lasso can only handle models with main effects, but while deterioration as p increases is still noteworthy, it is more limited, corresponding to squared error less than 50% higher than when modelling based on X 1 alone. This is consistent, of course, with the notion that putting more structure on the types of models allowed to be fit results in less of a problem from including unnecessary predictors.

The message from these results is apparent. Even in a world of ‘big data’, statistical models still matter. The less a data analyst is willing to assume, the bigger the potential price that could be paid by appealing only to automatic dimension reduction. The kinds of estimators discussed by TG can provide the important advantages that carefully reasoned appropriate assumptions about statistical structure offer while still being available for use with high-dimensional data. It is reasonable to think that this message will also apply far more widely than only when dealing with categorical data.