Nested cross-validation and selecting the best regression model - is this the right SKLearn process?

If I understand correctly, nested-CV can help me evaluate what model and hyperparameter tuning process is best. The inner loop (GridSearchCV) finds the best hyperparameters, and the outter loop (cross_val_score) evaluates the hyperparameter tuning algorithm. I then choose which tuning/model combo from the outer loop that minimizes mse (I'm looking at regression classifier) for my final model test.

I've read the questions/answers on nested-cross-validation, but haven't seen an example of a full pipeline that utilizes this. So, does my code below (please ignore the actual hyperparameter ranges - this is just for example) and thought process make sense?

from sklearn.cross_validation import cross_val_score, train_test_split
from sklearn.grid_search import GridSearchCV
from sklearn.metrics import mean_squared_error
from sklearn.ensemble import RandomForestRegressor
from sklearn.svm import SVR
from sklearn.datasets import make_regression

# create some regression data
X, y = make_regression(n_samples=1000, n_features=10)
params = [{'C':[0.01,0.05,0.1,1]},{'n_estimators':[10,100,1000]}]

# setup models, variables
mean_score = []
models = [SVR(), RandomForestRegressor()]

# split into train and test sets
X_train, X_test, y_train, y_test = train_test_split(X, y, train_size=0.3)

# estimate performance of hyperparameter tuning and model algorithm pipeline
for idx, model in enumerate(models):
clf = GridSearchCV(model, params[idx], scoring='mean_squared_error')

# this performs a nested CV in SKLearn
score = cross_val_score(clf, X_train, y_train, scoring='mean_squared_error')

# get the mean MSE across each fold
mean_score.append(np.mean(score))
print('Model:', model, 'MSE:', mean_score[-1])

# estimate generalization performance of the best model selection technique
best_idx = mean_score.index(max(mean_score)) # because SKLearn flips MSE signs, max works OK here
best_model = models[best_idx]

clf_final = GridSearchCV(best_model, params[best_idx])
clf_final.fit(X_train, y_train)

y_pred = clf_final.predict(X_test)
rmse = np.sqrt(mean_squared_error(y_test, y_pred))

print('Final Model': best_model, 'Final model RMSE:', rmse)


Yours is not an example of nested cross-validation.

Nested cross-validation is useful to figure out whether, say, a random forest or a SVM is better suited for your problem. Nested CV only outputs a score, it does not output a model like in your code.

This would be an example of nested cross validation:

from sklearn.datasets import load_boston
from sklearn.cross_validation import KFold
from sklearn.metrics import mean_squared_error
from sklearn.grid_search import GridSearchCV
from sklearn.ensemble import RandomForestRegressor
from sklearn.svm import SVR
import numpy as np

params = [{'C': [0.01, 0.05, 0.1, 1]}, {'n_estimators': [10, 100, 1000]}]
models = [SVR(), RandomForestRegressor()]

X = df['data']
y = df['target']

cv = [[] for _ in range(len(models))]
for tr, ts in KFold(len(X)):
for i, (model, param) in enumerate(zip(models, params)):
best_m = GridSearchCV(model, param)
best_m.fit(X[tr], y[tr])
s = mean_squared_error(y[ts], best_m.predict(X[ts]))
cv[i].append(s)
print(np.mean(cv, 1))


By the way, a couple of thoughts:

• I see no purpose to grid search for n_estimators for your random forest. Obviously, the more, the merrier. Things like max_depth is the kind of regularization that you want to optimize. The error for the nested CV of RandomForest was much higher because you did not optimize for the right hyperparameters, not necessarily because it is a worse model.
• You might also want to try gradient boosting trees.
• Thanks for that. My goal is to do exactly what you said - figure out what classifier algorithm would be best suited for my problem. I guess I'm confused in terms of the documentation of SKLearn: scikit-learn.org/stable/tutorial/statistical_inference/… (under 'nested cross-validation') Aug 5, 2016 at 11:26
• To test the performance of the best-selected model, would I do a final cross-validation on the whole dataset? Or should I split my dataset into train/test BEFORE nested CV, run nested CV on the train, and then fit the best model on the train data and test on test? Aug 5, 2016 at 11:38
• Sorry for the comment barrage. So my final model would be: best_idx = np.where(np.mean(cv,1).min())[0]; final_m = GridSearchCV(models[best_idx], params[best_idx]); final_m.fit(X,y) Aug 5, 2016 at 11:59
• Building off what you said, this was what I was going for with built-in SKLearn functions (gives the same as your answer): for model, param in zip(models, params): clf = GridSearchCV(model, param) my_score = cross_val_score(clf, X, y, scoring='mean_squared_error') my_scores.append(my_score) Aug 5, 2016 at 12:27

Nested cross validation estimates the generalization error of a model, so it is a good way to choose the best model from a list of candidate models and their associated parameter grids. The original post is close to doing nested CV: rather than doing a single train–test split, one should instead use a second cross-validation splitter. That is, one "nests" an "inner" cross-validation splitter inside an "outer" cross validation splitter.

The inner cross-validation splitter is used to choose hyperparameters. The outer cross-validation splitter averages the test error over multiple train–test splits. Averaging the generalization error over multiple train–test splits provides a more reliable estimate of the accuracy of the model on unseen data.

I modified the original post's code to update it to the latest version of sklearn (with sklearn.cross_validation superseded by sklearn.model_selection and with 'mean_squared_error' replaced by 'neg_mean_squared_error'), and I used two KFold cross-validation splitters to select the best model. To learn more about nested cross validation, see the sklearn's example on nested cross-validation.

from sklearn.model_selection import KFold, cross_val_score, GridSearchCV
from sklearn.datasets import make_regression
from sklearn.ensemble import RandomForestRegressor
from sklearn.svm import SVR
import numpy as np

# outer_cv creates 3 folds for estimating generalization error
outer_cv = KFold(3)

# when we train on a certain fold, we use a second cross-validation
# split in order to choose hyperparameters
inner_cv = KFold(3)

# create some regression data
X, y = make_regression(n_samples=1000, n_features=10)

# give shorthand names to models and use those as dictionary keys mapping
# to models and parameter grids for that model
models_and_parameters = {
'svr': (SVR(),
{'C': [0.01, 0.05, 0.1, 1]}),
'rf': (RandomForestRegressor(),
{'max_depth': [5, 10, 50, 100, 200, 500]})}

# we will collect the average of the scores on the 3 outer folds in this dictionary
# with keys given by the names of the models in models_and_parameters
average_scores_across_outer_folds_for_each_model = dict()

# find the model with the best generalization error
for name, (model, params) in models_and_parameters.items():
# this object is a regressor that also happens to choose
# its hyperparameters automatically using inner_cv
regressor_that_optimizes_its_hyperparams = GridSearchCV(
estimator=model, param_grid=params,
cv=inner_cv, scoring='neg_mean_squared_error')

# estimate generalization error on the 3-fold splits of the data
scores_across_outer_folds = cross_val_score(
regressor_that_optimizes_its_hyperparams,
X, y, cv=outer_cv, scoring='neg_mean_squared_error')

# get the mean MSE across each of outer_cv's 3 folds
average_scores_across_outer_folds_for_each_model[name] = np.mean(scores_across_outer_folds)
error_summary = 'Model: {name}\nMSE in the 3 outer folds: {scores}.\nAverage error: {avg}'
print(error_summary.format(
name=name, scores=scores_across_outer_folds,
avg=np.mean(scores_across_outer_folds)))
print()

print('Average score across the outer folds: ',
average_scores_across_outer_folds_for_each_model)

many_stars = '\n' + '*' * 100 + '\n'
print(many_stars + 'Now we choose the best model and refit on the whole dataset' + many_stars)

best_model_name, best_model_avg_score = max(
average_scores_across_outer_folds_for_each_model.items(),
key=(lambda name_averagescore: name_averagescore[1]))

best_model, best_model_params = models_and_parameters[best_model_name]

# now we refit this best model on the whole dataset so that we can start
# making predictions on other data, and now we have a reliable estimate of
# this model's generalization error and we are confident this is the best model
# among the ones we have tried
final_regressor = GridSearchCV(best_model, best_model_params, cv=inner_cv)
final_regressor.fit(X, y)

print('Best model: \n\t{}'.format(best_model), end='\n\n')
print('Estimation of its generalization error (negative mean squared error):\n\t{}'.format(
best_model_avg_score), end='\n\n')
print('Best parameter choice for this model: \n\t{params}'
'\n(according to cross-validation {cv} on the whole dataset).'.format(
params=final_regressor.best_params_, cv=inner_cv))

• At the very last comment, you say you "...refit this best model on the whole training set" but you actually do it on the whole data set (X and y). As far as I understand this is the right thing to do, but then the comment has to be corrected. What do you think? Jul 24, 2017 at 7:40
• Thanks @DrorAtariah for catching that. You're right. I fixed it. Jul 24, 2017 at 13:46

You do not need

# this performs a nested CV in SKLearn
score = cross_val_score(clf, X_train, y_train, scoring='mean_squared_error')


GridSearchCV does this for you. To get intuition of the grid-search process, try to use GridSearchCV(... , verbose=3)

To extract scores for each fold see this example in scikit-learn documentation

• I thought grid search was solely for optimizing hyper parameters? How would I use gridsearch in conjunction with something else to figure out the best classifier algorithm (i.e. SVR vs. RandomForest)? Aug 5, 2016 at 11:31
• Yes. For each combinations of hyper-parameters GridSearchCV makes folds and calculates the scores (mean squared error in your case) on the left-out data. Thus each combination of hyper-parameters gets its own mean score. The "optimization" is just choosing the combination with the best mean score. You can extract those mean scores and compare them directly for various models. Aug 5, 2016 at 11:47