Régression linéaire simple#

La corrélation décrit une relation ; la régression la modélise. Régresse-t-on la surface habitable sur le prix d’un logement, la dose d’un médicament sur la réponse biologique, ou les heures de travail sur la performance ? Dans chaque cas, on cherche non seulement à quantifier la relation, mais à prédire, à expliquer, et à quantifier l’incertitude de ces prédictions. Ce chapitre présente la régression linéaire simple — une variable explicative — en insistant sur les diagnostics et l’interprétation.

Hide code cell source

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.diagnostic import het_breuschpagan, het_white
from statsmodels.graphics.gofplots import ProbPlot
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')

Motivation : prédire et expliquer#

La régression linéaire simple modélise la relation entre une variable explicative (prédicteur) \(X\) et une variable réponse \(Y\) :

\[Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i \qquad \varepsilon_i \overset{i.i.d.}{\sim} \mathcal{N}(0, \sigma^2)\]
  • \(\beta_0\) : ordonnée à l’origine (valeur de \(Y\) quand \(X = 0\))

  • \(\beta_1\) : pente — changement moyen de \(Y\) pour une augmentation d’une unité de \(X\)

  • \(\varepsilon_i\) : terme d’erreur, non observable, capturant tout ce que le modèle ne prédit pas

Hide code cell source

# Exemple : surface habitable (m²) → prix de vente (k€)
n = 80
surface = rng.uniform(30, 150, n)
# Prix = 50 + 3.5 * surface + bruit gaussien
beta_0_vrai, beta_1_vrai = 50, 3.5
sigma_vrai = 25
prix = beta_0_vrai + beta_1_vrai * surface + rng.normal(0, sigma_vrai, n)

df_immo = pd.DataFrame({'surface': surface, 'prix': prix})

# Nuage de points exploratoire
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.scatter(surface, prix, alpha=0.6, color='steelblue', s=50)
ax.set_xlabel('Surface (m²)')
ax.set_ylabel('Prix (k€)')
ax.set_title('Prix de vente en fonction de la surface\n(données simulées, marché local)')
plt.tight_layout()
plt.savefig('_static/09_intro.png', bbox_inches='tight')
plt.show()
_images/4419b884568bf2eb70d433434dad150bcd08e929b46c7c59ee912029677d562d.png

Moindres Carrés Ordinaires (MCO)#

Intuition géométrique#

Les estimateurs MCO minimisent la somme des carrés des résidus : $\(\min_{\beta_0, \beta_1} \sum_{i=1}^n (y_i - \beta_0 - \beta_1 x_i)^2 = \min_{\beta_0, \beta_1} \|y - X\beta\|^2\)$

Les résidus \(\hat{\varepsilon}_i = y_i - \hat{y}_i\) sont les distances verticales entre les points et la droite ajustée.

Formules et propriétés BLUE#

Les estimateurs MCO sont : $\(\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{S_{XY}}{S_{XX}} \qquad \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}\)$

Le théorème de Gauss-Markov établit que sous les hypothèses classiques (linéarité, homoscédasticité, erreurs non corrélées, \(X\) non stochastique ou exogène), les estimateurs MCO sont BLUE (Best Linear Unbiased Estimators) : ils ont la variance minimale parmi tous les estimateurs linéaires sans biais.

# Ajustement MCO avec statsmodels
model = smf.ols('prix ~ surface', data=df_immo).fit()
print(model.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                   prix   R-squared:                       0.962
Model:                            OLS   Adj. R-squared:                  0.962
Method:                 Least Squares   F-statistic:                     1987.
Date:                Wed, 01 Apr 2026   Prob (F-statistic):           2.99e-57
Time:                        22:06:21   Log-Likelihood:                -364.72
No. Observations:                  80   AIC:                             733.4
Df Residuals:                      78   BIC:                             738.2
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     38.1162      7.780      4.899      0.000      22.628      53.604
surface        3.5970      0.081     44.571      0.000       3.436       3.758
==============================================================================
Omnibus:                        3.950   Durbin-Watson:                   1.522
Prob(Omnibus):                  0.139   Jarque-Bera (JB):                3.155
Skew:                           0.411   Prob(JB):                        0.207
Kurtosis:                       3.520   Cond. No.                         287.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
# Extraire les coefficients et leurs intervalles de confiance
beta0_est = model.params['Intercept']
beta1_est = model.params['surface']
ic_95 = model.conf_int(alpha=0.05)

print(f"\nEstimations MCO :")
print(f"  β₀ = {beta0_est:.2f}  (vrai : {beta_0_vrai})")
print(f"  β₁ = {beta1_est:.2f}  (vrai : {beta_1_vrai})")
print(f"  σ̂  = {np.sqrt(model.mse_resid):.2f}  (vrai : {sigma_vrai})")
print(f"\nIntervalles de confiance à 95% :")
print(f"  β₀ ∈ [{ic_95.loc['Intercept', 0]:.2f}, {ic_95.loc['Intercept', 1]:.2f}]")
print(f"  β₁ ∈ [{ic_95.loc['surface', 0]:.2f}, {ic_95.loc['surface', 1]:.2f}]")
Estimations MCO :
  β₀ = 38.12  (vrai : 50)
  β₁ = 3.60  (vrai : 3.5)
  σ̂  = 23.40  (vrai : 25)

Intervalles de confiance à 95% :
  β₀ ∈ [22.63, 53.60]
  β₁ ∈ [3.44, 3.76]

Coefficient de détermination R²#

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

\[R^2 = 1 - \frac{SS_{\text{résidus}}}{SS_{\text{total}}} = 1 - \frac{\sum (y_i - \hat{y}_i)^2}{\sum (y_i - \bar{y})^2}\]
ss_total = np.sum((prix - prix.mean())**2)
ss_res = np.sum(model.resid**2)
ss_reg = ss_total - ss_res

r2 = model.rsquared
r2_adj = model.rsquared_adj

print(f"Décomposition de la variance :")
print(f"  SS total     = {ss_total:.1f}")
print(f"  SS régression= {ss_reg:.1f}")
print(f"  SS résidus   = {ss_res:.1f}")
print(f"\nR² = {r2:.4f}  → le modèle explique {r2*100:.1f}% de la variance du prix")
print(f"R² ajusté = {r2_adj:.4f}  (pénalise les paramètres superflus)")
Décomposition de la variance :
  SS total     = 1130398.0
  SS régression= 1087691.2
  SS résidus   = 42706.9

R² = 0.9622  → le modèle explique 96.2% de la variance du prix
R² ajusté = 0.9617  (pénalise les paramètres superflus)

Limites du R²

  • Un R² élevé ne garantit pas que le modèle est correct (les résidus peuvent montrer des patterns)

  • Un R² faible ne signifie pas que la relation est sans intérêt pratique (en sciences sociales, R² = 0,10 peut être remarquable)

  • Le R² augmente mécaniquement avec le nombre de prédicteurs — utiliser le R² ajusté pour comparer des modèles

  • Le R² ne mesure pas la qualité des prédictions hors-échantillon

Intervalles de confiance et de prédiction#

Deux types d’intervalles répondent à des questions différentes :

  • IC sur la moyenne \(E[Y|X=x^*]\) : où se situe la valeur moyenne pour \(X = x^*\) ?

  • IP (Intervalle de Prédiction) sur une nouvelle observation \(Y^*\) : quel est le prix d”une maison spécifique de surface \(x^*\) ?

L’IP est toujours plus large que l’IC, car il incorpore l’incertitude sur la valeur individuelle en plus de l’incertitude sur la moyenne.

Hide code cell source

# Calculer les IC et IP pour toute la plage de surface
surface_pred = np.linspace(surface.min() - 5, surface.max() + 5, 200)
df_pred = pd.DataFrame({'surface': surface_pred})

predictions = model.get_prediction(df_pred)
pred_summary = predictions.summary_frame(alpha=0.05)

fig, ax = plt.subplots(figsize=(10, 5))

# Données
ax.scatter(surface, prix, alpha=0.5, color='steelblue', s=40, label='Données', zorder=5)

# Droite ajustée
ax.plot(surface_pred, pred_summary['mean'], 'crimson', lw=2.5, label='Droite ajustée')

# IC sur la moyenne (bande étroite)
ax.fill_between(surface_pred,
                pred_summary['mean_ci_lower'],
                pred_summary['mean_ci_upper'],
                alpha=0.3, color='crimson',
                label='IC 95% sur la moyenne')

# IP (bande large)
ax.fill_between(surface_pred,
                pred_summary['obs_ci_lower'],
                pred_summary['obs_ci_upper'],
                alpha=0.15, color='navy',
                label='IP 95% (prédiction individuelle)')

# Annotation d'un point spécifique
x_ex = 100
y_ex = model.predict({'surface': [x_ex]})[0]
ic_ex = predictions.conf_int(alpha=0.05)[np.searchsorted(surface_pred, x_ex)]
ax.axvline(x_ex, color='gray', linestyle=':', lw=1.5, alpha=0.7)
ax.plot(x_ex, y_ex, 'k*', markersize=14, zorder=10, label=f'x = {x_ex} m² → ŷ = {y_ex:.0f} k€')

ax.set_xlabel('Surface (m²)')
ax.set_ylabel('Prix (k€)')
ax.set_title('Régression linéaire simple\navec intervalles de confiance et de prédiction')
ax.legend(fontsize=9, loc='upper left')
plt.tight_layout()
plt.savefig('_static/09_regression.png', bbox_inches='tight')
plt.show()

# Largeurs comparées
idx_mid = len(surface_pred) // 2
ic_width = pred_summary['mean_ci_upper'].iloc[idx_mid] - pred_summary['mean_ci_lower'].iloc[idx_mid]
ip_width = pred_summary['obs_ci_upper'].iloc[idx_mid] - pred_summary['obs_ci_lower'].iloc[idx_mid]
print(f"À surface = {surface_pred[idx_mid]:.0f} m² :")
print(f"  Largeur IC (moyenne) = {ic_width:.1f} k€")
print(f"  Largeur IP (individuel) = {ip_width:.1f} k€")
print(f"  L'IP est {ip_width/ic_width:.1f}× plus large que l'IC")
_images/63b2c83fafc8aba86b5e931a591d354b51ff2ff6e6dcc2059f5973071e842502.png
À surface = 89 m² :
  Largeur IC (moyenne) = 10.4 k€
  Largeur IP (individuel) = 93.8 k€
  L'IP est 9.0× plus large que l'IC

Hypothèses du modèle linéaire#

Les quatre hypothèses classiques (acronyme LINE) :

  1. Linéarité : \(E[Y|X] = \beta_0 + \beta_1 X\) (la relation est linéaire)

  2. Indépendance : les résidus sont indépendants

  3. Normalité : les résidus suivent \(\mathcal{N}(0, \sigma^2)\)

  4. Egalité des variances (homoscédasticité) : \(\text{Var}(\varepsilon_i) = \sigma^2\) constant

Ces hypothèses sont testables via les diagnostics graphiques et formels.

Diagnostics graphiques#

Les quatre graphiques de diagnostic classiques révèlent les violations des hypothèses :

Hide code cell source

def plot_diagnostics(model, titre=''):
    """Produit les 4 graphiques de diagnostic standard."""
    fig, axes = plt.subplots(2, 2, figsize=(12, 9))

    fitted = model.fittedvalues
    resid = model.resid
    resid_std = model.get_influence().resid_studentized_internal
    leverage = model.get_influence().hat_matrix_diag
    cooks = model.get_influence().cooks_distance[0]
    n = len(resid)

    # 1. Résidus vs valeurs ajustées
    ax = axes[0, 0]
    ax.scatter(fitted, resid, alpha=0.5, color='steelblue', s=40)
    ax.axhline(0, color='crimson', lw=2, linestyle='--')
    # Courbe de tendance locale
    order = np.argsort(fitted)
    z = np.polyfit(fitted, resid, 2)
    x_line = np.linspace(fitted.min(), fitted.max(), 200)
    ax.plot(x_line, np.polyval(z, x_line), 'orange', lw=2, linestyle='--', alpha=0.8)
    ax.set_xlabel('Valeurs ajustées (ŷ)')
    ax.set_ylabel('Résidus')
    ax.set_title('Résidus vs Valeurs ajustées\n(détecte non-linéarité, hétéroscédasticité)')

    # 2. QQ-plot des résidus standardisés
    ax = axes[0, 1]
    pp = ProbPlot(resid_std)
    theoretical_q = pp.theoretical_quantiles
    sample_q = pp.sample_quantiles
    ax.scatter(theoretical_q, sample_q, alpha=0.5, color='steelblue', s=40)
    lim = max(abs(theoretical_q).max(), abs(sample_q).max()) * 1.1
    ax.plot([-lim, lim], [-lim, lim], 'crimson', lw=2, label='Ligne de référence')
    ax.set_xlabel('Quantiles théoriques (normale)')
    ax.set_ylabel('Quantiles standardisés des résidus')
    ax.set_title('QQ-plot des résidus\n(détecte les écarts à la normalité)')
    ax.legend()

    # 3. Scale-location (racine carrée des résidus standardisés)
    ax = axes[1, 0]
    sqrt_resid_std = np.sqrt(np.abs(resid_std))
    ax.scatter(fitted, sqrt_resid_std, alpha=0.5, color='steelblue', s=40)
    z2 = np.polyfit(fitted, sqrt_resid_std, 1)
    ax.plot(x_line, np.polyval(z2, np.linspace(fitted.min(), fitted.max(), 200)),
            'orange', lw=2, linestyle='--')
    ax.set_xlabel('Valeurs ajustées (ŷ)')
    ax.set_ylabel('√|Résidus standardisés|')
    ax.set_title('Scale-Location\n(détecte l\'hétéroscédasticité)')

    # 4. Résidus vs Leverage (Cook's distance)
    ax = axes[1, 1]
    ax.scatter(leverage, resid_std, alpha=0.5, color='steelblue', s=40)
    ax.axhline(0, color='gray', lw=1, linestyle=':')
    ax.axhline(2, color='orange', lw=1.5, linestyle='--')
    ax.axhline(-2, color='orange', lw=1.5, linestyle='--')
    # Courbes de Cook's distance
    x_lev = np.linspace(leverage.min(), leverage.max(), 200)
    for d_cook in [0.5, 1.0]:
        p_params = model.df_model + 1
        cook_line = np.sqrt(d_cook * p_params * (1 - x_lev) / x_lev)
        ax.plot(x_lev, cook_line, 'crimson', lw=1.5, linestyle=':', alpha=0.7,
                label=f"Cook's d = {d_cook}")
        ax.plot(x_lev, -cook_line, 'crimson', lw=1.5, linestyle=':', alpha=0.7)
    # Identifier les points influents
    idx_influents = np.where(cooks > 4 / n)[0]
    if len(idx_influents) > 0:
        ax.scatter(leverage[idx_influents], resid_std[idx_influents],
                   color='crimson', s=80, zorder=5)
        for idx in idx_influents[:3]:
            ax.annotate(str(idx), (leverage[idx], resid_std[idx]),
                        textcoords='offset points', xytext=(5, 5), fontsize=8)
    ax.set_xlabel('Effet de levier (leverage)')
    ax.set_ylabel('Résidus standardisés')
    ax.set_title('Résidus vs Leverage\n(détecte les observations influentes)')
    ax.legend(fontsize=8)

    if titre:
        fig.suptitle(titre, fontsize=13, y=1.01)
    plt.tight_layout()
    return fig

fig = plot_diagnostics(model, titre='Diagnostics — Régression Prix ~ Surface (modèle valide)')
plt.savefig('_static/09_diagnostics.png', bbox_inches='tight')
plt.show()
_images/0b005a960bfa93dd454abef03c42d831d1e8c00f9079719a9d05e4faae1b930d.png
# Test formel d'homoscédasticité : Breusch-Pagan
bp_stat, bp_p, _, _ = het_breuschpagan(model.resid, model.model.exog)
print(f"Test de Breusch-Pagan (homoscédasticité) :")
print(f"  χ² = {bp_stat:.3f}, p = {bp_p:.4f}")
print(f"  → {'Hétéroscédasticité détectée' if bp_p < 0.05 else 'Pas d\'hétéroscédasticité détectée'}")
print()

# Test de normalité des résidus
stat_sw, p_sw = stats.shapiro(model.resid)
print(f"Test de Shapiro-Wilk sur les résidus :")
print(f"  W = {stat_sw:.4f}, p = {p_sw:.4f}")
print(f"  → {'Non-normalité détectée' if p_sw < 0.05 else 'Résidus compatible avec la normalité'}")
Test de Breusch-Pagan (homoscédasticité) :
  χ² = 0.689, p = 0.4066
  → Pas d'hétéroscédasticité détectée

Test de Shapiro-Wilk sur les résidus :
  W = 0.9868, p = 0.5850
  → Résidus compatible avec la normalité

Violations des hypothèses : exemples#

Hide code cell source

# Générer trois exemples pathologiques
rng_path = np.random.default_rng(99)
n_path = 60
x_p = np.linspace(1, 10, n_path)

# 1. Relation non linéaire (quadratique)
y_nonlin = 2 * x_p**2 - 5 * x_p + 10 + rng_path.normal(0, 5, n_path)
model_nonlin = smf.ols('y ~ x', data=pd.DataFrame({'y': y_nonlin, 'x': x_p})).fit()

# 2. Hétéroscédasticité (variance croissante)
y_het = 3 * x_p + rng_path.normal(0, x_p * 0.8, n_path)
model_het = smf.ols('y ~ x', data=pd.DataFrame({'y': y_het, 'x': x_p})).fit()

# 3. Outliers influents
y_out = 2 * x_p + 5 + rng_path.normal(0, 2, n_path)
y_out[55] = 60   # outlier en y
y_out_arr = y_out.copy()
x_arr = x_p.copy()
x_arr[55] = 9.5  # et avec levier élevé
model_out = smf.ols('y ~ x', data=pd.DataFrame({'y': y_out, 'x': x_p})).fit()
model_out2 = smf.ols('y ~ x', data=pd.DataFrame({'y': y_out_arr, 'x': x_arr})).fit()

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

# Ligne 1 : données + droite ajustée
for i, (x_d, y_d, mod, title) in enumerate([
    (x_p, y_nonlin, model_nonlin, 'Non-linéarité'),
    (x_p, y_het, model_het, 'Hétéroscédasticité'),
    (x_arr, y_out_arr, model_out2, 'Valeur influente'),
]):
    ax = axes[0, i]
    ax.scatter(x_d, y_d, alpha=0.5, color='steelblue', s=40)
    x_line = np.linspace(x_d.min(), x_d.max(), 200)
    ax.plot(x_line, mod.params['Intercept'] + mod.params['x'] * x_line, 'crimson', lw=2.5)
    if i == 2:
        ax.scatter(x_arr[55], y_out_arr[55], color='gold', s=200, zorder=5,
                   edgecolors='black', linewidth=2, label='Outlier influent')
        ax.legend()
    ax.set_title(f'Violation : {title}')
    ax.set_xlabel('X')
    ax.set_ylabel('Y')

# Ligne 2 : résidus vs fitted
for i, (mod, title) in enumerate([
    (model_nonlin, 'Résidus : non-linéarité → pattern'),
    (model_het, 'Résidus : hétéroscédasticité → éventail'),
    (model_out2, 'Résidus : outlier influent'),
]):
    ax = axes[1, i]
    ax.scatter(mod.fittedvalues, mod.resid, alpha=0.5, color='steelblue', s=40)
    ax.axhline(0, color='crimson', lw=2, linestyle='--')
    z = np.polyfit(mod.fittedvalues, mod.resid, 2)
    xr = np.linspace(mod.fittedvalues.min(), mod.fittedvalues.max(), 200)
    ax.plot(xr, np.polyval(z, xr), 'orange', lw=2, linestyle='--')
    ax.set_xlabel('Valeurs ajustées')
    ax.set_ylabel('Résidus')
    ax.set_title(title)

plt.suptitle('Diagnostics de violations des hypothèses', fontsize=13, y=1.01)
plt.tight_layout()
plt.savefig('_static/09_violations.png', bbox_inches='tight')
plt.show()
_images/ce013a57ae6fc76765c515564e8cd65230666f22a26015051323dd635f2052b5.png

Valeurs influentes et outliers#

Effet de levier et distance de Cook#

  • Effet de levier \(h_{ii}\) : mesure à quel point \(x_i\) est éloigné du centre du nuage. Un \(h_{ii} > 2p/n\) est considéré comme élevé.

  • Distance de Cook \(D_i\) : combine le levier et le résidu standardisé. Elle mesure l’influence de l’observation \(i\) sur l’ensemble des coefficients estimés.

\[D_i = \frac{\hat{\varepsilon}_i^2}{p \cdot \hat{\sigma}^2} \cdot \frac{h_{ii}}{(1 - h_{ii})^2}\]
# Analyse des observations influentes
influence = model.get_influence()
summary_inf = influence.summary_frame()

# Identifier les cas remarquables
seuil_cook = 4 / n
seuil_leverage = 2 * 2 / n  # 2p/n

obs_influentes = summary_inf[summary_inf['cooks_d'] > seuil_cook]
obs_levier = summary_inf[summary_inf['hat_diag'] > seuil_leverage]

print(f"Seuil Cook's distance : {seuil_cook:.4f}")
print(f"Seuil levier : {seuil_leverage:.4f}")
print(f"\n{len(obs_influentes)} observations avec Cook's d > seuil :")
if len(obs_influentes) > 0:
    print(obs_influentes[['hat_diag', 'cooks_d', 'student_resid']].head(10))
else:
    print("  Aucune")

print(f"\n{len(obs_levier)} observations à fort levier :")
print(obs_levier[['hat_diag', 'cooks_d', 'student_resid']].head(5))
Seuil Cook's distance : 0.0500
Seuil levier : 0.0500

3 observations avec Cook's d > seuil :
    hat_diag   cooks_d  student_resid
31  0.048896  0.102039      -2.031949
61  0.012715  0.069424       3.514201
74  0.046916  0.053155      -1.480771

3 observations à fort levier :
    hat_diag   cooks_d  student_resid
5   0.050189  0.007220      -0.520298
51  0.055183  0.011504      -0.625181
68  0.051266  0.002325      -0.291637

Hide code cell source

# Influence plot
fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))

# Cook's distance
ax = axes[0]
cooks_d = influence.cooks_distance[0]
ax.stem(range(n), cooks_d, markerfmt='C0o', linefmt='C0-', basefmt='k-')
ax.axhline(seuil_cook, color='crimson', lw=2, linestyle='--',
           label=f'Seuil 4/n = {seuil_cook:.3f}')
ax.set_xlabel('Indice de l\'observation')
ax.set_ylabel('Distance de Cook')
ax.set_title('Distance de Cook\n(observations influentes sur les coefficients)')
ax.legend()
# Annoter les observations au-dessus du seuil
for i, d in enumerate(cooks_d):
    if d > seuil_cook:
        ax.annotate(str(i), (i, d), textcoords='offset points', xytext=(5, 5), fontsize=8)

# Levier vs résidu standardisé
ax = axes[1]
leverage_vals = influence.hat_matrix_diag
resid_std_vals = influence.resid_studentized_internal
colors_inf = ['crimson' if d > seuil_cook else 'steelblue' for d in cooks_d]
sc = ax.scatter(leverage_vals, resid_std_vals, c=colors_inf, s=60, alpha=0.7)
ax.axhline(2, color='gray', linestyle=':', lw=1.5)
ax.axhline(-2, color='gray', linestyle=':', lw=1.5)
ax.axvline(seuil_leverage, color='gray', linestyle=':', lw=1.5)
ax.set_xlabel('Levier (h_ii)')
ax.set_ylabel('Résidu studentisé')
ax.set_title('Graphique Levier-Résidu\n(rouge = Cook\'s d > seuil)')
ax.set_xlim(0, leverage_vals.max() * 1.2)

# Quadrants
ax.text(ax.get_xlim()[1] * 0.7, 3, 'Fort levier\nGrand résidu', ha='center', fontsize=8, color='crimson')
ax.text(ax.get_xlim()[1] * 0.7, -3, 'Fort levier\nPetit résidu', ha='center', fontsize=8, color='gray')

plt.tight_layout()
plt.savefig('_static/09_influence.png', bbox_inches='tight')
plt.show()
_images/079f32b36f6587c2b907e735054ad435b5d8a1e8451b808a6a9088e7d528cf75.png

Analyse complète : exemple biologique#

Appliquons l’ensemble de la démarche sur un exemple biologique : concentration d’un médicament en fonction du temps après injection.

Hide code cell source

# Exemple : pharmacocinétique — concentration vs temps (après log-transformation)
rng_bio = np.random.default_rng(77)
n_bio = 45
temps = rng_bio.uniform(0.5, 8, n_bio)
# Relation log-linéaire : ln(conc) = 4.5 - 0.6 * temps + bruit
log_conc = 4.5 - 0.6 * temps + rng_bio.normal(0, 0.3, n_bio)
conc = np.exp(log_conc)  # concentration en ng/mL

df_pharma = pd.DataFrame({'temps': temps, 'conc': conc, 'log_conc': log_conc})

# Modèle sur les données brutes (non transformées) — mauvais
model_brut = smf.ols('conc ~ temps', data=df_pharma).fit()
# Modèle sur les données transformées — correct
model_log = smf.ols('log_conc ~ temps', data=df_pharma).fit()

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

# Données brutes
ax = axes[0, 0]
ax.scatter(temps, conc, alpha=0.6, color='steelblue', s=50)
x_l = np.linspace(temps.min(), temps.max(), 200)
ax.plot(x_l, model_brut.params['Intercept'] + model_brut.params['temps'] * x_l, 'crimson', lw=2)
ax.set_title(f'Données brutes\n(R² = {model_brut.rsquared:.3f} — mal ajusté)')
ax.set_xlabel('Temps (h)')
ax.set_ylabel('Concentration (ng/mL)')

# Données log-transformées
ax = axes[0, 1]
ax.scatter(temps, log_conc, alpha=0.6, color='steelblue', s=50)
ax.plot(x_l, model_log.params['Intercept'] + model_log.params['temps'] * x_l, 'crimson', lw=2)
ax.set_title(f'Log-transformation\n(R² = {model_log.rsquared:.3f} — bien ajusté)')
ax.set_xlabel('Temps (h)')
ax.set_ylabel('ln(Concentration)')

# Tableau des résultats
ax = axes[0, 2]
ax.axis('off')
results_data = [
    ['Paramètre', 'Modèle brut', 'Modèle log'],
    ['β₀', f"{model_brut.params['Intercept']:.2f}", f"{model_log.params['Intercept']:.2f}"],
    ['β₁', f"{model_brut.params['temps']:.2f}", f"{model_log.params['temps']:.2f}"],
    ['R²', f"{model_brut.rsquared:.3f}", f"{model_log.rsquared:.3f}"],
    ['AIC', f"{model_brut.aic:.1f}", f"{model_log.aic:.1f}"],
    ['Shapiro p', f"{stats.shapiro(model_brut.resid)[1]:.3f}",
     f"{stats.shapiro(model_log.resid)[1]:.3f}"],
]
table = ax.table(cellText=results_data[1:], colLabels=results_data[0],
                 loc='center', cellLoc='center')
table.auto_set_font_size(False)
table.set_fontsize(10)
table.scale(1.2, 2)
for j in range(3):
    table[0, j].set_facecolor('#2c3e50')
    table[0, j].set_text_props(color='white', fontweight='bold')
ax.set_title('Comparaison des modèles', pad=15)

# Diagnostics du modèle log
for i, (mod, title) in enumerate([
    (model_brut, 'Résidus bruts : non-linéarité'),
    (model_log, 'Résidus log-transformés : OK'),
]):
    ax = axes[1, i]
    ax.scatter(mod.fittedvalues, mod.resid, alpha=0.5, color='steelblue', s=40)
    ax.axhline(0, color='crimson', lw=2, linestyle='--')
    z = np.polyfit(mod.fittedvalues, mod.resid, 2)
    xr = np.linspace(mod.fittedvalues.min(), mod.fittedvalues.max(), 200)
    ax.plot(xr, np.polyval(z, xr), 'orange', lw=2, linestyle='--', alpha=0.8)
    ax.set_xlabel('Valeurs ajustées')
    ax.set_ylabel('Résidus')
    ax.set_title(title)

# QQ-plot du modèle log
ax = axes[1, 2]
resid_std_log = model_log.get_influence().resid_studentized_internal
pp = ProbPlot(resid_std_log)
ax.scatter(pp.theoretical_quantiles, pp.sample_quantiles,
           alpha=0.5, color='steelblue', s=40)
lim_qq = max(abs(pp.theoretical_quantiles).max(), abs(pp.sample_quantiles).max()) * 1.1
ax.plot([-lim_qq, lim_qq], [-lim_qq, lim_qq], 'crimson', lw=2)
ax.set_xlabel('Quantiles théoriques')
ax.set_ylabel('Quantiles observés')
ax.set_title('QQ-plot résidus — modèle log')

plt.suptitle('Pharmacocinétique : importance de la transformation', fontsize=13, y=1.01)
plt.tight_layout()
plt.savefig('_static/09_pharma.png', bbox_inches='tight')
plt.show()

print(f"Interprétation du modèle log :")
print(f"  β₁ = {model_log.params['temps']:.3f}")
print(f"  → La concentration diminue de {abs(model_log.params['temps'])*100:.1f}% par heure")
print(f"  (demi-vie estimée : {np.log(2)/abs(model_log.params['temps']):.2f}h)")
_images/793b5a84b58a88cad346d4e8646415c57b1584a5fd75b39c090a21f0e9258609.png
Interprétation du modèle log :
  β₁ = -0.575
  → La concentration diminue de 57.5% par heure
  (demi-vie estimée : 1.21h)

Checklist d’une régression linéaire simple

  1. Explorer les données : nuage de points, corrélation, outliers évidents

  2. Transformer si nécessaire (log, racine carrée) pour linéariser la relation

  3. Ajuster le modèle MCO

  4. Examiner le summary : coefficients, IC, R², test F global

  5. Vérifier les diagnostics : résidus vs fitted, QQ-plot, scale-location, leverage

  6. Tester formellement : Shapiro-Wilk (normalité), Breusch-Pagan (homoscédasticité)

  7. Identifier les observations influentes (Cook’s distance, levier)

  8. Interpréter avec IC et tailles d’effet, pas seulement p-valeurs

  9. Distinguer IC sur la moyenne et IP sur une nouvelle observation