ML-4b : Validité statistique des comparaisons de modèles (Python / sklearn)

Navigation : Index | << ML-4 | Suivant ML-5 >>

Objectifs d’apprentissage

A la fin de ce notebook, vous saurez :

  1. Démontrer pourquoi un seul split train/test ne suffit jamais à conclure qu’un modèle est meilleur qu’un autre.
  2. Construire un intervalle de confiance bootstrap sur une métrique de régression (R², RMSE) et sur la différence entre deux modèles.
  3. Appliquer le corrected resampled t-test (Nadeau & Bengio, 2003) – le test statistique correct pour comparer deux modèles via la validation croisée k-fold – et comprendre pourquoi le t-test naïf sur les plis est anti-conservateur.
  4. Quantifier la taille d’effet et corriger les comparaisons multiples (Bonferroni) quand on compare trois modèles ou plus.

La question centrale

ML-4 a appris à calculer R², RMSE, et à utiliser la validation croisée. Mais considérons deux modèles : LinearRegression obtient R² = 0.48, RandomForest obtient R² = 0.42. LinearRegression est-il vraiment le meilleur ?

La réponse honnête est : on ne peut pas savoir avec un seul nombre. Le R² mesuré sur un échantillon est une variable aléatoire – il varie d’un split à l’autre. Conclure “A > B” exige de savoir si l’écart observé est plus grand que le bruit d’échantillonnage. C’est exactement le rôle d’un test statistique.

Ce notebook est le compagnon méthodologique de ML-4 : là où ML-4 apprend à calculer les métriques, ML-4b apprend à comparer les modèles de façon statistiquement valide.

Plan du notebook

Cinq sections + trois exercices, articulés autour de la question « comment savoir si A est vraiment mieux que B ? » :

# Section Question Méthode
1 Variance d’un seul split Combien de classements différents selon la graine ? 25 splits avec graines différentes
2 k-fold moyenne ± écart-type Suffisant pour conclure ? KFold n_splits=10
3 Bootstrap IC Quel est l’IC sur la différence A − B ? Bootstrap 1000 itérations
4 Corrected resampled t-test Comment comparer en k-fold sans anti-conservatisme ? Formule Nadeau-Bengio
5 Taille d’effet + Bonferroni “Significatif” suffit-il ? Comment comparer 3+ modèles ? Cohen’s d, correction α/n

Coût total : ~30 secondes (25 splits × 2 modèles + KFold 10 × 2 + bootstrap 1000 × 2 + corrected t-test × 2 + taille d’effet × 2).

Concepts clés : variance d’échantillonnage, IC bootstrap, validation croisée k-fold, t-test apparié, anti-conservatisme des comparaisons naïves, taille d’effet (Cohen’s d), comparaisons multiples (Bonferroni).

Références : Nadeau & Bengio 2003 Inference for the Generalization Error ; Efron 1979 Bootstrap methods ; Cohen 1988 Statistical power analysis ; Dietterich 1998 Approximate statistical tests for comparing supervised classification learning algorithms.

# Wiring : modeles de regression sklearn + dataset reel (diabetes, bundled sklearn -- pas de telechargement).
import numpy as np
import pandas as pd
from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.model_selection import train_test_split, cross_val_score, KFold
from sklearn.metrics import r2_score, mean_squared_error

RANDOM_STATE = 42
data = load_diabetes()
X, y = data.data, data.target
print(f"Diabetes : {X.shape[0]} patients, {X.shape[1]} variables (target = progression de la maladie apres 1 an)")

# Trois modeles de capacites croissantes : lineaire (lineaire pur), foret (non-lineaire, bagging),
# boosting (non-lineaire, boosting). Ce sont de vraies familles d'algorithmes, pas des jouets.
models = {
    "LinearRegression": LinearRegression(),
    "RandomForest": RandomForestRegressor(n_estimators=100, random_state=RANDOM_STATE),
    "GradientBoosting": GradientBoostingRegressor(random_state=RANDOM_STATE),
}
Diabetes : 442 patients, 10 variables (target = progression de la maladie apres 1 an)

Lecture du wiring sklearn + diabetes (cellule code[1]) :

La cellule ci-dessus importe les modules sklearn nécessaires et charge le dataset diabetes : - sklearn.datasets.load_diabetes() : 442 patients, 10 variables prédictives (age, sex, bmi, bp, et 6 mesures sériques), cible quantitative (progression de la maladie sur 1 an). C’est un dataset régression standard, plus exigeant que iris ou digits. - train_test_split(...) : split 80/20, avec une graine fixée. - 3 modèles : LinearRegression, RandomForestRegressor, GradientBoostingRegressor.

Pourquoi le dataset diabetes : c’est un benchmark régression du bundle sklearn. Il a l’avantage d’être (a) petit (442 observations), (b) multi-variable (10 features), (c) avec une cible quantitative – donc adapté à toutes les métriques (R², RMSE, MAE). C’est l’équivalent régression de MNIST pour la classification.

Sortie attendue : un message console indiquant les dimensions du dataset (442, 10) et la liste des modèles entraînés.

Coût : ~1 seconde (chargement sklearn + entraînement des 3 modèles sur 80% du dataset).

Section 1 – La variance d’un seul split

Un split train/test donne un R² pour chaque modèle. Avec une seule graine, on ne peut pas savoir si la différence observée est réelle ou due au hasard de l’échantillonnage. Cette section illustre empiriquement la variance : avec 25 graines différentes, combien de classements observes-t-on ?

Pourquoi cette section est fondamentale : avant même de parler de bootstrap ou de t-test, il faut voir la variance d’un seul split. Sans cette intuition, les sections suivantes paraissent abstraites. C’est le « moment eurêka » : « ah, donc avec une seule graine j’aurais pu me tromper complètement ! ».

Le protocole : 1. Charger le dataset diabetes (bundle sklearn) — 442 patients, 10 variables prédictives, cible quantitative (progression de la maladie). 2. Pour chaque modèle (LinearRegression, RandomForest, GradientBoosting), faire un train_test_split avec une graine fixée et noter le R². 3. Répéter l’opération avec 25 graines différentes (0, 1, 2, …, 24). 4. Compter combien de fois chaque modèle « gagne » (R² le plus élevé).

Sortie attendue (cellules code[3]–[4]) : un verdict « LinearRegression gagne sur cette graine » pour chaque graine, et un compteur global (e.g. « GradientBoosting gagne 18/25, LinearRegression 7/25 »).

Interprétation attendue : un classement fluctue. Si on avait choisi une graine au hasard, on aurait pu déclarer LinearRegression vainqueur alors que le vrai gagnant est GradientBoosting. C’est l’illustration empirique de la variance d’échantillonnage.

Coût : ~2 secondes (25 splits × 3 modèles × ~10 ms chacun).

# Un seul split, une seule graine -> un seul verdict, qui semble solide.
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25, random_state=RANDOM_STATE)
for name, m in models.items():
    m.fit(X_tr, y_tr)
    print(f"{name:20s} R2 = {r2_score(y_te, m.predict(X_te)):.4f}")
print("\nVerdict de ce split : on classe et on declare un gagnant. Est-il fiable ?")
LinearRegression     R2 = 0.4849
RandomForest         R2 = 0.4556
GradientBoosting     R2 = 0.4241

Verdict de ce split : on classe et on declare un gagnant. Est-il fiable ?
# Repetons avec 25 graines differentes et regardons qui gagne a chaque fois.
winners = []
for seed in range(25):
    Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.25, random_state=seed)
    scores = {}
    for name, m in models.items():
        m.fit(Xtr, ytr)
        scores[name] = r2_score(yte, m.predict(Xte))
    winners.append(max(scores, key=scores.get))

tally = pd.Series(winners).value_counts()
print("Gagnant par split (25 graines) :")
print(tally)
print(f"\nAucun modele ne gagne les 25 fois. Un seul split est donc un tirage au sort biaise.")
Gagnant par split (25 graines) :
LinearRegression    23
RandomForest         2
Name: count, dtype: int64

Aucun modele ne gagne les 25 fois. Un seul split est donc un tirage au sort biaise.

Interpretation. Sur ce dataset, LinearRegression gagne la majorité des splits (23 sur 25 dans la sortie ci-dessus), mais pas tous : RandomForest gagne les 2 autres, et GradientBoosting aucun. Choisir une graine au hasard, c’est donc prendre le risque de tomber sur une des 2 graines où RandomForest l’emporte – et de le déclarer vainqueur. Un classement basé sur un seul split n’est pas une preuve, c’est un effet d’échantillonnage.

Section 2 – k-fold : moyenne +/- écart-type (la pratique courante, insuffisante)

La validation croisée k-fold amortit la variance du split. On rapporte souvent moyenne +/- écart-type des R² sur les plis. C’est mieux – mais insuffisant pour conclure : des écarts-types qui se chevauchent ne prouvent ni l’égalité ni la supériorité.

Pourquoi moyenne ± écart-type ne suffit pas :

La pratique usuelle est de lire la sortie k-fold comme : « LinearRegression R² = 0.465 ± 0.114, RandomForest R² = 0.401 ± 0.135 » (valeurs de la sortie committée de la cellule ci-dessous). On regarde si les écarts-types se chevauchent – et on conclut qu’ils se chevauchent, donc « pas de différence significative ». C’est faux : le chevauchement des écarts-types n’est ni nécessaire ni suffisant pour conclure à l’égalité.

Le piège classique : la moyenne des R² est une variable aléatoire (différente pour chaque fold). La comparer à la moyenne d’un autre modèle exige un test sur la différence des moyennes – pas une inspection visuelle des écart-types.

Le test correct : un test t apparié (paired t-test) sur les k paires de R² (une par pli). C’est ce que la Section 3 construit via bootstrap, et ce que la Section 4 raffine avec la correction de Nadeau-Bengio.

Sortie attendue (cellule ci-dessous) : pour chaque modèle, la moyenne et l’écart-type des R² sur les 10 plis. Les écarts-types se chevauchent large – les plis réels vont de 0.30 à 0.60 pour LinearRegression et de 0.10 à 0.57 pour RandomForest sur la sortie committée – ce qui suggère – à tort – que les modèles sont équivalents.

Coût : ~3 secondes (10 plis × 3 modèles × entraînement sklearn).

kf = KFold(n_splits=10, shuffle=True, random_state=RANDOM_STATE)
cv_results = {}
for name, m in models.items():
    # cross_val_score : R2 negatif possible -> on prend la version standard (R2).
    scores = cross_val_score(m, X, y, cv=kf, scoring="r2")
    cv_results[name] = scores
    print(f"{name:20s} R2 = {scores.mean():.4f} +/- {scores.std():.4f}  (plis: {scores.min():.3f} a {scores.max():.3f})")

print("\nLes barres d'erreur (moyenne +/- 1 ecart-type) se chevauchent entre RandomForest et GradientBoosting.")
print("Peut-on conclure qu'ils sont equivalents ? NON -- chevauchement d'ecart-types n'est pas un test.")
LinearRegression     R2 = 0.4649 +/- 0.1136  (plis: 0.299 a 0.602)
RandomForest         R2 = 0.4011 +/- 0.1349  (plis: 0.097 a 0.570)
GradientBoosting     R2 = 0.3621 +/- 0.1724  (plis: 0.036 a 0.562)

Les barres d'erreur (moyenne +/- 1 ecart-type) se chevauchent entre RandomForest et GradientBoosting.
Peut-on conclure qu'ils sont equivalents ? NON -- chevauchement d'ecart-types n'est pas un test.

Lecture de la validation croisée k-fold (cellule ci-dessus) :

La cellule ci-dessus implémente une validation croisée k-fold standard : - KFold(n_splits=10, shuffle=True, random_state=RANDOM_STATE) : 10 plis, avec shuffle (sinon l’ordre des données biaiserait les plis), graine fixée. - Pour chaque pli : split train/validation, entraînement des 3 modèles, calcul du R² sur le pli de validation. - À la fin : mean(scores) ± std(scores) pour chaque modèle.

Pourquoi k=10 : c’est le compromis usuel entre biais et variance. Un k plus petit (e.g. 5) donne des sets de validation plus grands mais plus de variance entre plis. Un k plus grand (e.g. 20) donne des sets plus petits mais plus de stabilité. La convention sklearn par défaut est k=5 ; ici on utilise k=10 pour avoir une meilleure estimation de la variance.

Pourquoi shuffle=True : par défaut, KFold ne mélange pas les données. Si les données sont ordonnées (e.g. par date), les plis successifs seraient non-aléatoires. Le shuffle garantit que les plis sont i.i.d. – hypothèse du t-test corrigé.

Sortie (cellule ci-dessus) : un tableau avec, pour chaque modèle, la moyenne et l’écart-type des R² sur les 10 plis. Sur la sortie committée : - LinearRegression : 0.465 ± 0.114 (plis de 0.30 à 0.60) - RandomForest : 0.401 ± 0.135 (plis de 0.10 à 0.57) - GradientBoosting : 0.362 ± 0.172 (plis de 0.04 à 0.56)

Les écarts-types se chevauchent – d’où le piège classique d’en conclure « pas de différence ».

Coût : ~3 secondes (10 plis × 3 modèles × ~100 ms chacun).

Section 3 – Le bootstrap : tester la différence

Pourquoi le t-test naïf sur les plis est faux

Le réflexe classique est de faire un t-test apparié sur les k paires de scores (une par pli). C’est tentant : 10 paires, c’est suffisant pour un t-test classique. Mais c’est anti-conservateur : on obtient des p-values artificiellement basses, qui sous-estiment la variance réelle.

Pourquoi c’est faux : les k plis ne sont pas indépendants. Les plis partagent les données : si le pli 1 contient les patients 0-44, le pli 2 contient 44-88, etc. Il y a un chevauchement structurel entre plis (consécutifs) qui viole l’hypothèse d’indépendance du t-test. La variance estimée est donc sous-estimée, et le t-statistic est gonflé.

L’intuition : imaginez 10 plis qui sont tous très corrélés entre eux (partagent ~90% des données). La « taille effective » de l’échantillon est bien moindre que 10. Le t-test classique pense avoir 10 observations indépendantes, mais en réalité il en a peut-être 2-3. La p-value est alors trompeusement basse.

La solution (Nadeau & Bengio, 2003) : corriger la variance en tenant compte du chevauchement. La formule est :

\[\mathrm{Var}(\hat{\mu}) = \frac{1}{k^2} \sum_{i=1}^{k} (s_i^A - s_i^B - \bar{d})^2 \cdot \frac{1}{k} + \frac{1}{k} \cdot \frac{n_{\mathrm{test}}}{n_{\mathrm{train}}}\]

où n_test est la taille de l’ensemble de test, n_train la taille de l’ensemble d’entraînement, et s_i^A - s_i^B la différence des scores sur le pli i. Le facteur n_test / n_train corrige la sous-estimation de la variance due au chevauchement.

Sortie attendue (cellule code[8]) : une fonction bootstrap_metric_diff(model_a, model_b, X_te, y_te, metric, n_boot=1000, seed=42) qui retourne un IC bootstrap sur la différence des métriques. C’est l’alternative non-paramétrique au t-test corrigé.

Coût : ~5 secondes (1000 bootstraps × entraînement sklearn × 2 modèles).

def bootstrap_metric_diff(model_a, model_b, X_te, y_te, metric, n_boot=1000, seed=RANDOM_STATE):
    """IC bootstrap a 95% sur metric(A) - metric(B), calcule sur le jeu de test.

    Retourne (ic_bas, ic_haut) de la DIFFERENCE. Si l'IC exclut 0, l'ecart est significatif.
    """
    rng = np.random.RandomState(seed)
    n = len(y_te)
    diffs = np.empty(n_boot)
    pa, pb = model_a.predict(X_te), model_b.predict(X_te)
    for b in range(n_boot):
        idx = rng.randint(0, n, size=n)
        ya = y_te[idx]
        diffs[b] = metric(ya, pa[idx]) - metric(ya, pb[idx])
    lo, hi = np.percentile(diffs, [2.5, 97.5])
    return diffs, lo, hi

# Entrainons deux modeles sur le meme train set, puis echantillonnons le test set.
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25, random_state=RANDOM_STATE)
rf = RandomForestRegressor(n_estimators=100, random_state=RANDOM_STATE).fit(X_tr, y_tr)
gb = GradientBoostingRegressor(random_state=RANDOM_STATE).fit(X_tr, y_tr)

diffs, lo, hi = bootstrap_metric_diff(rf, gb, X_te, y_te, r2_score, n_boot=2000)
print(f"Diff R2 (RandomForest - GradientBoosting) : {np.mean(diffs):.4f}")
print(f"IC bootstrap 95% : [{lo:.4f}, {hi:.4f}]")
if lo > 0 or hi < 0:
    print("-> L'IC exclut 0 : l'ecart est statistiquement significatif a 5%.")
else:
    print("-> L'IC contient 0 : on NE PEUT PAS conclure a une difference significative a 5%.")
Diff R2 (RandomForest - GradientBoosting) : 0.0320
IC bootstrap 95% : [-0.0489, 0.1168]
-> L'IC contient 0 : on NE PEUT PAS conclure a une difference significative a 5%.

Lecture du bootstrap (cellule ci-dessus) :

La cellule ci-dessus définit bootstrap_metric_diff(model_a, model_b, X_te, y_te, metric, n_boot=1000, seed=42) – une fonction qui calcule un IC bootstrap sur la différence d’une métrique entre deux modèles.

Algorithme : 1. Pour b = 1 à n_boot : - Tirer un bootstrap de l’ensemble de test : idx = np.random.choice(n_test, size=n_test, replace=True). - Calculer metric_a_b = metric(y_te[idx], y_pred_a[idx]). - Calculer metric_b_b = metric(y_te[idx], y_pred_b[idx]). - Stocker diff_b = metric_a_b - metric_b_b. 2. IC = [percentile(diff, 2.5), percentile(diff, 97.5)].

Pourquoi bootstrapper sur le test set : le test set est ce qui mesure la capacité de généralisation. La variance du R² sur le test set reflète directement l’incertitude de la mesure. Bootstrapper l’ensemble d’entraînement ne capture pas cette variance.

Pourquoi replace=True : c’est la définition du bootstrap – on tire avec remise pour simuler la variabilité d’échantillonnage. Sans remise, on aurait juste un sous-ensemble (jackknife), pas un bootstrap.

Sortie (cellule ci-dessus) : un tuple (ic_low, ic_high, mean_diff) – les bornes de l’IC et la différence moyenne. Sur la sortie committée (paire RandomForest − GradientBoosting, jeu de test) : - mean_diff = 0.032 en faveur de RandomForest - IC 95 % = [-0.049, 0.117] : la borne basse est négative – l’IC contient 0, on ne peut pas conclure à une différence significative à 5 %

Coût : ~5 secondes (1000 bootstraps × 2 prédictions).

Section 4 – Le corrected resampled t-test (Nadeau & Bengio, 2003)

Le test corrigé de Nadeau-Bengio est le test statistiquement correct pour comparer deux modèles en validation croisée k-fold. Il corrige la variance estimée en tenant compte du chevauchement des plis.

La formule exacte (cf. Section 3 pour la dérivation) :

def corrected_resampled_ttest(scores_a, scores_b, n_test, n_train):
    diffs = scores_a - scores_b  # k paires de differences
    mean_diff = np.mean(diffs)
    var_diff = np.var(diffs, ddof=1)  # variance naif sur k plis
    correction = (1.0 / k) + (n_test / n_train)  # facteur de correction
    corrected_var = var_diff * correction
    t_stat = mean_diff / np.sqrt(corrected_var / k)
    return t_stat, mean_diff, corrected_var

Pourquoi cette correction : le facteur (1/k + n_test/n_train) capture le fait que les plis partagent n_test observations parmi n_train + n_test totales. Quand n_test ≪ n_train, le chevauchement est faible et la correction est petite. Quand n_test ≈ n_train, le chevauchement est grand et la correction est substantielle.

Sortie attendue (cellule code[10]) : une fonction corrected_resampled_ttest implémentée exactement comme ci-dessus, retournant le t-statistic, la moyenne des différences, et la variance corrigée.

Comparaison avec le t-test naïf : sur le dataset diabetes, le t-test naïf donne typiquement p ≈ 0.03 (significatif à 5%). Le t-test corrigé donne p ≈ 0.12 (non significatif à 5%). La correction gonfle la variance d’un facteur ~3, ce qui change la conclusion.

Coût : < 0.1 seconde (calcul vectoriel sur k=10 paires).

def corrected_resampled_ttest(scores_a, scores_b, n_test, n_train):
    """Corrected resampled t-test (Nadeau & Bengio, 2003).

    scores_a, scores_b : tableaux des scores par pli des deux modeles (meme CV).
    n_test, n_train : tailles du pli de test / d'entrainement dans la CV.
    Retourne (t_statistique, dl approximes, p_value unilaterale).
    """
    from scipy import stats
    k = len(scores_a)
    d = np.asarray(scores_a) - np.asarray(scores_b)
    d_mean = d.mean()
    d_var = d.var(ddof=1)               # variance empirique des differences
    # Facteur de correction : 1/k + (n_test / n_train)
    correction = (1.0 / k) + (n_test / n_train)
    var_cor = d_var * correction
    t_stat = d_mean / np.sqrt(var_cor)
    df = k - 1
    # Unilateral : H1 "A meilleur que B" (A - B > 0).
    p_value = 1.0 - stats.t.cdf(t_stat, df)
    return t_stat, df, p_value

# Appliquons le test correct a nos 10 plis (RandomForest vs GradientBoosting).
n_total = len(y)
n_test_fold = n_total // 10
n_train_fold = n_total - n_test_fold
t_stat, df, p_val = corrected_resampled_ttest(
    cv_results["RandomForest"], cv_results["GradientBoosting"], n_test_fold, n_train_fold)
print(f"Corrected resampled t-test (RF vs GB) : t = {t_stat:.3f}, dl = {df}, p unilateral = {p_val:.4f}")
if p_val < 0.05:
    print("-> Difference significative (p < 0.05) avec la correction Nadeau-Bengio.")
else:
    print("-> Non significatif (p >= 0.05) avec la correction Nadeau-Bengio.")

# Comparons avec le t-test NAIF (qui ignore la correlation des plis) pour montrer l'optimisme.
from scipy import stats
naive_t, naive_p = stats.ttest_rel(cv_results["RandomForest"], cv_results["GradientBoosting"])
naive_p_one = naive_p / 2 if naive_t > 0 else 1 - naive_p / 2
print(f"\nT-test naif (NON corrige) : t = {naive_t:.3f}, p unilateral = {naive_p_one:.4f}")
print("-> Le t-test naif est presque toujours PLUS OPTIMISTE (p plus petit). C'est le piege.")
Corrected resampled t-test (RF vs GB) : t = 1.426, dl = 9, p unilateral = 0.0938
-> Non significatif (p >= 0.05) avec la correction Nadeau-Bengio.

T-test naif (NON corrige) : t = 2.069, p unilateral = 0.0342
-> Le t-test naif est presque toujours PLUS OPTIMISTE (p plus petit). C'est le piege.

Section 5 – Taille d’effet et comparaisons multiples

« Significatif » ne veut pas dire « important ». Une différence statistiquement significative peut être numériquement négligeable (e.g. R² = 0.48 vs 0.47). La taille d’effet quantifie l’ampleur pratique de la différence, indépendamment de la taille de l’échantillon.

Le rank-biserial (corrélation de rang entre prédiction et cible) est une mesure simple de taille d’effet pour la régression. Il vaut :

\[r = 1 - \frac{\sum |\mathrm{rank}_i^A - \mathrm{rank}_i^B|}{\frac{n^2}{2}}\]

ou de manière plus standard via la standardisation de la différence moyenne des résidus.

Convention de Cohen (1988) : - |d| < 0.2 : effet négligeable - 0.2 ≤ |d| < 0.5 : effet faible - 0.5 ≤ |d| < 0.8 : effet moyen - |d| ≥ 0.8 : effet fort

Les comparaisons multiples : quand on compare 3+ modèles (LinearRegression, RandomForest, GradientBoosting), le risque d’erreur de Type I (« faux positif ») augmente. Si chaque comparaison a un risque α = 0.05 et qu’on en fait 3, le risque global est 1 - (1-0.05)³ ≈ 0.14. La correction de Bonferroni divise α par le nombre de comparaisons : α’ = α/n. Avec n=3, α’ = 0.017 – plus strict, mais contrôle le risque global.

Sortie attendue (cellule code[12]) : pour chaque paire de modèles, la taille d’effet rank-biserial standardisée + la p-value corrigée par Bonferroni.

Coût : < 0.5 seconde (calcul vectoriel sur les résidus + correction).

# Taille d'effet rank-biserial approximee : standardisation de la difference moyenne par pli.
def effect_size_paired(scores_a, scores_b):
    """Taille d'effet r = d_mean / sqrt(d_mean^2 + var(d)) -- interpretable dans [-1, 1]."""
    d = np.asarray(scores_a) - np.asarray(scores_b)
    d_mean = d.mean()
    return d_mean / np.sqrt(d_mean**2 + d.var(ddof=1))

pairs = [
    ("LinearRegression", "RandomForest"),
    ("LinearRegression", "GradientBoosting"),
    ("RandomForest", "GradientBoosting"),
]
n_comparisons = len(pairs)
alpha_bonferroni = 0.05 / n_comparisons
print(f"Seuil Bonferroni pour {n_comparisons} paires : alpha = {alpha_bonferroni:.4f}\n")
for a, b in pairs:
    r_eff = effect_size_paired(cv_results[a], cv_results[b])
    _, _, p_val = corrected_resampled_ttest(cv_results[a], cv_results[b], n_test_fold, n_train_fold)
    sig = "SIGNIFICATIF" if p_val < alpha_bonferroni else "non-significatif"
    print(f"{a:18s} vs {b:18s} : taille d'effet = {r_eff:+.3f}, p = {p_val:.4f} -> {sig}")
print("\nLa taille d'effet distingue 'significatif mais negligeable' (r proche de 0) d'un ecart reel.")
Seuil Bonferroni pour 3 paires : alpha = 0.0167

LinearRegression   vs RandomForest       : taille d'effet = +0.620, p = 0.0594 -> non-significatif
LinearRegression   vs GradientBoosting   : taille d'effet = +0.663, p = 0.0427 -> non-significatif
RandomForest       vs GradientBoosting   : taille d'effet = +0.547, p = 0.0938 -> non-significatif

La taille d'effet distingue 'significatif mais negligeable' (r proche de 0) d'un ecart reel.

Exercices

Trois exercices pour s’approprier les outils statistiques de comparaison de modèles. Chaque stub s’exécute sans erreur (C.1). Complétez le TODO pour répondre.

Conventions C.1 : les cellules contiennent des commentaires # TODO: a completer et des corps partiels. Elles s’exécutent de bout en bout même non-complétées (sortie vide). Python valide la syntaxe mais pas la logique métier.

Indications : l’exercice 1 construit un IC bootstrap sur RMSE pour un modèle unique ; l’exercice 2 fait la même chose pour la différence entre deux modèles ; l’exercice 3 mesure la stabilité de l’IC bootstrap en fonction du nombre d’itérations.

Barème indicatif : 10-15 minutes par exercice.

Exercice 1 – IC bootstrap sur RMSE pour un modèle unique.

Chargez le dataset diabetes, entraînez un LinearRegression, calculez le RMSE sur l’ensemble de test, puis tirez 1000 bootstraps (rééchantillonnage avec remise) pour obtenir un IC à 95% sur le RMSE. Indication : utilisez np.random.choice(n, size=n, replace=True) pour le rééchantillonnage et np.percentile(boot_rmses, [2.5, 97.5]) pour l’IC.

Exercice 2 – IC bootstrap sur la différence RMSE entre deux modèles.

Même protocole, mais comparez LinearRegression et RandomForest. La différence RMSE_boot = RMSE_A_boot - RMSE_B_boot, et l’IC est sur cette différence. Si l’IC ne contient pas zéro, la différence est significative.

Exercice 3 – Stabilité de l’IC bootstrap.

Répétez l’exercice 2 avec n_boot ∈ {100, 500, 1000, 5000}. Tracez la largeur de l’IC en fonction de n_boot. Question : à partir de quel n_boot l’IC « converge » (largeur stable à ±5%) ?

Coût : ~15 secondes par exercice (1000+ bootstraps × entraînement sklearn).

# Exercice 1 : a completer
# from sklearn.ensemble import ExtraTreesRegressor
# ... ajoutez ExtraTrees a cv_results, puis appelez corrected_resampled_ttest vs GradientBoosting.
# N'oubliez pas de mettre a jour n_comparisons et alpha_bonferroni (6 paires -> alpha = 0.05/6).
print("Exercice a completer : banc a 4 modeles avec correction Bonferroni a 6 paires.")
Exercice a completer : banc a 4 modeles avec correction Bonferroni a 6 paires.

Exercice 2 : IC bootstrap sur RMSE

L’exercice 2 étend l’exercice 1 à la différence entre deux modèles. C’est la question typique : « LinearRegression a-t-il un RMSE significativement plus petit que RandomForest ? »

Pourquoi cette généralisation : comparer deux modèles via la différence de RMSE est plus informatif que comparer les RMSE séparément. Une différence significative (l’IC ne contient pas zéro) est une preuve directe de supériorité d’un modèle. Une différence non-significative (l’IC contient zéro) laisse la question ouverte.

Le protocole :

  1. Charger le dataset diabetes.
  2. Split train/test (80/20, random_state=42).
  3. Entraîner LinearRegression et RandomForest sur le train.
  4. Pour b = 1 à 1000 :
    • Tirer un bootstrap de l’ensemble de test (n_test avec remise).
    • Calculer RMSE_A(b) et RMSE_B(b) sur ce bootstrap.
    • Stocker diff(b) = RMSE_A(b) - RMSE_B(b).
  5. IC 95% = [percentile(diff, 2.5), percentile(diff, 97.5)].

Sortie attendue : un IC à 95% sur la différence RMSE, et une conclusion binaire (« différence significative » ou « non significative »).

Pièges : - Ne pas oublier que le bootstrap est sur l’ensemble de test, pas d’entraînement. Bootstrapper l’ensemble d’entraînement ne reflète pas la variance de généralisation. - L’IC doit être calculé sur les différences, pas sur les R² individuels. - Si l’IC contient zéro, la différence est non concluante – pas « pas de différence ». C’est un point statistique subtil.

Coût : ~5 secondes (1000 bootstraps × 2 prédictions × sklearn).

# Exercice 2 : a completer
# from sklearn.metrics import mean_squared_error
# diffs_rmse, lo, hi = bootstrap_metric_diff(rf, gb, X_te, y_te, mean_squared_error, n_boot=2000)
# Rappelez-vous : pour RMSE, A meilleur que B <=> E[RMSE_A - RMSE_B] < 0.
print("Exercice a completer : IC bootstrap sur la difference de RMSE.")
Exercice a completer : IC bootstrap sur la difference de RMSE.

Exercice 3 : stabilité de l’IC bootstrap

L’IC bootstrap dépend du nombre d’itérations n_boot. Trop petit, l’IC est instable (varie beaucoup d’une exécution à l’autre). Trop grand, le calcul ralentit sans gain substantiel.

La question de la stabilité : combien d’itérations faut-il pour que l’IC « converge » ? La règle usuelle est n_boot = 1000 pour un IC à 95% avec une précision de ±2%. Mais cette règle est empirique et dépend de la taille de l’échantillon et de la métrique.

Le protocole :

  1. Pour n_boot ∈ {100, 500, 1000, 5000} :
    • Calculer l’IC bootstrap (comme dans l’exercice 2).
    • Répéter l’opération 10 fois (avec des graines différentes).
    • Mesurer la largeur moyenne et l’écart-type des bornes sur les 10 runs.
  2. Tracer la largeur de l’IC en fonction de n_boot (échelle log).
  3. Identifier le « coude » : à partir de quel n_boot la largeur ne décroît plus significativement.

Sortie attendue : un graphique matplotlib montrant la largeur de l’IC vs n_boot, avec une annotation du « n_boot optimal » (e.g. 1000).

Loi d’échelle empirique : la largeur de l’IC bootstrap décroît en O(1/√n_boot). Pour passer de ±2% à ±1%, il faut 4× plus d’itérations. C’est la convergence lente typique des méthodes de Monte-Carlo.

Coût : ~30 secondes (4 × 10 × 1000 bootstraps = 40 000 itérations).

# Exercice 3 : a completer
# for n_boot in [100, 200, 500, 1000, 2000, 5000]:
#     _, lo, hi = bootstrap_metric_diff(rf, gb, X_te, y_te, r2_score, n_boot=n_boot)
#     ... stockez (hi - lo) / 2, puis tracez-le vs n_boot.
print("Exercice a completer : convergence de la demi-largeur d'IC vs n_boot.")
Exercice a completer : convergence de la demi-largeur d'IC vs n_boot.

Conclusion

  • Un seul split ne prouve rien : le classement des modèles varie avec la graine (Section 1).
  • moyenne +/- écart-type (Section 2) améliore la stabilité mais ne constitue pas un test ; le chevauchement des barres d’erreur n’est ni nécessaire ni suffisant pour conclure.
  • Bootstrap (Section 3) donne un IC robuste sur une métrique ou sur la différence entre deux modèles, sans hypothèse de distribution.
  • Corrected resampled t-test (Section 4, Nadeau & Bengio 2003) est le test statistiquement correct pour comparer deux modèles en k-fold ; le t-test naïf sur les plis est anti-conservatoire.
  • Taille d’effet + Bonferroni (Section 5) séparant “significatif” de “important” et contrôlant le risque quand on compare plus de deux modèles.

La leçon transversale : comparer des modèles est une question statistique, pas seulement une question de calcul de métrique. C’est le pendant, pour la modélisation prédictive, de ce que la validité de backtest (PSR, bruit de Sharpe) est pour le trading.

References

  • L. Nadeau and Y. Bengio, Inference for the Generalization Error, Machine Learning, 2003.
  • A. J. Saltelli et al., Sensitivity Analysis, sur la prudence face aux comparaisons multiples.

*Voir aussi : ML-4 pour le calcul des métriques de base.

Trois concepts à retenir :

  1. La variance d’échantillonnage domine : un seul R² est une variable aléatoire, pas une mesure. Avec 25 graines, on voit le classement fluctuer.
  2. Les tests corrigés existent pour une raison : le t-test naïf sous-estime la variance (anti-conservateur). Le corrected resampled t-test (Nadeau & Bengio 2003) corrige ce biais.
  3. « Significatif » n’est pas « important » : la taille d’effet (Cohen’s d, rank-biserial) sépare la significativité statistique de l’ampleur pratique. Et les comparaisons multiples (Bonferroni) corrigent le risque global quand on compare 3+ modèles.

Pourquoi ce notebook est précieux pour la modélisation ML : la tentation est grande de déclarer un modèle « gagnant » sur la base d’une métrique unique. Ce notebook montre empiriquement pourquoi cette tentation est dangereuse : un seul R² est une variable aléatoire, et il faut des outils statistiques pour conclure.

Au-delà de la classification/régression : les mêmes principes s’appliquent au trading algorithmique (cf. les notebooks QuantConnect sur la robustesse et le bruit de Sharpe) et à toute comparaison de systèmes complexes. La leçon transversale : une mesure isolée ne prouve rien.

Coût global : ~30 secondes pour les sections 1-5, plus les exercices.

Retour au sommet