Safer fitting through regularization

Esben Jannik Bjerrum/ January 6, 2016/ Blog, Cheminformatics, Machine Learning, RDkit/ 1 comments

Last time a simple multiple linear regression (MLR) model was seriously overfitted to molecular solubility data. This time the concept of regularization will be tested.
Recall that the MLR model was of the form:
y = b1*x1 + b2*x2 + … + bi*xi + b0
Where each x correspond to the absence or presence of a particular bit in the morgan fingerprint and the coefficients b1 to bi was fit to the solubility data in y. With 1024 features and only 1311 molecular samples the model was seriously over-fit and had no predictive power in cross validation.
Regularization is the concept of penalizing “complex” solutions. So what is a complex solution and what is a simple solution? A complex solution may be one that has a lot of large coefficients with opposing signs, whereas a simple solution could be seen as one where there was a lot of small contributions with low coefficients. This can be build into the fit itself.
The coefficients in the MLR is fit so that the sum of squared errors between the predicted y and the measured y is lowest possible. The average squared error is usually used as the cost function that is optimized in the machine learning model. By adding an extra term to the cost function that penalizes “complex” solutions it is possible to regularize the solution so that over-fitting is avoided. For MLR the sum of squared coefficients are usually added to the cost-function multiplied by a parameter alpha that thus controls the amount of regularization.
This regularized MLR is available in Scikit-learn as Ridge Regression. It will be explored what happens in cross-validation of the previously generated data when varying the regularization parameter alpha.
The X matrix and y vector were produced in the same way as in the last blog post. In python a tuning of the regularization could be implemented in this way.

#prepare scoring lists
fitscores = []
predictscores = []
#prepare a log spaced list of alpha values to test
alphas = np.logspace(-1, 4, num=25)
#Iterate through alphas and fit with Ridge Regression
for alpha in alphas:
	estimator = linear_model.Ridge(alpha = alpha)
	estimator.fit(X_train,y_train)
	fitscores.append(estimator.score(X_train,y_train))
	predictscores.append(estimator.score(X_test,y_test))

#show a plot
import matplotlib.pyplot as plt
ax = plt.gca()
ax.set_xscale('log')
ax.plot(alphas, fitscores,'g')
ax.plot(alphas, predictscores,'b')
plt.xlabel('alpha')
plt.ylabel('Correlation Coefficient')
plt.show()

estimator= linear_model.Ridge(alpha = 10)
estimator.fit(X_train,y_train)
#Predict the test set
y_pred = estimator.predict(X_test)
rms = (np.mean((y_test - y_pred)**2))**0.5
s = np.std(y_test -y_pred)
print rms, s