Cross-validation for ridge regression is selecting too low value of the regularization parameter

Question

I perform ridge regression for classification. To find regularization parameter I do K-fold cross-validation with classification accuracy as a measure. This gives me some $\lambda$, which I then use in training of a final model on the whole available training data. The problem is that when I take 10*$\lambda$ my test accuracy on separate dataset is much better than with $\lambda$. I cannot see a reason for that. Tell me please, why this might happen?

The lambda I get is 10^4 and beta coefficients are about 10^(-3). I have about 15000 features, which I standardize before doing regression.

Answer 1

1. Cross validation estimates are known to have high variance, so there is no guarantee you always get exactly the best result

2. As Theja has hinted, search for the best lambda adapts to the training sample and so can introduce "second level overfit". One recommendation is "One Standard Error Rule": the ideas is to increase lambda taking variability into account. Please see books: CART, ESL

3. A more reliable way to select lambda is through marginal likelihood or "evidence", see this book. Can get computationally intense, I don't see it used very often.

Answer 2

I tend not to use classification accuracy as a model selection criterion because it is discrete and a small change in lambda can produce a larger change in the error rate. I find that Allen's PRESS (i.e. the cross-validation estimate of the sum of squares error) works well for setting lambda in classification problems. It also has a certain symmetry to it, if we use least squares for fitting the ridge regression model, then using it for tuning lambda as well ought to be reasonable. As it is a classification task, a regularised logistic regression model might be a better option, with lambda chosen to minimise the cross-validated negative log-likelihood.

Answer 3

Based on the values you provide, here is what I think is happening. Your $\lambda$ value is incredibly high and your $\beta$s are very small. In essence, you are overfitting the data and modeling the noise. Not sure how many observations you have, but 15000 features is a lot and hopefully your ratio of $p/n$ is not astronomical. I am assuming that in your data, $p>n$ so here is one recommendation or steps to try.

Split the data into a training and testing set (80-20$\%$ ratio).

Standardize the training set. Use the same parameters to standardize the testing set. So the mean of the columns of $X_{test} $ won't be zero or the std. won't be exactly 1, but that is okay. You are looking to generalize your results. It is important to note that there is no column of 1's in your training or testing set (not evaluating the intercept). Since you are doing classification, I assume that this is logisitic ridge regression with class labels within $(0,1)$.

Within the CV steps for the $k^{th}$ fold, for different values of $\lambda$, evaluate the generalized cross validation error. $GCVE_{\lambda,k} =$ $\frac{\frac{1}{j}\sum_{i=1}^{j}(y_i-\hat{y_i})^2}{ (n_{cv} - df)} $. Here, $n_{cv}$ is the total number of samples used across the other $(k-1)$ CV folds in training the model, $df$ is the degree of freedom of the model which is the $trace$ of $X(X^\top X + \lambda I)^{-1} X^\top$ (where $X$ is constructed from the data in the $(k-1)$ folds), $j$ is the number of samples in the $k^{th}$ fold, and ($\hat{y_i}$) are the predicted class labels. Note that if you add a 1 to the denominator of $GCVE_{\lambda,k}$ and set $\lambda=0$, the entire model reduces to straightforward logistic regression with no regularization. Without a loss of generalizability, we can assume that $\lambda>0$ in your CV-folds and compute the average $GCVE_\lambda$ across the $k$ folds. Also compute the mean $\hat{\beta_\lambda}$'s from the CV folds for each $\lambda$. It is very important that you have an equal proportion of both classes in the training set, otherwise the solution can sometimes degenerate and simply just predict one class label all the time.

A plot of $\lambda$ (x-axis) vs. mean $GCVE_\lambda$ (y-axis) should give you a montonically decreasing trace. As you increase $\lambda$, you are reducing the $df$ in the model or its ability to do anything useful. However, at a certain $\lambda$, the $GCVE_\lambda$ will taper and not change much; this decay happens usually for small values of $\lambda$. At the same time, the ridge trace ($\lambda$ vs mean $\hat{\beta_\lambda}$) should hopefully tell you that the above value of $\lambda$ coincides with stability in the ridge trace computed from the CV procedure (no wild oscillations). Basically you will have to eye ball the value of $\lambda$ that is in 'steady-state' mode in both plots and smaller the $\lambda$, the better. It is important to note what $\lambda$ actually does; you are basically correcting for multicollinearity or non-independence of the predictors. Hopefully the above traces converge to steady state quickly; larger the values of $\lambda$ and more imperfect is your assumption that there is truly a relationship between the predictors and the class variable.

For this visually determined value of $\lambda$, train the entire set with all training data to obtain the $\hat{\beta}$ parameters.

Now evaluate the model by trying to predict the $y_{test}$ values; you can compute $R^2$ from this result or you can test the null hypothesis of obtaining the test classification accuracy due to chance using the binomial distribution (for 2 class case). If your testing set is imbalanced, use balanced accuracy as a metric rather than simple accuracy. Here, importance is also placed for misclassification and is given by $BA = 0.5(TP/(TP+FN) + TN/(TN+FP))$. TP-True Positive. TN-True Negative. FN-False negative. FP-fasle positive.

You can repeat the above procedure for different partitions of the training and testing set to get a distribution of $R^2$ or $p$ values to run some simple statistics at the global level.

Personally, I would recommend running a permutation test to evaluate the significance of the $\beta$'s. Here, across many iterations (say 10000 times), the elements of the $y$ are shuffled so that the elements in row $j$ of $X$ do not correspond to the $j^{th}$ element of $y$. This procedure obtains a distribution of each $\beta_i, i=i...p$. Evaluate p-value as the proportion of the magnitude of those values obtained from the permutation > $|\beta_i| $. FDR (False discovery rate) can be used for correction as you are now running 15000 tests and you want to protect against a Type I error. It is important to note that if you are running this latter procedure for the significance test of the betas, partitioning into training and test set is unnecessary and you can use the entire dataset as the 'training set', with the cross-validation procedure. For example looking at your obtained betas (all of which are $10^{-3}$), none of them would be significant from a permutation test as the magnitude of those obtained from the permutations would likely be of the same value. Try tuning your model using the above procedure and if you still don't get any results, then it looks like there is really no relationship between your predictors and the class labels. Hope this helps...

Answer 4

Consider the following steps:

1. Pick $\lambda$ using training set where it was the best in some sense (best average performance on validation folds, which are subsets of the training set).

2. Fix this $\lambda$ and get a new regression model using the training set.

3. Evaluate the model on a held out test set (which we hope was drawn from the same distribution as the training set).

Think of what we did in steps 1 and 2. All we did was choose the coefficients of the ridge model (say $\beta$) and the regularization coefficient $\lambda$ using the training set only. In other words, Steps 1 and 2 can be seen as a black box which optimized all the parameters of our full model, a.k.a $(\beta,\lambda)$, using the training data.

In the above point of view, there is no reason why this full model $(\beta,\lambda)$ has to perform better than a new full model $(\beta,10*\lambda)$ we come up with after getting feedback from the held out test set.

Although we expect that the best performing model on the training set translates to great performance on the held out test set, it cannot beat a model which we come up with after looking at/taking into account the held out test set.