---
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 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.

```{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.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

```{code-cell} python3
:tags: [hide-input]

# 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()
```

## 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.

```{code-cell} python3
# Ajustement MCO avec statsmodels
model = smf.ols('prix ~ surface', data=df_immo).fit()
print(model.summary())
```

```{code-cell} python3
# 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}]")
```

## 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}$$

```{code-cell} python3
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)")
```

```{admonition} Limites du R²
:class: warning

- 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.

```{code-cell} python3
:tags: [hide-input]

# 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")
```

## Hypothèses du modèle linéaire

Les quatre hypothèses classiques (acronyme LINE) :

1. **L**inéarité : $E[Y|X] = \beta_0 + \beta_1 X$ (la relation est linéaire)
2. **I**ndépendance : les résidus sont indépendants
3. **N**ormalité : les résidus suivent $\mathcal{N}(0, \sigma^2)$
4. **E**galité 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 :

```{code-cell} python3
:tags: [hide-input]

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()
```

```{code-cell} python3
# 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é'}")
```

## Violations des hypothèses : exemples

```{code-cell} python3
:tags: [hide-input]

# 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()
```

## 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}$$

```{code-cell} python3
# 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))
```

```{code-cell} python3
:tags: [hide-input]

# 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()
```

## 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.

```{code-cell} python3
:tags: [hide-input]

# 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)")
```

```{admonition} Checklist d'une régression linéaire simple
:class: tip

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
```