TP++ - Régression linéaire

(les TP++ sont des compléments à l’enseignement de base et ne font pas partie des séances de travaux pratiques avec encadrement)

_images/jupyter_logo.png

Cahier Jupyter

Références externes utiles :

L’objectif de cette séance de TP est de présenter l’utilisation des fonctionnalités de Scikit-learn concernant la régression linéaire. Des transparents de présentation de la régression linéaire sont accessibles ici.

Régression linéaire - travaux pratiques

La régression linéaire avec Scikit-learn

Dans Scikit-learn, la régression linéaire est mise en œuvre dans la classe sklearn.linear_model.LinearRegression (voir https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html).

Pour s’en servir, on appelle LinearRegression(*, fit_intercept=True, copy_X=True). Les principaux paramètres d’appel sont (pour des explications plus détaillées voir le lien ci-dessus) :

  • fit_intercept : si False, l’intercept (paramètre \(w_0\)) n’est pas calculé ; par défaut True.

  • copy_X : si False, la matrice de données X est écrasée par les données transformées ; par défaut True.

Les principaux attributs accessibles sont les suivants :

  • coef_ : array de forme (n_features, ) ou (n_targets, n_features) correspondant les coefficients estimés pour le problème de régression linéaire ; n_features est le nombre de variables explicatives, n_targets le nombre de variables expliquées.

  • intercept_ : float si une seule variable expliquée ou array de forme (n_targets,) si n_targets variables expliquées ; valeur(s) d’intercept (paramètre(s) \(w_0\)).

  • n_features_in_ : nombre de variables explicatives trouvées (features).

  • rank_ : le rang de la matrice X dont les colonnes sont les variables explicatives et les lignes les observations ; disponible seulement pour X dense.

  • singular_ : valeurs propres de X ; disponible seulement pour X dense.

Les principales méthodes qui peuvent être employées :

  • fit(X, y, sample_weight=None) : estimer les paramètres du modèle linéaire à partir des variables explicatives de X et de la (ou des) variable(s) expliquée(s) de y ; sample_weight permet d’employer des pondérations différentes pour les différentes observations.

  • get_params() : retourner les paramètres du modèle linéaire estimé.

  • predict(X) : prédire les valeurs de la (des) variables expliquée(s) à partir des valeurs des variables explicatives fournies dans X.

  • score(X, y, sample_weight=None) : calculer le coefficient de détermination \(R^2\) pour les observations définies par les variables explicatives de X et la (ou des) variable(s) expliquée(s) de y ; sample_weight permet d’employer des pondérations différentes pour les observations.

Régression linéaire sur des données synthétiques

Importation de bibliothèques :

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

Une seule variable explicative

Nous générons d’abord des données bidimensionnelles autour d’une droite de pente pente qui passe par l’origine. La dispersion des données autour de cette droite est réglée à l’aide de ry :

rndn2d = np.random.randn(100,2)  # données 2D suivant une loi normale de moyenne 0 et variance 1
pente = 0.8 # paramètre qui règle la pente de la droite
ry = 3      # paramètre qui règle la dispersion autour de la droite
donnees = np.column_stack((10*rndn2d[:,0], 10*pente*rndn2d[:,0]+ry*rndn2d[:,1]))
print(donnees[:3,:])

%matplotlib widget

# Affichage des données générées
fig = plt.figure()
plt.scatter(donnees[:,0], donnees[:,1])
plt.xlabel("x")
plt.ylabel("y")
plt.title("Données avec une variable explicative")
plt.show()

Nous pouvons maintenant estimer les paramètres de la droite de régression :

X,y = np.hsplit(donnees, 2)
regresseur = LinearRegression().fit(X, y)
print("Pente :", pente)
print("Pente estimée :", regresseur.coef_)
print("Intercept :", regresseur.intercept_)
print("R2 :", regresseur.score(X, y))

Nous pouvons afficher la droite de régression avec les données :

xmin = X.min()
xmax = X.max()
ymin = (xmin * regresseur.coef_[0][0]) + regresseur.intercept_[0]
ymax = (xmax * regresseur.coef_[0][0]) + regresseur.intercept_[0]
projectionsVerticales = (donnees[:,0] * regresseur.coef_[0][0]) + regresseur.intercept_[0]
# Affichage des données générées et de la droite de régression
fig = plt.figure()
plt.scatter(donnees[:,0], donnees[:,1])
plt.plot([xmin, xmax], [ymin, ymax], 'r-')
plt.vlines(donnees[:,0], projectionsVerticales, donnees[:,1], linestyles='dashed')
plt.xlabel("x")
plt.ylabel("y")
plt.title("Données et droite de régression")
plt.show()

Question :

Sans modifier la pente (pente), réduire d’abord et ensuite augmenter la dispersion autour de la droite (ry), générez de nouveau des données et apliquez sur ces données la régression linéaire. Comment évolue l’écart entre la pente estimée et la vraie pente ? Et l’intercept ? Et le coefficient de détermination (R2) ?

Correction :

rndn2d = np.random.randn(100,2)
ry = 1      # paramètre qui règle la dispersion autour de la droite
donnees1 = np.column_stack((10*rndn2d[:,0], 10*pente*rndn2d[:,0]+ry*rndn2d[:,1]))
# Affichage des données générées
fig = plt.figure()
plt.scatter(donnees1[:,0], donnees1[:,1])
plt.xlabel("x")
plt.ylabel("y")
plt.title("ry = 1")
plt.show()
X,y = np.hsplit(donnees1, 2)
regresseur1 = LinearRegression().fit(X, y)
print("Pente :", pente)
print("Pente estimée :", regresseur1.coef_)
print("Intercept :", regresseur1.intercept_)
print("R2 :", regresseur1.score(X, y))

Deux variables explicatives

Nous générons maintenant des données 3D, avec deux variables explicatives et une variable expliquée. La variable expliquée dépend seulement de la première des variables explicatives :

rndn2d = np.random.randn(100,2)  # données 2D suivant une loi normale de moyenne 0 et variance 1
rndn1d = np.random.randn(100,1)  # données 1D suivant une loi normale de moyenne 0 et variance 1
pente2 = 0.8 # paramètre qui règle la pente de la droite par rapport à la 1ère var. explicative
rzx = 6      # paramètre qui règle la dispersion autour de la dépendance à la var.1
rzy = 3      # paramètre qui règle la variance de la var.2
# génération des données ; la seconde variable 'explicative' ne contribue pas à la variable expliquée
donnees3d = np.column_stack((10*rndn2d[:,0], rzy*rndn1d, 10*pente*rndn2d[:,0]+rzx*rndn2d[:,1]))

# Activer le mode interactif de Matplotlib pour les notebooks (permet de faire tourner les graphiques 3D)
%matplotlib widget

# Affichage des données générées
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(donnees3d[:,0], donnees3d[:,1], donnees3d[:,2])
plt.xlabel("ve1")
plt.ylabel("ve2")
plt.title("Données avec 2 variables explicatives")
plt.show()

Nous pouvons maintenant estimer les paramètres du plan de régression :

print("donnees3d :\n", donnees3d[:2,:])  # affichage 2 premières lignes
X,y = np.array_split(donnees3d, 2, axis=1)
print("X :\n", X[:2])  # affichage 2 premières lignes
print("y :\n", y[:2])  # affichage 2 premières lignes
regresseur3d = LinearRegression().fit(X, y)
print("Pente par rapport à la 1ère var. explicative :", pente)
print("Paramètres estimés par rapport aux 2 var. explicatives :", regresseur3d.coef_)
print("Intercept :", regresseur3d.intercept_)
print("R2 :", regresseur3d.score(X, y))

Question :

Modifier rzx, puis rzy, générer de nouveau les données et regarder l’impact de ces modifications sur les paramètres estimés et sur le coefficient de détermination R2.

Correction :

rzx = 6      # paramètre qui règle la dispersion autour de la dépendance à la var.1
rzy = 30      # paramètre qui règle la variance de la var.2
# génération des données ; la seconde variable 'explicative' ne contribue pas à la variable expliquée
donnees3d = np.column_stack((10*rndn2d[:,0], rzy*rndn1d, 10*pente*rndn2d[:,0]+rzx*rndn2d[:,1]))

X,y = np.array_split(donnees3d, 2, axis=1)
print("X :\n", X[:2])  # affichage 2 premières lignes
print("y :\n", y[:2])  # affichage 2 premières lignes
regresseur3d = LinearRegression().fit(X, y)
print("Pente par rapport à la 1ère var. explicative :", pente)
print("Paramètres estimés par rapport aux 2 var. explicatives :", regresseur3d.coef_)
print("Intercept :", regresseur3d.intercept_)
print("R2 :", regresseur3d.score(X, y))

La régularisation \(L_1\) (LASSO) permet d’obtenir un modèle linéaire dont sont exclues les variables explicatives qui n’ont pas un impact significatif sur la variable expliquée : les paramètres associés à ces variables (pas vraiment) explicatives sont 0.0 à l’issue de l’estimation du modèle. Nous pouvons appliquer la régression linéaire avec régularisation \(L_1\) grâce à la classe sklearn.linear_model.Lasso (voir https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Lasso.html). L’appel contient un paramètre supplémentaire alpha qui est la pondération du terme de régularisation ; la valeur par défaut est \(\alpha = 1\).

from sklearn.linear_model import Lasso
regresseur3dLasso = Lasso(alpha=1.0).fit(X, y)  # par défaut alpha=1.0
print("Pente par rapport à la 1ère var. explicative :", pente)
print("Paramètres estimés par rapport aux 2 var. explicatives :", regresseur3dLasso.coef_)
print("Intercept :", regresseur3dLasso.intercept_)
print("R2 :", regresseur3dLasso.score(X, y))

Question :

Faire varier \(\alpha\) et examiner les résultats. A partir de quelle valeur de \(\alpha\) (approximativement) le paramètre associé à la seconde variable explicative est égal à 0 ?

Correction :

regresseur3dLasso = Lasso(alpha=1.5).fit(X, y)  # par défaut alpha=1.0
print("Pente par rapport à la 1ère var. explicative :", pente)
print("Paramètres estimés par rapport aux 2 var. explicatives :", regresseur3dLasso.coef_)
print("Intercept :", regresseur3dLasso.intercept_)
print("R2 :", regresseur3dLasso.score(X, y))

Régression linéaire sur les données ozone

Lecture des données en format CSV (comma separated values) dans un dataframe ozonedf. Affichage des projections des observations sur des paires de variables quantitatives, permettant de voir la distribution conjointe des variables de chaque paire. Sur la diagonale sont représentés les histogrammes des valeurs des variables correspondantes.

ozonedf = pd.read_csv('Dataset_ozone.csv', delimiter=';', decimal=',')
print(ozonedf.dtypes)  # afficher le codage de chaque variable
selected_cols = ['maxO3','T9','T12','T15','Vx9','Vx12','Vx15','maxO3v']
donneesOzoneSel1 = ozonedf[selected_cols]

%matplotlib inline

plt.figure()
pairplot = pd.plotting.scatter_matrix(donneesOzoneSel1, figsize=(20,20), marker = 'o', hist_kwds = {'bins': 14}, s = 50, alpha = 0.8)

Nous souhaitons d’abord savoir dans quelle mesure il est possible de prédire la variable maxO3 à partir de toutes les autres variables quantitatives. Les paramètres de ce modèle de régression linéaire peuvent être obtenus facilement :

npDonneesOzoneSel1 = donneesOzoneSel1.to_numpy()
X = npDonneesOzoneSel1[:, 1:8]
y = npDonneesOzoneSel1[:, 0]
regresseurOzone1 = LinearRegression().fit(X, y)
print("Paramètres estimés par rapport aux 7 var. explicatives :", regresseurOzone1.coef_)
print("Intercept :", regresseurOzone1.intercept_)
print("R2 :", regresseurOzone1.score(X, y))

Le modèle a d’assez bonnes performances, mais il n’est pas utilisable pour prédire la valeur de maxO3 car plusieurs des variables explicatives employées sont mesurées en même temps ou après que le pic d’ozone soit survenu.

Pour obtenir un odèle permettant éventuellement de faire des prédictions plusieurs heures avant le pic d’ozone, il faudrait exclure au moins les variables T15 et Vx15 :

#print(X[:3,:])
X = npDonneesOzoneSel1[:, np.r_[1, 2, 4, 5, 7]]  # np.r_[] permet de sélectionner une liste de colonnes
#print(X[:3,:])
y = npDonneesOzoneSel1[:, 0]
regresseurOzone2 = LinearRegression().fit(X, y)
print("Paramètres estimés par rapport aux 5 var. explicatives :", regresseurOzone2.coef_)
print("Intercept :", regresseurOzone2.intercept_)
print("R2 :", regresseurOzone2.score(X, y))

Le coefficient de détermination du nouveau modèle est presque aussi élevé que celui du modèle précédent. Pour évaluer plus directement les erreurs de prédiction, nous séparons les données disponibles en données d’apprentissage (sur lesquelles sont estimés les paramètres du modèle) et données de test (sur lesquelles le modèle est évalué) :

from sklearn.model_selection import train_test_split
# X_train, y_train : données sur lesquelles estimés les paramètres du modèle
# X_test, y_test : échantillon-test sur lequel on évalue les prédictions faites par le modèle
np.random.seed(1234567890)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.3) # test_size: proportion de l'échantillon-test
regresseurOzone3 = LinearRegression().fit(X_train, y_train)
print("R2 train :", regresseurOzone3.score(X_train, y_train))  # coefficient de détermination pour comparer
print("R2 test :", regresseurOzone3.score(X_test, y_test))  # coefficient de détermination sur échantillon-test

from sklearn.metrics import mean_squared_error
print("EQM test :", mean_squared_error(y_test, regresseurOzone3.predict(X_test), squared=False)) # erreur quadratique moyenne de prédiction sur échantillon-test

Diagnostic du modèle

Nous pouvons apliquer les techniques de diagnostic pour voir dans quelle mesure les hypothèses de modélisation sont respectées :

# Calcul des résidus

predictions = regresseurOzone2.predict(X)
residus = y - predictions

# Histogramme des résidus

fig = plt.figure()
plt.hist(residus, bins=10)
plt.title("Histogramme des résidus")
plt.show()

# Diagramme quantile-quantile des résidus

import statsmodels.api as sm

residusNormalises = (residus - np.mean(residus))/np.std(residus,ddof=1)
fig = sm.qqplot(residusNormalises, line='45')  # par défaut compare à loi normale d'espérance 0 et variance 1
plt.title("Diagramme quantile-quantile des résidus")
plt.show()

Question :

Affichez le graphique des résidus par rapport à la valeur prédite.

Correction :

# Résidus versus valeur prédite

fig = plt.figure()
plt.axis('equal')
plt.scatter(predictions, residus)
plt.axhline(y=0, color='r', linestyle='-')
plt.title("Résidus versus valeur prédite")
plt.show()

Ajout d’une variable nominale comme variable explicative

Jusqu’ici nous avons pris en compte dans le problème de régression uniquement les variables quantitatives, mais d’autres variables peuvent s’avérer utiles pour la prédiction. En effet, nous avions vu (analyse bivariée) que la présence ou l’absence de pluie avait un impact sur la concentration d’ozone :

data = [ozonedf[ozonedf['pluie']=='Pluie']["maxO3"], ozonedf[ozonedf['pluie']=='Sec']["maxO3"]]
fig, ax = plt.subplots()
ax.set_title('Distributions des concentrations O3 par temps de pluie et par temps sec :')
ax.boxplot(data, labels=['Pluie','Sec'])

Pour inclure les variables nominales, pluie et vent dans notre cas, dans la construction du modèle de régression, nous devons coder leurs valeurs de façon cohérente avec la nature du modèle recherché. Nous nous servirons ici d’une fonction de pandas qui permet d’encoder chaque variable nominale à \(k\) modalités sous la forme de \(k\) variables binaires (encodage one-hot) :

ozonedfEncode = pd.get_dummies(ozonedf, columns=["pluie","vent"])
ozonedfEncode.head()

Nous pouvons maintenant estimer le modèle en incluant les colonnes ajoutées (et en excluant les colonnes T15 et Vx15) :

selected_cols2 = ['maxO3','T9','T12','Vx9','Vx12','maxO3v','pluie_Pluie','pluie_Sec','vent_Est','vent_Nord','vent_Ouest','vent_Sud']
npDonneesOzoneSel2 = ozonedfEncode[selected_cols2].to_numpy()
X2 = npDonneesOzoneSel2[:, 1:]
y2 = npDonneesOzoneSel2[:, 0]

np.random.seed(1234567890)
X_train2, X_test2, y_train2, y_test2 = train_test_split(X2, y2, test_size = 0.3) # test_size: proportion de l'échantillon-test

regresseurOzone4 = LinearRegression().fit(X_train2, y_train2)
print("R2 train :", regresseurOzone4.score(X_train2, y_train2))  # coefficient de détermination sur train
print("R2 test :", regresseurOzone4.score(X_test2, y_test2))  # coefficient de détermination sur échantillon-test
print("EQM test :", mean_squared_error(y_test2, regresseurOzone4.predict(X_test2), squared=False)) # erreur quadratique moyenne de prédiction sur échantillon-test

Question :

Les résultats s’améliorent si la direction du vent est ignorée dans la construction du modèle et dans la prédiction ?

Correction :

selected_cols3 = ['maxO3','T9','T12','Vx9','Vx12','maxO3v','pluie_Pluie','pluie_Sec']
npDonneesOzoneSel3 = ozonedfEncode[selected_cols3].to_numpy()
X3 = npDonneesOzoneSel3[:, 1:]
y3 = npDonneesOzoneSel3[:, 0]

np.random.seed(1234567890)
X_train3, X_test3, y_train3, y_test3 = train_test_split(X3, y3, test_size = 0.3) # test_size: proportion de l'échantillon-test

regresseurOzone4 = LinearRegression().fit(X_train3, y_train3)
print("R2 train :", regresseurOzone4.score(X_train3, y_train3))  # coefficient de détermination sur train
print("R2 test :", regresseurOzone4.score(X_test3, y_test3))  # coefficient de détermination sur échantillon-test
print("EQM test :", mean_squared_error(y_test3, regresseurOzone4.predict(X_test3), squared=False)) # erreur quadratique moyenne de prédiction sur échantillon-test

Question :

Quels résultats sont obtenus avec une régression linéaire LASSO ?

Correction :

from sklearn.linear_model import Lasso

selected_cols2 = ['maxO3','T9','T12','Vx9','Vx12','maxO3v','pluie_Pluie','pluie_Sec','vent_Est','vent_Nord','vent_Ouest','vent_Sud']
npDonneesOzoneSel2 = ozonedfEncode[selected_cols2].to_numpy()
X2 = npDonneesOzoneSel2[:, 1:]
y2 = npDonneesOzoneSel2[:, 0]

np.random.seed(1234567890)
X_train2, X_test2, y_train2, y_test2 = train_test_split(X2, y2, test_size = 0.3) # test_size: proportion de l'échantillon-test

regresseur3dLasso2 = Lasso(alpha=0.5).fit(X_train2, y_train2)  # par défaut alpha=1.0
print("R2 train :", regresseur3dLasso2.score(X_train2, y_train2))  # coefficient de détermination sur train
print("R2 test :", regresseur3dLasso2.score(X_test2, y_test2))  # coefficient de détermination sur échantillon-test
print("EQM test :", mean_squared_error(y_test2, regresseur3dLasso2.predict(X_test2), squared=False)) # erreur quadratique moyenne de prédiction sur échantillon-test
print("Paramètres estimés par rapport aux 12 var. explicatives :", regresseur3dLasso2.coef_)