# Import necessary modules and set options
import pandas as pd
import numpy as np
import itertools
import matplotlib.pyplot as plt
import statsmodels.api as sm
from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet, RidgeCV, LassoCV, ElasticNetCV, LarsCV, lars_path
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import mean_squared_error
import warnings
warnings.filterwarnings("ignore")
prostate = pd.read_csv('../data/Prostate_Cancer.txt', delimiter=',')
prostate.head(10)
prostate = prostate.drop(['id'], axis=1)
prostate.shape
# Train-test split
y_train = np.array(prostate[prostate.train == "T"]['lpsa'])
y_test = np.array(prostate[prostate.train == "F"]['lpsa'])
X_train = np.array(prostate[prostate.train == "T"].drop(['lpsa', 'train'], axis=1))
X_test = np.array(prostate[prostate.train == "F"].drop(['lpsa', 'train'], axis=1))
linreg_model = LinearRegression(normalize=True).fit(X_train, y_train)
linreg_prediction = linreg_model.predict(X_test)
linreg_mse = np.mean((y_test - linreg_prediction)**2)
linreg_coefs = dict(
zip(['Intercept'] + prostate.columns.tolist()[:-1],
np.round(np.concatenate((linreg_model.intercept_, linreg_model.coef_),
axis=None), 3))
)
print('Linear Regression MSE: {}'.format(np.round(linreg_mse, 3)))
print('Linear Regression coefficients:')
linreg_coefs
On est habitué à ça
X1_train = sm.add_constant(X_train) # adding a constant
mod = sm.OLS(y_train, X1_train)
res = mod.fit()
print(res.summary())
Une approche simple pour choisir un sous-ensemble de variables pour la régression linéaire consiste à essayer toutes les combinaisons possibles et à en choisir une qui minimise un critère certain. C'est ce que vise la recherche du meilleur sous-ensemble de variables pour la régression. Pour chaque $k \in \{1, 2, \ldots, p\}$, où $p$ est le nombre total de variables disponibles, il sélectionne le sous-ensemble de taille $k$ qui donne la plus petite somme résiduelle de carrés. Cependant, la somme des carrés ne peut pas être utilisée comme critère pour déterminer $k$ elle-même, car elle diminue nécessairement avec $k$: plus le nombre de variables incluses dans le modèle est grand, plus ses résidus sont petits. Cela ne garantit toutefois pas de meilleures performances prédictives. C’est pourquoi un autre critère devrait être utilisé pour sélectionner le modèle final. Pour les modèles axés sur la prédiction, une erreur (éventuellement croisée) sur les données de test est un choix courant. Comme la recherche du meilleur sous-ensemble de régression n'est implémentée dans aucun package Python, nous devons effectuer une boucle sur $k$ et tous les sous-ensembles de taille $k$ manuellement. Le bout de code suivant implémente la méthode :
results = pd.DataFrame(columns=['num_features', 'features', 'MSE'])
# Loop over all possible numbers of features to be included
for k in range(1, X_train.shape[1] + 1):
# Loop over all possible subsets of size k
for subset in itertools.combinations(range(X_train.shape[1]), k):
subset = list(subset)
linreg_model = LinearRegression(normalize=True).fit(X_train[:, subset], y_train)
linreg_prediction = linreg_model.predict(X_test[:, subset])
linreg_mse = np.mean((y_test - linreg_prediction)**2)
results = results.append(pd.DataFrame([{'num_features': k,
'features': subset,
'MSE': linreg_mse}]))
# Inspect best combinations
results = results.sort_values('MSE').reset_index()
print(results.head())
# Fit best model
best_subset_model = LinearRegression(normalize=True).fit(X_train[:, results['features'][0]], y_train)
best_subset_coefs = dict(
zip(['Intercept'] + prostate.columns.tolist()[:-1],
np.round(np.concatenate((best_subset_model.intercept_, best_subset_model.coef_), axis=None), 3))
)
print('Best Subset Regression MSE: {}'.format(np.round(results['MSE'][0], 4)))
print('Best Subset Regression coefficients:')
best_subset_coefs
n_alphas = 300
alphas = np.logspace(-5, 4, n_alphas)
coefs = []
for a in alphas:
ridge = Ridge(alpha=a, fit_intercept=False)
ridge.fit(X_train, y_train)
coefs.append(ridge.coef_)
# #############################################################################
# Display results
ax = plt.gca()
ax.plot(alphas, coefs)
ax.set_xscale('log')
ax.set_xlim(ax.get_xlim()[::-1]) # reverse axis
plt.xlabel('alpha')
plt.ylabel('weights')
plt.title('Ridge coefficients as a function of the regularization')
plt.axis('tight')
plt.show()
ridge_cv = RidgeCV(normalize=True, alphas=np.logspace(-5, 4, n_alphas), store_cv_values=True)
ridge_model = ridge_cv.fit(X_train, y_train)
ridge_prediction = ridge_model.predict(X_test)
ridge_mse = np.mean((y_test - ridge_prediction)**2)
ridge_coefs = dict(
zip(['Intercept'] + prostate.columns.tolist()[:-1],
np.round(np.concatenate((ridge_model.intercept_, ridge_model.coef_),
axis=None), 3))
)
print('Ridge Regression MSE: {}'.format(np.round(ridge_mse, 3)))
print('Ridge Regression coefficients:')
ridge_coefs
errors = ridge_cv.cv_values_.mean(axis=0)
ax = plt.gca()
ax.plot(alphas, errors)
ax.set_xscale('log')
plt.xlabel('alpha')
plt.ylabel('error')
plt.axis('tight')
plt.show()
print("Computing regularization path using the LARS ...")
_, _, coefs = lars_path(X_train, y_train, method='lasso', verbose=True)
xx = np.sum(np.abs(coefs.T), axis=1)
xx /= xx[-1]
plt.plot(xx, coefs.T)
ymin, ymax = plt.ylim()
plt.vlines(xx, ymin, ymax, linestyle='dashed')
plt.xlabel('|coef| / max|coef|')
plt.ylabel('Coefficients')
plt.title('LASSO Path')
plt.axis('tight')
plt.show()
lasso_cv = LassoCV(normalize=True, alphas=np.logspace(-5, 4, n_alphas))
lasso_model = lasso_cv.fit(X_train, y_train)
lasso_prediction = lasso_model.predict(X_test)
lasso_mse = np.mean((y_test - lasso_prediction)**2)
lasso_coefs = dict(
zip(['Intercept'] + prostate.columns.tolist()[:-1],
np.round(np.concatenate((lasso_model.intercept_, lasso_model.coef_), axis=None), 3))
)
print('LASSO MSE: {}'.format(np.round(lasso_mse, 3)))
print('LASSO coefficients:')
lasso_coefs
Erreur quadratique moyenne calculée par validation croisée en fonction de la valeur de alpha
errors = lasso_cv.mse_path_.mean(axis=1)
ax = plt.gca()
ax.plot(alphas, errors)
ax.set_xscale('log')
plt.xlabel('alpha')
plt.ylabel('error')
plt.axis('tight')
plt.show()
elastic_net_cv = ElasticNetCV(normalize=False, alphas=np.logspace(-5, 4, 400),
l1_ratio=np.linspace(0, 1, 100), cv=3)
elastic_net_model = elastic_net_cv.fit(X_train, y_train)
elastic_net_prediction = elastic_net_model.predict(X_test)
elastic_net_mse = np.mean((y_test - elastic_net_prediction)**2)
elastic_net_coefs = dict(
zip(['Intercept'] + prostate.columns.tolist()[:-1],
np.round(np.concatenate((elastic_net_model.intercept_,
elastic_net_model.coef_), axis=None), 3))
)
print('Best alpha : {}'.format(np.round(elastic_net_model.alpha_, 3)))
print('Best l1_ratio : {}'.format(np.round(elastic_net_model.l1_ratio_, 3)))
print('Elastic Net MSE: {}'.format(np.round(elastic_net_mse, 3)))
print('Elastic Net coefficients:')
elastic_net_coefs
GridSearchCV¶# Create the hyperparameter grid
alphas = np.logspace(-5, 4, n_alphas)
l1_space = np.linspace(0, 1, 100)
param_grid = {'l1_ratio': l1_space, 'alpha':alphas}
# Instantiate the ElasticNet regressor: elastic_net
elastic_net = ElasticNet(normalize=True)
# Setup the GridSearchCV object: gm_cv
gm_cv = GridSearchCV(elastic_net, param_grid, cv=3, n_jobs=3) ## attention n_jobs=-1 va utiliser tous les CPUs
# Fit it to the training data
np.random.seed(11)
gm_cv.fit(X_train, y_train)
# Predict on the test set and compute metrics
y_pred = gm_cv.predict(X_test)
r2 = gm_cv.score(X_test, y_test)
mse = mean_squared_error(y_test, y_pred)
print("Tuned ElasticNet alpha and l1 ratio: {}".format(gm_cv.best_params_))
print("Tuned ElasticNet R squared: {}".format(r2))
print("Tuned ElasticNet MSE: {}".format(mse))
best_coefs = dict(zip(['Intercept'] + prostate.columns.tolist()[:-1],
np.round(np.concatenate((gm_cv.best_estimator_.intercept_,
gm_cv.best_estimator_.coef_), axis=None), 3)))
print("Best model:")
best_coefs