Régression linéaire

Régression linéaire#

Tous les modèles sont faux, mais certains sont utiles.

George E. P. Box

La régression linéaire est le modèle le plus fondamental de l’apprentissage supervisé. Malgré sa simplicité apparente, elle constitue la pierre angulaire sur laquelle reposent de nombreuses méthodes plus élaborées. Ce chapitre présente la régression linéaire simple et multiple, ses hypothèses, les méthodes de diagnostic, la régularisation et les métriques d’évaluation, le tout illustré par un exemple complet avec Scikit-learn.

Introduction#

En apprentissage supervisé, on distingue deux grandes familles de problèmes selon la nature de la variable cible \(y\) :

  • Régression : la variable cible est continue (prix d’un bien immobilier, température, chiffre d’affaires).

  • Classification : la variable cible est discrète (catégorie d’un email : spam ou non, espèce d’une fleur).

La régression cherche à modéliser la relation entre une ou plusieurs variables explicatives \(x_1, \ldots, x_p\) et une variable cible \(y \in \mathbb{R}\). Le modèle le plus simple pour exprimer cette relation est le modèle linéaire.

Remarque 47

La régression linéaire ne se limite pas aux relations « en ligne droite ». Grâce aux transformations polynomiales, aux interactions et aux fonctions de base, un modèle linéaire dans ses paramètres peut capturer des relations fortement non linéaires dans l’espace des features. La linéarité du modèle porte sur les coefficients \(\boldsymbol{\beta}\), pas sur les variables d’entrée.

Régression linéaire simple#

La régression linéaire simple modélise la relation entre une variable explicative \(x\) et une variable cible \(y\).

Le modèle#

Définition 53 (Modèle de régression linéaire simple)

Soit \((x_1, y_1), \ldots, (x_n, y_n)\) un échantillon de \(n\) observations. Le modèle de régression linéaire simple postule que

\[y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \quad i = 1, \ldots, n\]

où :

  • \(\beta_0 \in \mathbb{R}\) est l”ordonnée à l’origine (intercept) : la valeur prédite de \(y\) lorsque \(x = 0\)

  • \(\beta_1 \in \mathbb{R}\) est la pente : la variation de \(y\) associée à une augmentation unitaire de \(x\)

  • \(\varepsilon_i\) est le terme d’erreur (ou résidu théorique), supposé aléatoire, de moyenne nulle et de variance constante \(\sigma^2\)

Remarque 48

Le terme d’erreur \(\varepsilon_i\) capture tout ce que le modèle linéaire ne peut expliquer : l’effet de variables non incluses, le bruit de mesure, et toute non-linéarité résiduelle. C’est l’écart entre la valeur observée \(y_i\) et la valeur prédite \(\hat{y}_i = \beta_0 + \beta_1 x_i\).

Méthode des moindres carrés#

L’estimation des paramètres \(\beta_0\) et \(\beta_1\) se fait en minimisant la somme des carrés des résidus (Residual Sum of Squares, RSS).

Définition 54 (Moindres carrés ordinaires (MCO))

Les estimateurs des moindres carrés ordinaires (MCO, ou Ordinary Least Squares, OLS) sont les valeurs \(\hat{\beta}_0, \hat{\beta}_1\) qui minimisent

\[\text{RSS}(\beta_0, \beta_1) = \sum_{i=1}^n (y_i - \beta_0 - \beta_1 x_i)^2\]

En annulant les dérivées partielles par rapport à \(\beta_0\) et \(\beta_1\), on obtient les formules closes :

\[\hat{\beta}_1 = \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^n (x_i - \bar{x})^2} = \frac{\text{Cov}(x, y)}{\text{Var}(x)}\]
\[\hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}\]

où \(\bar{x} = \frac{1}{n}\sum_{i=1}^n x_i\) et \(\bar{y} = \frac{1}{n}\sum_{i=1}^n y_i\) sont les moyennes empiriques.

Proof. On pose \(f(\beta_0, \beta_1) = \sum_{i=1}^n (y_i - \beta_0 - \beta_1 x_i)^2\). Les conditions nécessaires d’optimalité (annulation du gradient) donnent le système :

\[\frac{\partial f}{\partial \beta_0} = -2 \sum_{i=1}^n (y_i - \beta_0 - \beta_1 x_i) = 0\]
\[\frac{\partial f}{\partial \beta_1} = -2 \sum_{i=1}^n x_i (y_i - \beta_0 - \beta_1 x_i) = 0\]

De la première équation, on tire \(\sum y_i = n\beta_0 + \beta_1 \sum x_i\), soit \(\hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}\).

En substituant dans la seconde équation et en développant, on obtient

\[\hat{\beta}_1 = \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^n (x_i - \bar{x})^2}\]

La matrice hessienne est définie positive (car \(\sum (x_i - \bar{x})^2 > 0\) dès que les \(x_i\) ne sont pas tous égaux), donc le point critique est bien un minimum global.

Interprétation géométrique#

Proposition 10 (Projection orthogonale)

La solution des moindres carrés \(\hat{\mathbf{y}} = \hat{\beta}_0 \mathbf{1} + \hat{\beta}_1 \mathbf{x}\) est la projection orthogonale du vecteur \(\mathbf{y} \in \mathbb{R}^n\) sur le sous-espace vectoriel \(\text{Im}(\mathbf{X})\) engendré par les colonnes de la matrice de design \(\mathbf{X} = [\mathbf{1} \mid \mathbf{x}]\).

Le vecteur des résidus \(\hat{\boldsymbol{\varepsilon}} = \mathbf{y} - \hat{\mathbf{y}}\) est orthogonal à \(\text{Im}(\mathbf{X})\) :

\[\mathbf{X}^T \hat{\boldsymbol{\varepsilon}} = \mathbf{X}^T (\mathbf{y} - \hat{\mathbf{y}}) = \mathbf{0}\]

Autrement dit, les résidus sont orthogonaux à chaque colonne de \(\mathbf{X}\). C’est l’interprétation géométrique du théorème de la projection.

Remarque 49

Cette interprétation géométrique est fondamentale. Parmi tous les vecteurs de \(\text{Im}(\mathbf{X})\), celui qui est le plus proche de \(\mathbf{y}\) (au sens de la norme euclidienne) est la projection orthogonale \(\hat{\mathbf{y}}\). La régression linéaire par moindres carrés est donc un problème d’approximation au sens des moindres carrés dans un espace vectoriel de dimension finie.

Hide code cell source

 
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats

sns.set_theme(style="whitegrid", palette="muted", font_scale=1.1)

Hide code cell source

 
# Illustration de la régression linéaire simple
np.random.seed(42)
n = 50
x = np.random.uniform(0, 10, n)
y = 2.5 * x + 3 + np.random.randn(n) * 3

# Estimation manuelle des coefficients
x_bar = x.mean()
y_bar = y.mean()
beta_1_hat = np.sum((x - x_bar) * (y - y_bar)) / np.sum((x - x_bar) ** 2)
beta_0_hat = y_bar - beta_1_hat * x_bar

print(f"Coefficients estimés (calcul manuel) :")
print(f"  β₀ = {beta_0_hat:.4f}  (valeur théorique : 3)")
print(f"  β₁ = {beta_1_hat:.4f}  (valeur théorique : 2.5)")

Coefficients estimés (calcul manuel) :
  β₀ = 3.2901  (valeur théorique : 3)
  β₁ = 2.4330  (valeur théorique : 2.5)

Hide code cell source

 
# Visualisation de la droite de régression
fig, axes = plt.subplots(2, 1, figsize=(9, 9))

# Graphique 1 : droite de régression avec résidus
x_line = np.linspace(0, 10, 100)
y_line = beta_0_hat + beta_1_hat * x_line
y_pred = beta_0_hat + beta_1_hat * x

axes[0].scatter(x, y, alpha=0.7, edgecolors='k', linewidths=0.5, zorder=3, label='Observations')
axes[0].plot(x_line, y_line, color='red', linewidth=2, label=f'$\\hat{{y}} = {beta_0_hat:.2f} + {beta_1_hat:.2f}x$')

# Résidus
for i in range(n):
    axes[0].plot([x[i], x[i]], [y[i], y_pred[i]], color='gray', alpha=0.4, linewidth=0.8)

axes[0].set_xlabel('$x$')
axes[0].set_ylabel('$y$')
axes[0].set_title('Régression linéaire simple — MCO')
axes[0].legend()

# Graphique 2 : résidus vs valeurs prédites
residus = y - y_pred
axes[1].scatter(y_pred, residus, alpha=0.7, edgecolors='k', linewidths=0.5)
axes[1].axhline(y=0, color='red', linestyle='--', linewidth=1)
axes[1].set_xlabel('Valeurs prédites $\\hat{y}$')
axes[1].set_ylabel('Résidus $y - \\hat{y}$')
axes[1].set_title('Résidus vs valeurs prédites')

plt.tight_layout()
plt.show()
_images/b96a3ad7738a6eae8ba5cdb7a963d83ae55e407190d5c43f5a267030d0ff1798.png

Régression linéaire multiple#

La régression linéaire multiple généralise le modèle simple à \(p\) variables explicatives.

Le modèle#

Définition 55 (Modèle de régression linéaire multiple)

Soit un échantillon de \(n\) observations avec \(p\) variables explicatives. Le modèle de régression linéaire multiple s’écrit

\[y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_p x_{ip} + \varepsilon_i, \quad i = 1, \ldots, n\]

Sous forme matricielle :

\[\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\]

où :

  • \(\mathbf{y} = (y_1, \ldots, y_n)^T \in \mathbb{R}^n\) est le vecteur des réponses

  • \(\mathbf{X} \in \mathbb{R}^{n \times (p+1)}\) est la matrice de design, dont la première colonne est un vecteur de \(1\) (pour l’intercept) et les colonnes suivantes contiennent les valeurs des variables explicatives

  • \(\boldsymbol{\beta} = (\beta_0, \beta_1, \ldots, \beta_p)^T \in \mathbb{R}^{p+1}\) est le vecteur des paramètres

  • \(\boldsymbol{\varepsilon} = (\varepsilon_1, \ldots, \varepsilon_n)^T \in \mathbb{R}^n\) est le vecteur des erreurs

Solution des moindres carrés#

Proposition 11 (Estimateur des moindres carrés (cas multiple))

Si la matrice \(\mathbf{X}^T\mathbf{X}\) est inversible (c’est-à-dire si \(\mathbf{X}\) est de rang plein \(p + 1\)), l’estimateur MCO est donné par

\[\hat{\boldsymbol{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}\]

Le vecteur des valeurs prédites est

\[\hat{\mathbf{y}} = \mathbf{X}\hat{\boldsymbol{\beta}} = \mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y} = \mathbf{H}\mathbf{y}\]

où \(\mathbf{H} = \mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\) est la matrice chapeau (hat matrix), qui projette orthogonalement \(\mathbf{y}\) sur l’espace colonne de \(\mathbf{X}\).

Proof. On minimise \(\|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 = (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})\). En développant et en dérivant par rapport à \(\boldsymbol{\beta}\) :

\[\frac{\partial}{\partial \boldsymbol{\beta}} \|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 = -2\mathbf{X}^T(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})\]

L’annulation du gradient donne les équations normales :

\[\mathbf{X}^T\mathbf{X}\boldsymbol{\beta} = \mathbf{X}^T\mathbf{y}\]

Si \(\mathbf{X}^T\mathbf{X}\) est inversible, on obtient \(\hat{\boldsymbol{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}\). La matrice \(\mathbf{X}^T\mathbf{X}\) est semi-définie positive ; elle est définie positive si et seulement si \(\mathbf{X}\) est de rang plein, ce qui garantit l’unicité de la solution.

Remarque 50

En pratique, on n’inverse jamais directement \(\mathbf{X}^T\mathbf{X}\) (numériquement instable et coûteux). On résout les équations normales par décomposition de Cholesky, ou mieux, on utilise la décomposition QR de \(\mathbf{X}\) : si \(\mathbf{X} = \mathbf{Q}\mathbf{R}\), alors \(\hat{\boldsymbol{\beta}} = \mathbf{R}^{-1}\mathbf{Q}^T\mathbf{y}\). Scikit-learn utilise cette approche en interne.

Hide code cell source

 
# Illustration de la régression multiple avec NumPy
np.random.seed(42)
n = 100

x1 = np.random.uniform(0, 10, n)
x2 = np.random.uniform(0, 5, n)
y_multi = 3 + 2 * x1 - 1.5 * x2 + np.random.randn(n) * 2

# Matrice de design (avec colonne de 1 pour l'intercept)
X_design = np.column_stack([np.ones(n), x1, x2])

# Solution MCO : β = (X^T X)^{-1} X^T y
beta_hat = np.linalg.solve(X_design.T @ X_design, X_design.T @ y_multi)

print("Coefficients estimés (calcul matriciel) :")
print(f"  β₀ = {beta_hat[0]:.4f}  (valeur théorique : 3)")
print(f"  β₁ = {beta_hat[1]:.4f}  (valeur théorique : 2)")
print(f"  β₂ = {beta_hat[2]:.4f}  (valeur théorique : -1.5)")

Coefficients estimés (calcul matriciel) :
  β₀ = 2.8212  (valeur théorique : 3)
  β₁ = 1.9317  (valeur théorique : 2)
  β₂ = -1.2123  (valeur théorique : -1.5)

Hide code cell source

 
# Vérification avec Scikit-learn
from sklearn.linear_model import LinearRegression

X_multi = np.column_stack([x1, x2])
lr = LinearRegression()
lr.fit(X_multi, y_multi)

print("Coefficients estimés (Scikit-learn) :")
print(f"  intercept = {lr.intercept_:.4f}")
print(f"  coef      = {lr.coef_}")

Coefficients estimés (Scikit-learn) :
  intercept = 2.8212
  coef      = [ 1.93165494 -1.21227544]

Hypothèses du modèle#

Le modèle de régression linéaire repose sur un ensemble d’hypothèses. Leur vérification est essentielle pour s’assurer de la validité des inférences (intervalles de confiance, tests d’hypothèses) et de la qualité des prédictions.

Définition 56 (Hypothèses de Gauss-Markov)

Le modèle \(\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\) repose sur les hypothèses suivantes :

  1. Linéarité : la relation entre \(\mathbf{y}\) et \(\mathbf{X}\) est linéaire dans les paramètres \(\boldsymbol{\beta}\)

  2. Exogénéité : \(\mathbb{E}[\boldsymbol{\varepsilon} \mid \mathbf{X}] = \mathbf{0}\) (les erreurs sont de moyenne nulle conditionnellement aux régresseurs)

  3. Homoscédasticité : \(\text{Var}(\varepsilon_i \mid \mathbf{X}) = \sigma^2\) pour tout \(i\) (variance constante des erreurs)

  4. Absence d’autocorrélation : \(\text{Cov}(\varepsilon_i, \varepsilon_j \mid \mathbf{X}) = 0\) pour \(i \neq j\) (indépendance des erreurs)

  5. Rang plein : \(\text{rang}(\mathbf{X}) = p + 1\) (pas de multicolinéarité parfaite)

Sous ces hypothèses, le théorème de Gauss-Markov garantit que l’estimateur MCO est le meilleur estimateur linéaire non biaisé (Best Linear Unbiased Estimator, BLUE).

Remarque 51

On ajoute parfois une sixième hypothèse, la normalité des erreurs : \(\boldsymbol{\varepsilon} \sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I}_n)\). Cette hypothèse n’est pas nécessaire pour la propriété BLUE, mais elle est requise pour construire des intervalles de confiance exacts et des tests de Student sur les coefficients.

Proposition 12 (Conséquences de la violation des hypothèses)

Hypothèse violée

Conséquence

Non-linéarité

Biais systématique, sous-apprentissage

Hétéroscédasticité

MCO non optimal, erreurs standard incorrectes

Autocorrélation

Erreurs standard sous-estimées, tests invalides

Multicolinéarité

Coefficients instables, grande variance des estimateurs

Non-normalité des résidus

Tests et intervalles de confiance approximatifs

Diagnostic des résidus#

Le diagnostic des résidus est la vérification a posteriori des hypothèses du modèle. On analyse les résidus \(\hat{\varepsilon}_i = y_i - \hat{y}_i\) pour détecter d’éventuelles violations.

Définition 57 (Résidus standardisés)

Les résidus standardisés (ou résidus studentisés internes) sont définis par

\[r_i = \frac{\hat{\varepsilon}_i}{\hat{\sigma}\sqrt{1 - h_{ii}}}\]

où \(\hat{\sigma}\) est l’estimateur de l’écart-type résiduel et \(h_{ii}\) est le \(i\)-ème élément diagonal de la matrice chapeau \(\mathbf{H}\). La correction par \(\sqrt{1 - h_{ii}}\) tient compte du fait que les résidus n’ont pas tous la même variance, même sous homoscédasticité.

Hide code cell source

 
# Génération de données pour le diagnostic
np.random.seed(42)
n = 200
x = np.random.uniform(0, 10, n)

# Cas 1 : hypothèses respectées
y_ok = 2 * x + 5 + np.random.randn(n) * 2

# Cas 2 : hétéroscédasticité
y_hetero = 2 * x + 5 + np.random.randn(n) * (0.5 + 0.5 * x)

# Cas 3 : non-linéarité
y_nonlin = 0.3 * x**2 - x + 5 + np.random.randn(n) * 2

Hide code cell source

 
from sklearn.linear_model import LinearRegression

fig, axes = plt.subplots(2, 3, figsize=(16, 10))

cas = [("Hypothèses respectées", y_ok),
       ("Hétéroscédasticité", y_hetero),
       ("Non-linéarité", y_nonlin)]

for j, (titre, y_cas) in enumerate(cas):
    lr = LinearRegression()
    lr.fit(x.reshape(-1, 1), y_cas)
    y_pred = lr.predict(x.reshape(-1, 1))
    residus = y_cas - y_pred

    # Résidus vs valeurs prédites
    axes[0, j].scatter(y_pred, residus, alpha=0.5, s=15, edgecolors='k', linewidths=0.3)
    axes[0, j].axhline(y=0, color='red', linestyle='--', linewidth=1)
    axes[0, j].set_xlabel('Valeurs prédites $\\hat{y}$')
    axes[0, j].set_ylabel('Résidus')
    axes[0, j].set_title(f'{titre}\nRésidus vs prédites')

    # QQ-plot
    stats.probplot(residus, dist="norm", plot=axes[1, j])
    axes[1, j].set_title(f'{titre}\nQQ-plot')
    axes[1, j].get_lines()[0].set_color("steelblue")

plt.tight_layout()
plt.show()
_images/8f2a3f2884abf0d5cc3b7559f878090b7124a99f136f7d00901f5564fd1b6a10.png

Remarque 52

Les graphiques de diagnostic se lisent comme suit :

  • Résidus vs prédites : on attend un nuage sans structure visible, centré sur zéro, de largeur constante. Un motif en « entonnoir » indique de l’hétéroscédasticité ; une forme en « U » ou « arc » indique de la non-linéarité.

  • QQ-plot : on attend un alignement sur la diagonale. Des écarts dans les queues signalent une distribution non normale (queues lourdes ou légères).

Hide code cell source

 
# Tests statistiques de diagnostic
from scipy.stats import shapiro, normaltest

for titre, y_cas in cas:
    lr = LinearRegression()
    lr.fit(x.reshape(-1, 1), y_cas)
    residus = y_cas - lr.predict(x.reshape(-1, 1))

    stat_sw, p_sw = shapiro(residus)
    stat_dp, p_dp = normaltest(residus)

    print(f"{titre}")
    print(f"  Shapiro-Wilk       : W = {stat_sw:.4f}, p = {p_sw:.4f}")
    print(f"  D'Agostino-Pearson : S = {stat_dp:.4f}, p = {p_dp:.4f}")
    print(f"  → {'Normal' if p_sw > 0.05 else 'Non normal'} (Shapiro-Wilk, α = 0.05)")
    print()

Hypothèses respectées
  Shapiro-Wilk       : W = 0.9879, p = 0.0862
  D'Agostino-Pearson : S = 7.0282, p = 0.0298
  → Normal (Shapiro-Wilk, α = 0.05)

Hétéroscédasticité
  Shapiro-Wilk       : W = 0.9669, p = 0.0001
  D'Agostino-Pearson : S = 19.4711, p = 0.0001
  → Non normal (Shapiro-Wilk, α = 0.05)

Non-linéarité
  Shapiro-Wilk       : W = 0.9860, p = 0.0445
  D'Agostino-Pearson : S = 4.5133, p = 0.1047
  → Non normal (Shapiro-Wilk, α = 0.05)

Régularisation#

Lorsque le nombre de variables est élevé, que les variables sont corrélées entre elles (multicolinéarité) ou que l’on souhaite prévenir le surapprentissage, on peut régulariser le modèle en ajoutant une pénalité sur la norme des coefficients à la fonction de coût.

Régression Ridge (\(L_2\))#

Définition 58 (Régression Ridge)

La régression Ridge (ou régression de Tikhonov) ajoute une pénalité \(L_2\) au problème des moindres carrés :

\[\hat{\boldsymbol{\beta}}_{\text{Ridge}} = \arg\min_{\boldsymbol{\beta}} \left\{ \|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|_2^2 + \alpha \|\boldsymbol{\beta}\|_2^2 \right\}\]

où \(\alpha > 0\) est le paramètre de régularisation. La solution admet une forme close :

\[\hat{\boldsymbol{\beta}}_{\text{Ridge}} = (\mathbf{X}^T\mathbf{X} + \alpha \mathbf{I})^{-1}\mathbf{X}^T\mathbf{y}\]

L’ajout de \(\alpha \mathbf{I}\) rend la matrice toujours inversible, même en cas de multicolinéarité.

Remarque 53

Ridge ne met jamais de coefficient exactement à zéro : il les réduit (shrinks) vers zéro. Plus \(\alpha\) est grand, plus les coefficients sont petits, mais aucun n’est annulé. Ridge est donc un outil de régularisation, pas de sélection de variables.

Régression Lasso (\(L_1\))#

Définition 59 (Régression Lasso)

La régression Lasso (Least Absolute Shrinkage and Selection Operator) ajoute une pénalité \(L_1\) :

\[\hat{\boldsymbol{\beta}}_{\text{Lasso}} = \arg\min_{\boldsymbol{\beta}} \left\{ \frac{1}{2n}\|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|_2^2 + \alpha \|\boldsymbol{\beta}\|_1 \right\}\]

où \(\|\boldsymbol{\beta}\|_1 = \sum_{j=1}^p |\beta_j|\) est la norme \(L_1\).

Contrairement à Ridge, Lasso peut mettre certains coefficients exactement à zéro, réalisant ainsi une sélection automatique de variables.

Remarque 54

La propriété de parcimonie du Lasso provient de la géométrie de la boule \(L_1\) : ses « coins » sont alignés sur les axes, ce qui favorise les solutions où certaines coordonnées sont nulles. Géométriquement, la contrainte \(\|\boldsymbol{\beta}\|_1 \leq t\) est un polyèdre (losange en 2D) dont les sommets se trouvent sur les axes.

ElasticNet#

Définition 60 (ElasticNet)

Le ElasticNet combine les pénalités \(L_1\) et \(L_2\) :

\[\hat{\boldsymbol{\beta}}_{\text{EN}} = \arg\min_{\boldsymbol{\beta}} \left\{ \frac{1}{2n}\|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|_2^2 + \alpha \left( \rho \|\boldsymbol{\beta}\|_1 + \frac{1 - \rho}{2} \|\boldsymbol{\beta}\|_2^2 \right) \right\}\]

où \(\rho \in [0, 1]\) est le ratio de mélange (l1_ratio) :

  • \(\rho = 1\) : Lasso pur

  • \(\rho = 0\) : Ridge pur

  • \(0 < \rho < 1\) : compromis entre sélection de variables (L1) et stabilité (L2)

Remarque 55

ElasticNet corrige deux limitations du Lasso :

  1. Lorsque \(p > n\) (plus de variables que d’observations), Lasso sélectionne au plus \(n\) variables. ElasticNet peut en sélectionner davantage.

  2. En présence de variables fortement corrélées, Lasso a tendance à n’en sélectionner qu’une et à ignorer les autres. ElasticNet les conserve en groupe grâce à la composante Ridge.

Hide code cell source

 
# Comparaison Ridge, Lasso, ElasticNet
from sklearn.linear_model import Ridge, Lasso, ElasticNet, LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split

np.random.seed(42)
n, p = 100, 20
X_reg = np.random.randn(n, p)
# Seules les 5 premières variables sont informatives
beta_true = np.zeros(p)
beta_true[:5] = [3, -2, 1.5, -1, 0.5]
y_reg = X_reg @ beta_true + np.random.randn(n) * 2

X_train, X_test, y_train, y_test = train_test_split(X_reg, y_reg, test_size=0.3, random_state=42)

# Standardisation
scaler = StandardScaler()
X_train_s = scaler.fit_transform(X_train)
X_test_s = scaler.transform(X_test)

# Modèles
models = {
    'MCO': LinearRegression(),
    'Ridge (α=1)': Ridge(alpha=1),
    'Lasso (α=0.1)': Lasso(alpha=0.1),
    'ElasticNet (α=0.1, ρ=0.5)': ElasticNet(alpha=0.1, l1_ratio=0.5),
}

fig, ax = plt.subplots(figsize=(12, 6))
bar_width = 0.2
indices = np.arange(p)

for k, (name, model) in enumerate(models.items()):
    model.fit(X_train_s, y_train)
    coefs = model.coef_
    ax.bar(indices + k * bar_width, coefs, bar_width, label=name, alpha=0.8)

ax.axhline(y=0, color='black', linewidth=0.5)
ax.set_xlabel('Index du coefficient')
ax.set_ylabel('Valeur du coefficient')
ax.set_title('Comparaison des coefficients estimés')
ax.set_xticks(indices + 1.5 * bar_width)
ax.set_xticklabels([f'$\\beta_{{{j}}}$' for j in range(p)], fontsize=8)
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()
_images/d3d817aa3dd15cc56e9f4b346aeb6a5d1d971e234d29c51d15b8e557ec951869.png

Hide code cell source

 
# Chemin de régularisation du Lasso
from sklearn.linear_model import lasso_path

alphas, coefs_path, _ = lasso_path(X_train_s, y_train, alphas=np.logspace(-3, 1, 100))

fig, ax = plt.subplots(figsize=(10, 6))
for j in range(p):
    label = f'$\\beta_{{{j}}}$' if j < 5 else None
    style = '-' if j < 5 else '--'
    alpha_val = 0.9 if j < 5 else 0.2
    ax.plot(np.log10(alphas), coefs_path[j], linestyle=style, alpha=alpha_val, label=label)

ax.set_xlabel('$\\log_{10}(\\alpha)$')
ax.set_ylabel('Coefficient')
ax.set_title('Chemin de régularisation — Lasso')
ax.legend(loc='upper right', fontsize=9)
ax.axhline(y=0, color='black', linewidth=0.5)
plt.tight_layout()
plt.show()
_images/f08eee3b89dd46ea8342f2573e6feac178da818888b3496259ff17e00d5f0f22.png

Métriques d’évaluation#

L’évaluation d’un modèle de régression repose sur des métriques quantifiant l’écart entre les valeurs observées \(y_i\) et les valeurs prédites \(\hat{y}_i\).

Définition 61 (Erreur quadratique moyenne (MSE))

L”erreur quadratique moyenne (Mean Squared Error) est

\[\text{MSE} = \frac{1}{n}\sum_{i=1}^n (y_i - \hat{y}_i)^2\]

Elle pénalise fortement les grandes erreurs (à cause du carré). Son unité est le carré de l’unité de \(y\).

Définition 62 (Racine de l’erreur quadratique moyenne (RMSE))

La RMSE (Root Mean Squared Error) est la racine carrée de la MSE :

\[\text{RMSE} = \sqrt{\text{MSE}} = \sqrt{\frac{1}{n}\sum_{i=1}^n (y_i - \hat{y}_i)^2}\]

Elle a la même unité que \(y\), ce qui la rend plus interprétable que la MSE.

Définition 63 (Erreur absolue moyenne (MAE))

L”erreur absolue moyenne (Mean Absolute Error) est

\[\text{MAE} = \frac{1}{n}\sum_{i=1}^n |y_i - \hat{y}_i|\]

Elle est plus robuste aux valeurs aberrantes que la MSE, car elle ne met pas les erreurs au carré.

Définition 64 (Coefficient de détermination (\(R^2\)))

Le coefficient de détermination \(R^2\) mesure la proportion de la variance de \(y\) expliquée par le modèle :

\[R^2 = 1 - \frac{\text{RSS}}{\text{TSS}} = 1 - \frac{\sum_{i=1}^n (y_i - \hat{y}_i)^2}{\sum_{i=1}^n (y_i - \bar{y})^2}\]

où \(\text{TSS} = \sum(y_i - \bar{y})^2\) est la somme totale des carrés.

  • \(R^2 = 1\) : le modèle explique parfaitement la variance de \(y\)

  • \(R^2 = 0\) : le modèle ne fait pas mieux que prédire la moyenne \(\bar{y}\)

  • \(R^2 < 0\) : le modèle est pire que la prédiction par la moyenne (possible sur le test set)

Définition 65 (\(R^2\) ajusté)

Le \(R^2\) ajusté corrige le \(R^2\) pour le nombre de variables \(p\) :

\[R^2_{\text{adj}} = 1 - \frac{n - 1}{n - p - 1}(1 - R^2)\]

Contrairement au \(R^2\) brut, le \(R^2\) ajusté peut diminuer lorsqu’on ajoute une variable non informative, car la pénalisation \(\frac{n-1}{n-p-1}\) augmente avec \(p\).

Remarque 56

Le \(R^2\) augmente mécaniquement lorsqu’on ajoute des variables au modèle, même si elles sont inutiles. C’est pourquoi il ne doit jamais être utilisé seul pour comparer des modèles avec des nombres de variables différents. Le \(R^2\) ajusté ou les critères d’information (AIC, BIC) sont plus appropriés.

Hide code cell source

 
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score

# Exemple avec données synthétiques
np.random.seed(42)
n = 100
x = np.random.uniform(0, 10, n)
y_true = 2 * x + 5 + np.random.randn(n) * 3

lr = LinearRegression()
lr.fit(x.reshape(-1, 1), y_true)
y_pred = lr.predict(x.reshape(-1, 1))

mse = mean_squared_error(y_true, y_pred)
rmse = np.sqrt(mse)
mae = mean_absolute_error(y_true, y_pred)
r2 = r2_score(y_true, y_pred)

# R² ajusté
p = 1
r2_adj = 1 - (n - 1) / (n - p - 1) * (1 - r2)

print("Métriques d'évaluation :")
print(f"  MSE      = {mse:.4f}")
print(f"  RMSE     = {rmse:.4f}")
print(f"  MAE      = {mae:.4f}")
print(f"  R²       = {r2:.4f}")
print(f"  R² ajusté = {r2_adj:.4f}")

Métriques d'évaluation :
  MSE      = 7.2593
  RMSE     = 2.6943
  MAE      = 2.1031
  R²       = 0.8071
  R² ajusté = 0.8052

Hide code cell source

 
# Visualisation des métriques
fig, axes = plt.subplots(2, 1, figsize=(9, 9))

# Prédictions vs observations
axes[0].scatter(y_true, y_pred, alpha=0.6, edgecolors='k', linewidths=0.5)
lim_min = min(y_true.min(), y_pred.min()) - 1
lim_max = max(y_true.max(), y_pred.max()) + 1
axes[0].plot([lim_min, lim_max], [lim_min, lim_max], 'r--', linewidth=1.5, label='Prédiction parfaite')
axes[0].set_xlabel('Valeurs observées $y$')
axes[0].set_ylabel('Valeurs prédites $\\hat{y}$')
axes[0].set_title(f'Prédictions vs observations ($R^2 = {r2:.3f}$)')
axes[0].legend()
axes[0].set_aspect('equal')

# Distribution des erreurs
erreurs = y_true - y_pred
axes[1].hist(erreurs, bins=20, edgecolor='white', alpha=0.8, color='steelblue')
axes[1].axvline(x=0, color='red', linestyle='--', linewidth=1.5)
axes[1].set_xlabel('Erreur ($y - \\hat{y}$)')
axes[1].set_ylabel('Fréquence')
axes[1].set_title(f'Distribution des erreurs (MAE = {mae:.2f}, RMSE = {rmse:.2f})')

plt.tight_layout()
plt.show()
_images/37124fd011412d4497664cd4542d141ad94ce10cbebc376f5cdb6420c8e5ea4f.png

Exemple complet avec Scikit-learn#

Mettons en pratique l’ensemble des concepts sur un jeu de données réel : California Housing. L’objectif est de prédire la valeur médiane des maisons dans un district californien à partir de caractéristiques sociodémographiques et géographiques.

Hide code cell source

 
from sklearn.datasets import fetch_california_housing

housing = fetch_california_housing(as_frame=True)
df = housing.frame

print(f"Dimensions : {df.shape}")
print(f"\nVariables : {list(df.columns)}")
print(f"\nAperçu :")
df.head()

Dimensions : (20640, 9)

Variables : ['MedInc', 'HouseAge', 'AveRooms', 'AveBedrms', 'Population', 'AveOccup', 'Latitude', 'Longitude', 'MedHouseVal']

Aperçu :
MedInc HouseAge AveRooms AveBedrms Population AveOccup Latitude Longitude MedHouseVal
0 8.3252 41.0 6.984127 1.023810 322.0 2.555556 37.88 -122.23 4.526
1 8.3014 21.0 6.238137 0.971880 2401.0 2.109842 37.86 -122.22 3.585
2 7.2574 52.0 8.288136 1.073446 496.0 2.802260 37.85 -122.24 3.521
3 5.6431 52.0 5.817352 1.073059 558.0 2.547945 37.85 -122.25 3.413
4 3.8462 52.0 6.281853 1.081081 565.0 2.181467 37.85 -122.25 3.422

Hide code cell source

 
df.describe().round(3)
MedInc HouseAge AveRooms AveBedrms Population AveOccup Latitude Longitude MedHouseVal
count 20640.000 20640.000 20640.000 20640.000 20640.000 20640.000 20640.000 20640.000 20640.000
mean 3.871 28.639 5.429 1.097 1425.477 3.071 35.632 -119.570 2.069
std 1.900 12.586 2.474 0.474 1132.462 10.386 2.136 2.004 1.154
min 0.500 1.000 0.846 0.333 3.000 0.692 32.540 -124.350 0.150
25% 2.563 18.000 4.441 1.006 787.000 2.430 33.930 -121.800 1.196
50% 3.535 29.000 5.229 1.049 1166.000 2.818 34.260 -118.490 1.797
75% 4.743 37.000 6.052 1.100 1725.000 3.282 37.710 -118.010 2.647
max 15.000 52.000 141.909 34.067 35682.000 1243.333 41.950 -114.310 5.000

Analyse exploratoire#

Hide code cell source

 
fig, axes = plt.subplots(2, 4, figsize=(18, 8))

for ax, col in zip(axes.flat, df.columns):
    ax.hist(df[col], bins=40, edgecolor='white', alpha=0.8, color='steelblue')
    ax.set_title(col, fontsize=11)
    ax.set_ylabel('Fréquence')

plt.suptitle('Distribution des variables — California Housing', fontsize=14, y=1.01)
plt.tight_layout()
plt.show()
_images/3dfedd87db5e8ab478233d0a9a971fc0334a32a24a351db575c32e05ad9daeab.png

Hide code cell source

 
# Matrice de corrélation
corr = df.corr()

fig, ax = plt.subplots(figsize=(10, 8))
sns.heatmap(corr, annot=True, fmt=".2f", cmap="coolwarm", center=0,
            square=True, linewidths=0.5, ax=ax)
ax.set_title('Matrice de corrélation — California Housing')
plt.tight_layout()
plt.show()
_images/8253db880f26b15a2cf03dd60e3da1ced3f8474b407805095d5a94545a23fc30.png

Exemple 5

On observe que MedInc (revenu médian) est la variable la plus corrélée à MedHouseVal (\(r \approx 0.69\)). Les variables AveRooms et AveBedrms sont fortement corrélées entre elles (\(r \approx 0.85\)), signe de multicolinéarité. La variable HouseAge est faiblement corrélée à la cible (\(r \approx 0.11\)), mais pourrait néanmoins contribuer en combinaison avec d’autres features.

Préparation des données#

Hide code cell source

 
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

X = df.drop(columns='MedHouseVal')
y = df['MedHouseVal']

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

scaler = StandardScaler()
X_train_s = scaler.fit_transform(X_train)
X_test_s = scaler.transform(X_test)

print(f"Entraînement : {X_train_s.shape[0]} observations")
print(f"Test         : {X_test_s.shape[0]} observations")

Entraînement : 16512 observations
Test         : 4128 observations

Entraînement et comparaison des modèles#

Hide code cell source

 
from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score

models = {
    'MCO': LinearRegression(),
    'Ridge (α=1)': Ridge(alpha=1),
    'Lasso (α=0.01)': Lasso(alpha=0.01),
    'ElasticNet (α=0.01, ρ=0.5)': ElasticNet(alpha=0.01, l1_ratio=0.5),
}

results = []
for name, model in models.items():
    model.fit(X_train_s, y_train)
    y_pred = model.predict(X_test_s)

    mse = mean_squared_error(y_test, y_pred)
    rmse = np.sqrt(mse)
    mae = mean_absolute_error(y_test, y_pred)
    r2 = r2_score(y_test, y_pred)

    results.append({'Modèle': name, 'MSE': mse, 'RMSE': rmse, 'MAE': mae, 'R²': r2})
    print(f"{name:35s}  RMSE = {rmse:.4f}  MAE = {mae:.4f}  R² = {r2:.4f}")

df_results = pd.DataFrame(results)

MCO                                  RMSE = 0.7456  MAE = 0.5332  R² = 0.5758
Ridge (α=1)                          RMSE = 0.7456  MAE = 0.5332  R² = 0.5758
Lasso (α=0.01)                       RMSE = 0.7404  MAE = 0.5353  R² = 0.5816
ElasticNet (α=0.01, ρ=0.5)           RMSE = 0.7416  MAE = 0.5341  R² = 0.5803

Hide code cell source

 
# Coefficients des différents modèles
fig, ax = plt.subplots(figsize=(12, 6))

feature_names = X.columns
x_pos = np.arange(len(feature_names))
bar_width = 0.2

for k, (name, model) in enumerate(models.items()):
    coefs = model.coef_
    ax.bar(x_pos + k * bar_width, coefs, bar_width, label=name, alpha=0.85)

ax.set_xlabel('Variable')
ax.set_ylabel('Coefficient')
ax.set_title('Coefficients estimés — California Housing')
ax.set_xticks(x_pos + 1.5 * bar_width)
ax.set_xticklabels(feature_names, rotation=45, ha='right')
ax.legend(fontsize=9, loc='lower left')
ax.axhline(y=0, color='black', linewidth=0.5)
plt.tight_layout()
plt.show()
_images/6f1cacb17f1df9050ebc4413b8d588d9abf0725f32dcfd125104a0052599276f.png

Diagnostic du modèle MCO#

Hide code cell source

 
# Diagnostic des résidus du modèle MCO
lr = models['MCO']
y_pred_lr = lr.predict(X_test_s)
residus = y_test.values - y_pred_lr

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. Résidus vs prédites
axes[0, 0].scatter(y_pred_lr, residus, alpha=0.3, s=8, c='steelblue')
axes[0, 0].axhline(y=0, color='red', linestyle='--', linewidth=1.5)
axes[0, 0].set_xlabel('Valeurs prédites $\\hat{y}$')
axes[0, 0].set_ylabel('Résidus')
axes[0, 0].set_title('Résidus vs valeurs prédites')

# 2. QQ-plot
stats.probplot(residus, dist="norm", plot=axes[0, 1])
axes[0, 1].set_title('QQ-plot des résidus')
axes[0, 1].get_lines()[0].set_color("steelblue")

# 3. Distribution des résidus
axes[1, 0].hist(residus, bins=50, edgecolor='white', alpha=0.8, color='steelblue', density=True)
x_kde = np.linspace(residus.min(), residus.max(), 200)
from scipy.stats import norm
mu, sigma = norm.fit(residus)
axes[1, 0].plot(x_kde, norm.pdf(x_kde, mu, sigma), 'r-', linewidth=2, label=f'$\\mathcal{{N}}({mu:.2f}, {sigma:.2f}^2)$')
axes[1, 0].set_xlabel('Résidu')
axes[1, 0].set_ylabel('Densité')
axes[1, 0].set_title('Distribution des résidus')
axes[1, 0].legend()

# 4. Prédictions vs observations
axes[1, 1].scatter(y_test, y_pred_lr, alpha=0.3, s=8, c='steelblue')
lim = [y_test.min() - 0.5, y_test.max() + 0.5]
axes[1, 1].plot(lim, lim, 'r--', linewidth=1.5)
axes[1, 1].set_xlabel('Valeurs observées $y$')
axes[1, 1].set_ylabel('Valeurs prédites $\\hat{y}$')
axes[1, 1].set_title('Prédictions vs observations')

plt.suptitle('Diagnostic du modèle MCO — California Housing', fontsize=14, y=1.01)
plt.tight_layout()
plt.show()
_images/8d91ccf10ac7fc0525a6474b651f1a783d7a6436aae441cbcb619299b5a68b53.png

Remarque 57

Le diagnostic sur California Housing révèle plusieurs violations des hypothèses :

  • Le graphique des résidus vs prédites montre une structure non aléatoire, avec une accumulation de résidus négatifs pour les grandes valeurs prédites. Cela s’explique par le plafonnement de la variable cible à \(5.0\) (les prix supérieurs à 500 000 $ ont été censurés).

  • Le QQ-plot montre des queues lourdes, indiquant une non-normalité des résidus.

  • Ces observations suggèrent qu’un modèle linéaire simple est insuffisant pour capturer pleinement la structure des données. Des transformations non linéaires ou des modèles plus flexibles seraient à envisager.

Sélection du paramètre de régularisation par validation croisée#

Hide code cell source

 
from sklearn.linear_model import RidgeCV, LassoCV

# Ridge avec validation croisée
alphas_ridge = np.logspace(-3, 3, 100)
ridge_cv = RidgeCV(alphas=alphas_ridge, cv=5)
ridge_cv.fit(X_train_s, y_train)
print(f"Ridge — meilleur α : {ridge_cv.alpha_:.4f}")
print(f"Ridge — R² test    : {ridge_cv.score(X_test_s, y_test):.4f}")

# Lasso avec validation croisée
lasso_cv = LassoCV(cv=5, random_state=42, max_iter=10000)
lasso_cv.fit(X_train_s, y_train)
print(f"\nLasso — meilleur α : {lasso_cv.alpha_:.4f}")
print(f"Lasso — R² test    : {lasso_cv.score(X_test_s, y_test):.4f}")
print(f"Lasso — coefficients non nuls : {np.sum(lasso_cv.coef_ != 0)} / {len(lasso_cv.coef_)}")

Ridge — meilleur α : 0.0010
Ridge — R² test    : 0.5758

Lasso — meilleur α : 0.0008
Lasso — R² test    : 0.5766
Lasso — coefficients non nuls : 8 / 8

Hide code cell source

 
# Visualisation du chemin de régularisation Lasso avec CV
fig, ax = plt.subplots(figsize=(10, 6))

alphas_lasso = lasso_cv.alphas_
mse_mean = lasso_cv.mse_path_.mean(axis=1)
mse_std = lasso_cv.mse_path_.std(axis=1)

ax.semilogx(alphas_lasso, mse_mean, color='steelblue', linewidth=2)
ax.fill_between(alphas_lasso, mse_mean - mse_std, mse_mean + mse_std,
                alpha=0.2, color='steelblue')
ax.axvline(lasso_cv.alpha_, color='red', linestyle='--',
           label=f'$\\alpha^* = {lasso_cv.alpha_:.4f}$')
ax.set_xlabel('$\\alpha$ (paramètre de régularisation)')
ax.set_ylabel('MSE (validation croisée)')
ax.set_title('Sélection de $\\alpha$ par validation croisée — Lasso')
ax.legend()
plt.tight_layout()
plt.show()
_images/273ff880b02b2b6168178ad0290339867087246bc25b5fe748c90fc4d0cfdbf3.png

Résumé des résultats#

Hide code cell source

 
print(df_results.to_string(index=False))

                    Modèle      MSE     RMSE      MAE       R²
                       MCO 0.555892 0.745581 0.533200 0.575788
               Ridge (α=1) 0.555855 0.745557 0.533193 0.575816
            Lasso (α=0.01) 0.548255 0.740442 0.535326 0.581615
ElasticNet (α=0.01, ρ=0.5) 0.549953 0.741588 0.534077 0.580319

Remarque 58

Sur le jeu California Housing, les quatre modèles donnent des performances très proches. Cela s’explique par le fait que toutes les variables sont informatives et que le nombre d’observations (\(n \approx 20\,000\)) est largement supérieur au nombre de variables (\(p = 8\)). La régularisation apporte un bénéfice principalement lorsque \(p\) est grand relativement à \(n\), ou en présence de forte multicolinéarité. Dans le cas présent, la régularisation ne dégrade pas les performances, ce qui confirme la robustesse des modèles pénalisés.

Résumé#

Ce chapitre a couvert les fondements de la régression linéaire :

Concept

Formule / idée clé

Modèle simple

\(y = \beta_0 + \beta_1 x + \varepsilon\)

Modèle multiple

\(\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\)

Solution MCO

\(\hat{\boldsymbol{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}\)

Ridge (\(L_2\))

\(\hat{\boldsymbol{\beta}} = (\mathbf{X}^T\mathbf{X} + \alpha\mathbf{I})^{-1}\mathbf{X}^T\mathbf{y}\)

Lasso (\(L_1\))

Sélection de variables par parcimonie

ElasticNet

Compromis \(L_1\) / \(L_2\)

MSE, RMSE

Sensibles aux grandes erreurs

MAE

Robuste aux valeurs aberrantes

\(R^2\)

Proportion de variance expliquée

Remarque 59

La régression linéaire est rarement le modèle le plus performant sur un problème complexe, mais elle offre plusieurs avantages irremplaçables : une interprétabilité totale des coefficients, une rapidité d’entraînement, une solution analytique en forme close, et un cadre théorique bien établi (inférence, intervalles de confiance, tests). C’est le modèle de référence (baseline) avec lequel tout modèle plus complexe doit être comparé.

Contenu