# Computing Cross-Validation Errors for Subset Selection: error in standard code in the literature?

I am currently trying to understand how to use cross-validation in order to choose among the "best" subsets of different sizes returned by the R function regsubsets (regsubsets returns the "best" submodel of each size, based on fit).

I found in p250 of Introduction to Statistical Learning by James, Witten, Hastie and Tibshirani the following bit of code (cited by many people online):

for(j in 1:k){
best.fit=regsubsets(Salary~.,data=Hitters[folds!=j,], nvmax=19)
for(i in 1:19){
pred=predict(best.fit, Hitters[folds==j,],id=i)
cv.errors[j,i]=mean((Hitters$Salary[folds==j]-pred)^2) } }  My question is the following: for each fold j we use regsubsets with different training data (all observations except the one in the j-fold). How do we know that we don't have a different "best" model of size i for each different fold j? Apologies if this has been asked before (could not locate it) or if the answer is trivial. ## 2 Answers Note: This post is for educational purposes. The author does not condone the usage of best subset selection for any other purpose. With that out of the way... How do we know that we don't have a different 'best' model of size i for each different fold j? You don't. If you did you'd just pick that subset and be done with it, there'd be no need for cross validation at all. The idea behind cross validation is to get an estimate of the hold out performance of a model trained on each subset size, because that's what's really important. To do so, you "fold" your data set into many different train, validate pairs, and then train and validate on each in turn. Each subset size gives you a collection of estimated validation errors, one for each fold. Each of these validation errors is a point estimate of the true hold out error, so their average should be a lower variance estimate. Now you have a single estimate of the hold out error for each subset size, so you can choose the subset size with the most desirable estimate. # For each fold for(j in 1:k){ # Fit the model with each subset of predictors on the training part of the fold best.fit=regsubsets(Salary~.,data=Hitters[folds!=j,], nvmax=19) # For each subset for(i in 1:19){ # Predict on the hold out part of the fold for that subset pred=predict(best.fit, Hitters[folds==j,],id=i) # Get the mean squared error for the model trained on the fold with the subset cv.errors[j,i]=mean((Hitters$Salary[folds==j]-pred)^2)
}
}


The missing part of this code is taking the row means of this matrix, and using the result to choose the "best" subset size.

The answer to my question is that it is possible to have a different "best" subset of fixed sized i, at each fold j. So using cross validation we only choose the optimal number of variables, not a particular subset. Having chosen the optimal size say i0, we then need to use regsubsets again on the whole training set and choose among the different "best" subsets of size i=1,...,19 the one of size i0.

To understand the reason for applying cross validation in this context in this way, see: https://www.youtube.com/watch?v=S06JpVoNaA0

• Correct, but isn't it exactly what Matthew wrote in his answer from Feb 2016? – amoeba Feb 23 '18 at 8:14
• yes you are absolutely right! – Sergios Agapiou Feb 23 '18 at 19:31