--- jupytext: text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.16.0 kernelspec: name: python3 display_name: Python 3 --- # Régression linéaire multiple La régression linéaire simple est un outil puissant, mais le monde réel est rarement aussi simple qu'une seule cause. Le prix d'un logement dépend non seulement de la surface, mais aussi du quartier, du nombre de pièces, de l'état du bien. La consommation énergétique d'un bâtiment dépend des températures, de l'isolation, du nombre d'occupants. La régression multiple permet de modéliser ces relations à plusieurs dimensions, d'estimer l'effet *propre* de chaque variable (toutes choses égales par ailleurs) et de construire des modèles prédictifs robustes. ```{code-cell} python3 :tags: [hide-input] import numpy as np import pandas as pd import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec import seaborn as sns from scipy import stats import statsmodels.api as sm import statsmodels.formula.api as smf from statsmodels.stats.outliers_influence import variance_inflation_factor from sklearn.linear_model import Ridge, Lasso, ElasticNet from sklearn.preprocessing import StandardScaler from sklearn.model_selection import cross_val_score, KFold from sklearn.metrics import mean_squared_error, r2_score from sklearn.pipeline import Pipeline import pingouin as pg import warnings warnings.filterwarnings('ignore') rng = np.random.default_rng(42) plt.rcParams.update({ 'figure.dpi': 110, 'axes.spines.top': False, 'axes.spines.right': False, 'font.size': 11, }) sns.set_palette('husl') ``` ## Extension au cas multivarié Le modèle de régression linéaire multiple est : $$Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \cdots + \beta_p X_{ip} + \varepsilon_i$$ ou en notation matricielle : $Y = X\beta + \varepsilon$ où $X$ est la matrice de design $(n \times (p+1))$. Les estimateurs MCO sont : $\hat{\beta} = (X^TX)^{-1}X^TY$ **Interprétation des coefficients** : $\hat{\beta}_j$ est la variation moyenne de $Y$ pour une augmentation d'une unité de $X_j$, *toutes les autres variables maintenues constantes* — la notion de *ceteris paribus*. ```{code-cell} python3 # Données immobilières synthétiques avec plusieurs prédicteurs n = 200 surface = rng.uniform(30, 200, n) nb_pieces = np.round(surface / 30 + rng.normal(0, 0.8, n)).clip(1, 8) distance_centre = rng.uniform(1, 25, n) etage = rng.integers(0, 15, n) # Quartier (variable catégorielle) quartier = rng.choice(['Centre', 'Nord', 'Sud', 'Est'], n) quartier_coef = {'Centre': 30, 'Nord': 0, 'Sud': 10, 'Est': -5} prix = (80 + 2.8 * surface + 5.0 * nb_pieces - 1.5 * distance_centre + 0.5 * etage + np.array([quartier_coef[q] for q in quartier]) + rng.normal(0, 20, n)) df_immo = pd.DataFrame({ 'prix': prix, 'surface': surface, 'nb_pieces': nb_pieces, 'distance_centre': distance_centre, 'etage': etage, 'quartier': quartier }) # Modèle avec toutes les variables model_full = smf.ols( 'prix ~ surface + nb_pieces + distance_centre + etage + C(quartier)', data=df_immo ).fit() print(model_full.summary()) ``` ## Multicolinéarité ### Le problème Lorsque deux prédicteurs sont fortement corrélés, leur effet propre devient difficile à estimer : les estimateurs sont instables (forte variance), les intervalles de confiance s'élargissent, et les signes des coefficients peuvent s'inverser. ### VIF — Variance Inflation Factor Le **VIF** mesure le facteur par lequel la variance de $\hat{\beta}_j$ est multipliée par rapport à un modèle sans corrélation entre prédicteurs : $$\text{VIF}_j = \frac{1}{1 - R^2_j}$$ où $R^2_j$ est le $R^2$ de la régression de $X_j$ sur tous les autres prédicteurs. | VIF | Interprétation | |---|---| | 1 | Aucune multicolinéarité | | 1–5 | Modérée, généralement acceptable | | 5–10 | Forte, à surveiller | | > 10 | Sévère, problématique | ```{code-cell} python3 # Calculer les VIF X_vif = df_immo[['surface', 'nb_pieces', 'distance_centre', 'etage']].copy() X_vif_const = sm.add_constant(X_vif) vif_data = pd.DataFrame({ 'Variable': X_vif.columns, 'VIF': [variance_inflation_factor(X_vif_const.values, i + 1) for i in range(len(X_vif.columns))] }) print("VIF des prédicteurs continus :") print(vif_data.to_string(index=False)) print() print("Note : surface et nb_pieces sont corrélés → VIF > 1") print(f"Corrélation surface / nb_pieces : {np.corrcoef(surface, nb_pieces)[0,1]:.3f}") ``` ```{code-cell} python3 # Simulation : effet de la multicolinéarité sur les coefficients rng_mc = np.random.default_rng(55) n_mc = 100 # Cas 1 : prédicteurs indépendants x1_ind = rng_mc.normal(0, 1, n_mc) x2_ind = rng_mc.normal(0, 1, n_mc) y_mc = 2 * x1_ind + 3 * x2_ind + rng_mc.normal(0, 1, n_mc) model_ind = smf.ols('y ~ x1 + x2', data=pd.DataFrame({'y': y_mc, 'x1': x1_ind, 'x2': x2_ind})).fit() # Cas 2 : prédicteurs fortement corrélés (r = 0.95) x1_col = rng_mc.normal(0, 1, n_mc) x2_col = 0.95 * x1_col + np.sqrt(1 - 0.95**2) * rng_mc.normal(0, 1, n_mc) y_col = 2 * x1_col + 3 * x2_col + rng_mc.normal(0, 1, n_mc) model_col = smf.ols('y ~ x1 + x2', data=pd.DataFrame({'y': y_col, 'x1': x1_col, 'x2': x2_col})).fit() print("Cas 1 — Prédicteurs indépendants (r = 0) :") print(f" β₁ = {model_ind.params['x1']:.3f} ± {model_ind.bse['x1']:.3f} (vrai : 2)") print(f" β₂ = {model_ind.params['x2']:.3f} ± {model_ind.bse['x2']:.3f} (vrai : 3)") print() print("Cas 2 — Prédicteurs corrélés (r = 0.95) :") print(f" β₁ = {model_col.params['x1']:.3f} ± {model_col.bse['x1']:.3f} (vrai : 2)") print(f" β₂ = {model_col.params['x2']:.3f} ± {model_col.bse['x2']:.3f} (vrai : 3)") print() print(f"L'erreur-type est {model_col.bse['x1']/model_ind.bse['x1']:.1f}× plus grande avec multicolinéarité") ``` ```{code-cell} python3 :tags: [hide-input] # Visualisation : instabilité des coefficients sous multicolinéarité fig, axes = plt.subplots(1, 2, figsize=(12, 4.5)) # Distribution des estimateurs sur 500 simulations n_sim = 500 beta1_ind, beta1_col = [], [] beta2_ind, beta2_col = [], [] for _ in range(n_sim): # Indépendants x1i = rng.normal(0, 1, 80) x2i = rng.normal(0, 1, 80) yi = 2 * x1i + 3 * x2i + rng.normal(0, 1, 80) mod = smf.ols('y~x1+x2', data=pd.DataFrame({'y':yi,'x1':x1i,'x2':x2i})).fit() beta1_ind.append(mod.params['x1']) beta2_ind.append(mod.params['x2']) # Corrélés x1c = rng.normal(0, 1, 80) x2c = 0.95 * x1c + np.sqrt(1 - 0.95**2) * rng.normal(0, 1, 80) yc = 2 * x1c + 3 * x2c + rng.normal(0, 1, 80) mod2 = smf.ols('y~x1+x2', data=pd.DataFrame({'y':yc,'x1':x1c,'x2':x2c})).fit() beta1_col.append(mod2.params['x1']) beta2_col.append(mod2.params['x2']) ax = axes[0] ax.hist(beta1_ind, bins=30, density=True, alpha=0.6, color='steelblue', edgecolor='white', label='Sans multicolinéarité') ax.hist(beta1_col, bins=30, density=True, alpha=0.6, color='crimson', edgecolor='white', label='Avec multicolinéarité (r=0.95)') ax.axvline(2, color='black', lw=2.5, linestyle='--', label='Vrai β₁ = 2') ax.set_xlabel('Estimation de β₁') ax.set_ylabel('Densité') ax.set_title('Instabilité des estimateurs\nsous multicolinéarité') ax.legend(fontsize=8) ax = axes[1] ax.scatter(beta1_col, beta2_col, alpha=0.3, s=20, color='crimson', label='Avec multi. (r=0.95)') ax.scatter(beta1_ind, beta2_ind, alpha=0.3, s=20, color='steelblue', label='Sans multi.') ax.plot(2, 3, 'k*', markersize=15, label='Vraies valeurs', zorder=5) ax.set_xlabel('Estimation de β₁') ax.set_ylabel('Estimation de β₂') ax.set_title('Corrélation entre estimateurs\n(multicolinéarité = grande ellipse)') ax.legend(fontsize=8) plt.tight_layout() plt.savefig('_static/10_multicolinearite.png', bbox_inches='tight') plt.show() ``` ## Variables catégorielles : encodage dummy Les variables catégorielles sont représentées par des variables *indicatrices* (dummy) : pour $k$ catégories, on crée $k-1$ variables binaires. La catégorie omise devient la **variable de référence** — tous les coefficients s'interprètent par rapport à elle. ```{code-cell} python3 # Modèle avec variable catégorielle quartier # statsmodels / patsy gèrent automatiquement l'encodage avec C() model_cat = smf.ols('prix ~ surface + C(quartier)', data=df_immo).fit() print("Coefficients de la variable catégorielle 'quartier' :") print("(référence automatique = première catégorie alphabétique : 'Centre')") print() for param, val, se, p in zip(model_cat.params.index, model_cat.params.values, model_cat.bse.values, model_cat.pvalues.values): if 'quartier' in param or 'Intercept' in param: print(f" {param:<35} = {val:8.2f} ± {se:5.2f} (p={p:.4f})") print() print("Interprétation :") print(" Quartier Nord est la référence (β = 0 par convention)") print(" Quartier Centre : survaleur de ~30 k€ par rapport au Nord") print(" Quartier Sud : survaleur de ~10 k€ par rapport au Nord") # Encodage explicite pour illustration dummies = pd.get_dummies(df_immo['quartier'], drop_first=True) print(f"\nExemple d'encodage (drop_first=True, ref = {df_immo['quartier'].value_counts().index[-1]}) :") print(dummies.head(5)) ``` ## Interactions Un **terme d'interaction** $X_1 \times X_2$ capture l'effet de $X_1$ qui change en fonction du niveau de $X_2$. Sans interaction, l'effet de $X_1$ est le même quelle que soit la valeur de $X_2$ (hypothèse souvent trop simpliste). ```{code-cell} python3 :tags: [hide-input] # Exemple : effet de la surface sur le prix, varie selon le quartier model_inter = smf.ols( 'prix ~ surface * C(quartier)', # l'étoile inclut effets principaux + interaction data=df_immo ).fit() print("Modèle avec interaction surface × quartier :") print(f"R² = {model_inter.rsquared:.4f}") print(f"AIC = {model_inter.aic:.1f}") print() # Visualisation de l'interaction : pentes différentes par quartier fig, axes = plt.subplots(1, 2, figsize=(13, 4.5)) ax = axes[0] quartiers = ['Centre', 'Nord', 'Sud', 'Est'] colors_q = ['#4C72B0', '#DD8452', '#55A868', '#C44E52'] x_range = np.linspace(df_immo['surface'].min(), df_immo['surface'].max(), 100) for q, col in zip(quartiers, colors_q): sub = df_immo[df_immo['quartier'] == q] ax.scatter(sub['surface'], sub['prix'], alpha=0.3, color=col, s=20) # Prédiction du modèle avec interaction pour ce quartier df_pred_q = pd.DataFrame({'surface': x_range, 'quartier': q}) y_pred_q = model_inter.predict(df_pred_q) ax.plot(x_range, y_pred_q, color=col, lw=2.5, label=q) ax.set_xlabel('Surface (m²)') ax.set_ylabel('Prix (k€)') ax.set_title('Interaction surface × quartier\n(pentes différentes = interaction)') ax.legend(title='Quartier') # Comparaison AIC des modèles ax = axes[1] models = { 'Additif\n(sans inter.)': smf.ols('prix ~ surface + C(quartier)', data=df_immo).fit(), 'Interaction\nsurface × quartier': model_inter, 'Surface seule': smf.ols('prix ~ surface', data=df_immo).fit(), } aic_vals = [m.aic for m in models.values()] r2_vals = [m.rsquared for m in models.values()] noms = list(models.keys()) x_comp = np.arange(len(noms)) ax2_twin = ax.twinx() bars = ax.bar(x_comp, aic_vals, color=['#4C72B0', '#55A868', '#DD8452'], alpha=0.7) ax2_twin.plot(x_comp, r2_vals, 'k-o', lw=2, markersize=8, label='R²') ax.set_xticks(x_comp) ax.set_xticklabels(noms, fontsize=9) ax.set_ylabel('AIC (plus bas = meilleur)', color='steelblue') ax2_twin.set_ylabel('R²', color='black') ax.set_title('Comparaison des modèles\n(AIC et R²)') ax2_twin.legend() plt.tight_layout() plt.savefig('_static/10_interactions.png', bbox_inches='tight') plt.show() ``` ## Sélection de variables et régularisation ### Le problème de la sélection Plus on ajoute de prédicteurs, plus le $R^2$ augmente — mais le modèle peut *sur-ajuster* les données d'entraînement et mal se généraliser. L'objectif est de trouver le modèle le plus parcimonieux qui capture les vraies relations. ### AIC et BIC $$\text{AIC} = -2\ln(\hat{L}) + 2p \qquad \text{BIC} = -2\ln(\hat{L}) + p\ln(n)$$ Le BIC pénalise plus fortement les modèles complexes que l'AIC (facteur $\ln(n)$ vs 2). On préfère **le modèle avec l'AIC/BIC le plus bas**. ```{code-cell} python3 # Créer quelques prédicteurs supplémentaires (dont certains sans effet réel) rng_sel = np.random.default_rng(123) df_immo['bruit1'] = rng_sel.normal(0, 1, n) df_immo['bruit2'] = rng_sel.normal(0, 1, n) df_immo['bruit3'] = rng_sel.normal(0, 1, n) modeles_candidats = { 'M1: surface seulement': 'prix ~ surface', 'M2: + nb_pieces': 'prix ~ surface + nb_pieces', 'M3: + distance': 'prix ~ surface + nb_pieces + distance_centre', 'M4: + etage': 'prix ~ surface + nb_pieces + distance_centre + etage', 'M5: + quartier': 'prix ~ surface + nb_pieces + distance_centre + etage + C(quartier)', 'M6: + bruit1': 'prix ~ surface + nb_pieces + distance_centre + etage + C(quartier) + bruit1', 'M7: + bruit1,2,3': 'prix ~ surface + nb_pieces + distance_centre + etage + C(quartier) + bruit1 + bruit2 + bruit3', } comparaison = [] for nom, formule in modeles_candidats.items(): mod = smf.ols(formule, data=df_immo).fit() comparaison.append({ 'Modèle': nom, 'p': mod.df_model, 'R²': mod.rsquared, 'R² ajusté': mod.rsquared_adj, 'AIC': mod.aic, 'BIC': mod.bic, }) df_comp = pd.DataFrame(comparaison) print(df_comp.to_string(index=False)) print() print(f"Meilleur AIC : {df_comp.loc[df_comp['AIC'].idxmin(), 'Modèle']}") print(f"Meilleur BIC : {df_comp.loc[df_comp['BIC'].idxmin(), 'Modèle']}") ``` ## Régression Ridge (régularisation L2) La régression Ridge ajoute une pénalité proportionnelle au carré des coefficients : $$\min_{\beta} \left\{ \|y - X\beta\|^2 + \lambda \sum_{j=1}^p \beta_j^2 \right\}$$ La solution analytique est $\hat{\beta}_{\text{Ridge}} = (X^TX + \lambda I)^{-1}X^Ty$. Ridge **réduit** les coefficients vers zéro mais ne les annule pas — utile pour la multicolinéarité. ```{code-cell} python3 :tags: [hide-input] # Préparation des données pour sklearn feature_cols = ['surface', 'nb_pieces', 'distance_centre', 'etage', 'bruit1', 'bruit2', 'bruit3'] X_sk = df_immo[feature_cols].values y_sk = df_immo['prix'].values scaler = StandardScaler() X_scaled = scaler.fit_transform(X_sk) # Path Ridge : coefficients en fonction de lambda lambdas = np.logspace(-2, 4, 100) coefs_ridge = [] for lam in lambdas: ridge = Ridge(alpha=lam) ridge.fit(X_scaled, y_sk) coefs_ridge.append(ridge.coef_) coefs_ridge = np.array(coefs_ridge) # Path LASSO coefs_lasso = [] for lam in lambdas: try: lasso = Lasso(alpha=lam, max_iter=10000) lasso.fit(X_scaled, y_sk) coefs_lasso.append(lasso.coef_) except: coefs_lasso.append(np.zeros(len(feature_cols))) coefs_lasso = np.array(coefs_lasso) fig, axes = plt.subplots(1, 2, figsize=(14, 5)) # Path Ridge ax = axes[0] colors_feat = sns.color_palette('husl', len(feature_cols)) for j, (feat, col) in enumerate(zip(feature_cols, colors_feat)): ax.plot(np.log10(lambdas), coefs_ridge[:, j], color=col, lw=2, label=feat) ax.axvline(0, color='gray', linestyle=':', lw=1) ax.set_xlabel('log₁₀(λ)') ax.set_ylabel('Coefficient standardisé') ax.set_title('Path Ridge (L2)\n(coefficients réduits mais non nuls)') ax.legend(fontsize=8, bbox_to_anchor=(1.02, 1), loc='upper left') ax.axhline(0, color='black', lw=0.5) # Path LASSO ax = axes[1] for j, (feat, col) in enumerate(zip(feature_cols, colors_feat)): ax.plot(np.log10(lambdas), coefs_lasso[:, j], color=col, lw=2, label=feat) ax.set_xlabel('log₁₀(λ)') ax.set_ylabel('Coefficient standardisé') ax.set_title('Path LASSO (L1)\n(certains coefficients annulés = sélection)') ax.legend(fontsize=8, bbox_to_anchor=(1.02, 1), loc='upper left') ax.axhline(0, color='black', lw=0.5) plt.tight_layout() plt.savefig('_static/10_regularisation.png', bbox_inches='tight') plt.show() ``` ## Régression LASSO : sélection automatique La régression LASSO (Least Absolute Shrinkage and Selection Operator) utilise une pénalité L1 : $$\min_{\beta} \left\{ \|y - X\beta\|^2 + \lambda \sum_{j=1}^p |\beta_j| \right\}$$ La pénalité L1 produit des solutions *sparse* : certains coefficients sont exactement nuls. LASSO réalise simultanément la **régularisation** et la **sélection de variables**. ```{code-cell} python3 # Choix de lambda par validation croisée from sklearn.linear_model import LassoCV, RidgeCV lasso_cv = LassoCV(cv=5, max_iter=10000, random_state=42) lasso_cv.fit(X_scaled, y_sk) ridge_cv = RidgeCV(alphas=lambdas, cv=5) ridge_cv.fit(X_scaled, y_sk) print(f"LASSO — λ optimal (CV 5-fold) : {lasso_cv.alpha_:.4f}") print(f"\nCoefficients LASSO optimaux :") for feat, coef in zip(feature_cols, lasso_cv.coef_): status = '← annulé' if abs(coef) < 1e-10 else '' print(f" {feat:<20} = {coef:8.3f} {status}") print(f"\nRidge — λ optimal (CV 5-fold) : {ridge_cv.alpha_:.4f}") print(f"\nCoefficients Ridge optimaux :") for feat, coef in zip(feature_cols, ridge_cv.coef_): print(f" {feat:<20} = {coef:8.3f}") ``` ## Elastic Net : combiner Ridge et LASSO L'Elastic Net combine les pénalités L1 et L2 : $$\min_{\beta} \left\{ \|y - X\beta\|^2 + \lambda \left[ \rho \sum |\beta_j| + \frac{1-\rho}{2} \sum \beta_j^2 \right] \right\}$$ où $\rho \in [0,1]$ contrôle le mélange : $\rho = 1$ → LASSO, $\rho = 0$ → Ridge. ```{code-cell} python3 from sklearn.linear_model import ElasticNetCV enet_cv = ElasticNetCV(l1_ratio=[0.1, 0.5, 0.7, 0.9, 0.95, 1.0], cv=5, max_iter=10000, random_state=42) enet_cv.fit(X_scaled, y_sk) print(f"Elastic Net — λ optimal : {enet_cv.alpha_:.4f}") print(f"Elastic Net — ρ optimal (l1_ratio) : {enet_cv.l1_ratio_:.2f}") print(f"\nCoefficients Elastic Net :") for feat, coef in zip(feature_cols, enet_cv.coef_): status = '← annulé' if abs(coef) < 1e-10 else '' print(f" {feat:<20} = {coef:8.3f} {status}") ``` ## Coefficient plot avec intervalles de confiance ```{code-cell} python3 :tags: [hide-input] # Modèle statsmodels complet avec IC (sans les bruits pour clarté) model_final = smf.ols( 'prix ~ surface + nb_pieces + distance_centre + etage + C(quartier)', data=df_immo ).fit() # Extraire les coefficients et IC (hors intercept) params = model_final.params conf = model_final.conf_int() pvals = model_final.pvalues # Exclure l'intercept pour le graphe param_names = [p for p in params.index if p != 'Intercept'] coefs = params[param_names] ci_low = conf.loc[param_names, 0] ci_high = conf.loc[param_names, 1] pvals_plot = pvals[param_names] fig, axes = plt.subplots(1, 2, figsize=(14, 5)) # Coefficient plot ax = axes[0] y_pos = np.arange(len(param_names)) colors_coef = ['#55A868' if p < 0.05 else '#C44E52' for p in pvals_plot] ax.barh(y_pos, coefs, xerr=np.array([coefs - ci_low, ci_high - coefs]), color=colors_coef, alpha=0.7, height=0.6, capsize=5) ax.axvline(0, color='black', lw=1.5, linestyle='--') ax.set_yticks(y_pos) short_names = [p.replace('C(quartier)[T.', 'Q.').replace(']', '') .replace('_', ' ') for p in param_names] ax.set_yticklabels(short_names, fontsize=9) ax.set_xlabel('Coefficient ± IC 95%') ax.set_title('Coefficients standardisés\n(vert = significatif, rouge = non significatif)') from matplotlib.patches import Patch legend_elements = [Patch(facecolor='#55A868', alpha=0.7, label='p < 0.05'), Patch(facecolor='#C44E52', alpha=0.7, label='p ≥ 0.05')] ax.legend(handles=legend_elements, fontsize=9) # Validation croisée : erreur de prédiction ax = axes[1] # Comparaison MCO / Ridge / LASSO par validation croisée X_cont = df_immo[feature_cols].values X_s = StandardScaler().fit_transform(X_cont) y_v = df_immo['prix'].values kf = KFold(n_splits=5, shuffle=True, random_state=42) results_cv = {} for name, estimator in [ ('MCO (OLS)', Ridge(alpha=1e-10)), # Ridge avec alpha quasi-nul ≈ OLS ('Ridge (CV)', Ridge(alpha=ridge_cv.alpha_)), ('LASSO (CV)', Lasso(alpha=lasso_cv.alpha_, max_iter=10000)), ('Elastic Net', ElasticNet(alpha=enet_cv.alpha_, l1_ratio=enet_cv.l1_ratio_, max_iter=10000)), ]: scores = cross_val_score(estimator, X_s, y_v, cv=kf, scoring='neg_root_mean_squared_error') results_cv[name] = -scores positions = np.arange(len(results_cv)) bp = ax.boxplot(list(results_cv.values()), positions=positions, widths=0.5, patch_artist=True, medianprops=dict(color='black', lw=2)) for patch, color in zip(bp['boxes'], ['#4C72B0', '#DD8452', '#55A868', '#C44E52']): patch.set_facecolor(color) patch.set_alpha(0.7) ax.set_xticks(positions) ax.set_xticklabels(list(results_cv.keys()), fontsize=9, rotation=10) ax.set_ylabel('RMSE (validation croisée 5-fold)') ax.set_title('Comparaison des modèles\npar validation croisée') plt.tight_layout() plt.savefig('_static/10_coefficients.png', bbox_inches='tight') plt.show() ``` ## Validation croisée pour évaluer la régression La validation croisée donne une estimation honnête des performances hors-échantillon, sans nécessiter un jeu de test séparé. ```{code-cell} python3 # Validation croisée détaillée par fold from sklearn.linear_model import LinearRegression model_sk = LinearRegression() X_mod = df_immo[['surface', 'nb_pieces', 'distance_centre', 'etage']].values X_mod_s = StandardScaler().fit_transform(X_mod) kf5 = KFold(n_splits=5, shuffle=True, random_state=42) rmse_folds = [] r2_folds = [] mae_folds = [] fold_results = [] for fold, (train_idx, test_idx) in enumerate(kf5.split(X_mod_s)): X_train, X_test = X_mod_s[train_idx], X_mod_s[test_idx] y_train, y_test = y_sk[train_idx], y_sk[test_idx] model_sk.fit(X_train, y_train) y_pred = model_sk.predict(X_test) rmse = np.sqrt(mean_squared_error(y_test, y_pred)) r2 = r2_score(y_test, y_pred) mae = np.mean(np.abs(y_test - y_pred)) fold_results.append({'Fold': fold+1, 'RMSE': rmse, 'MAE': mae, 'R²': r2, 'n_train': len(train_idx), 'n_test': len(test_idx)}) df_cv = pd.DataFrame(fold_results) print("Résultats par fold (validation croisée 5-fold) :") print(df_cv.to_string(index=False)) print() print(f"RMSE moyen : {df_cv['RMSE'].mean():.2f} ± {df_cv['RMSE'].std():.2f} k€") print(f"MAE moyen : {df_cv['MAE'].mean():.2f} ± {df_cv['MAE'].std():.2f} k€") print(f"R² moyen : {df_cv['R²'].mean():.4f} ± {df_cv['R²'].std():.4f}") ``` ## Partial regression plots Les graphiques de régression partielle (*added variable plots*) montrent la relation nette entre un prédicteur et la réponse, *après avoir retiré l'effet de tous les autres prédicteurs*. ```{code-cell} python3 :tags: [hide-input] fig = plt.figure(figsize=(16, 5)) axes_pr = [] for i in range(4): axes_pr.append(fig.add_subplot(1, 4, i+1)) # Régression partielle pour chaque prédicteur continu predictors_pr = ['surface', 'nb_pieces', 'distance_centre', 'etage'] model_pr = smf.ols( 'prix ~ surface + nb_pieces + distance_centre + etage', data=df_immo ).fit() for i, pred in enumerate(predictors_pr): ax = axes_pr[i] # Résidu de y ~ autres prédicteurs autres = [p for p in predictors_pr if p != pred] formule_y = f'prix ~ {" + ".join(autres)}' formule_x = f'{pred} ~ {" + ".join(autres)}' resid_y = smf.ols(formule_y, data=df_immo).fit().resid resid_x = smf.ols(formule_x, data=df_immo).fit().resid ax.scatter(resid_x, resid_y, alpha=0.4, color='steelblue', s=30) # Droite de régression coef = np.polyfit(resid_x, resid_y, 1) x_r = np.linspace(resid_x.min(), resid_x.max(), 100) ax.plot(x_r, np.polyval(coef, x_r), 'crimson', lw=2.5) ax.axhline(0, color='gray', lw=0.5, linestyle=':') ax.axvline(0, color='gray', lw=0.5, linestyle=':') ax.set_xlabel(f'Résidu de {pred}') ax.set_ylabel('Résidu de prix' if i == 0 else '') # Corrélation partielle r_partial, p_partial = stats.pearsonr(resid_x, resid_y) ax.set_title(f'{pred}\n(r partiel = {r_partial:.2f}, p = {p_partial:.3f})') plt.suptitle('Graphiques de régression partielle\n(relation nette après contrôle des autres prédicteurs)', fontsize=12, y=1.01) plt.tight_layout() plt.savefig('_static/10_partial.png', bbox_inches='tight') plt.show() ``` ## Résumé : choisir la méthode de régression ```{code-cell} python3 :tags: [hide-input] # Tableau de synthèse recap = pd.DataFrame({ 'Méthode': ['MCO (OLS)', 'Ridge', 'LASSO', 'Elastic Net'], 'Pénalité': ['Aucune', 'L2 : λΣβ²', 'L1 : λΣ|β|', 'αρΣ|β| + α(1-ρ)Σβ²/2'], 'Sélection': ['Non', 'Non', 'Oui (sparse)', 'Oui (sparse)'], 'Multicolinéarité': ['Sensible', 'Robuste', 'Sélectionne 1 var.', 'Sélectionne groupe'], 'Hyperparamètre': ['—', 'λ (CV)', 'λ (CV)', 'λ, ρ (CV)'], 'Quand l\'utiliser': [ 'Peu de variables, hypothèses vérifiées', 'Multicolinéarité, beaucoup de variables', 'Sélection automatique souhaitée', 'Corrélation ET sélection', ] }) print(recap.to_string(index=False)) ``` ```{admonition} Checklist régression multiple :class: tip 1. **Explorer** les distributions et les corrélations entre prédicteurs (matrice de corrélation, pairplot) 2. **Calculer les VIF** pour détecter la multicolinéarité 3. **Encoder** correctement les variables catégorielles (dummy, variable de référence) 4. **Sélectionner** les variables (AIC/BIC, LASSO, théorie du domaine) 5. **Examiner les diagnostics** (résidus, QQ-plot, influence) 6. **Valider** hors-échantillon (validation croisée, R² out-of-sample) 7. **Interpréter** en termes de *ceteris paribus* et IC, pas seulement p-valeurs 8. **Tester les interactions** si des effets modificateurs sont attendus 9. **Communiquer** l'incertitude sur les prédictions (IC et IP) ```