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

> Quand la réponse est oui ou non, le modèle doit l'être aussi.

## Introduction

La régression linéaire suppose une variable réponse continue et non bornée. Mais de nombreux problèmes pratiques impliquent des résultats binaires : un patient est malade ou sain, un crédit est remboursé ou en défaut, un email est spam ou légitime. Appliquer une régression linéaire à une variable indicatrice $Y \in \{0, 1\}$ produit des prédictions en dehors de $[0, 1]$, viole l'hypothèse de normalité des résidus, et n'est pas interprétable comme une probabilité.

La **régression logistique** résout ce problème en modélisant directement $P(Y = 1 \mid X)$ via une transformation qui confine les prédictions dans $(0, 1)$.

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

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import seaborn as sns
import pandas as pd
from scipy import stats
from scipy.special import expit  # fonction logistique

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

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# 1. Problème de la régression linéaire sur binaire
ax = axes[0]
np.random.seed(42)
x = np.linspace(-3, 3, 200)
y_bin = (x + np.random.normal(0, 1, 200) > 0).astype(int)
# Régression linéaire
from numpy.polynomial import polynomial as P
coef = np.polyfit(x, y_bin, 1)
y_lin = np.polyval(coef, x)
ax.scatter(x, y_bin, alpha=0.3, color='steelblue', s=20, label='Données')
ax.plot(x, y_lin, 'tomato', lw=2, label='Régression linéaire')
ax.axhline(0, color='gray', lw=0.8, ls='--')
ax.axhline(1, color='gray', lw=0.8, ls='--')
ax.fill_between(x, 0, 1, alpha=0.08, color='green')
ax.set_xlabel('x'); ax.set_ylabel('Y')
ax.set_title('Problème : RLin sur binaire\n(prédictions hors [0,1])')
ax.legend(fontsize=8)

# 2. Fonction logistique (sigmoïde)
ax = axes[1]
z = np.linspace(-6, 6, 300)
sigma = expit(z)
ax.plot(z, sigma, 'steelblue', lw=2.5)
ax.axhline(0.5, color='tomato', ls='--', lw=1.5, label='p = 0.5 (z = 0)')
ax.axvline(0, color='tomato', ls='--', lw=1.5)
ax.fill_between(z, 0, sigma, alpha=0.12, color='steelblue')
ax.set_xlabel('z = β₀ + β₁x'); ax.set_ylabel('σ(z) = P(Y=1|x)')
ax.set_title('Fonction sigmoïde\nσ(z) = 1/(1+e⁻ᶻ)')
ax.set_ylim(-0.05, 1.05)
ax.legend(fontsize=9)

# 3. Régression logistique correcte
ax = axes[2]
y_log = expit(1.5 * x)
ax.scatter(x, y_bin, alpha=0.3, color='steelblue', s=20, label='Données')
ax.plot(x, y_log, 'green', lw=2.5, label='Régression logistique')
ax.axhline(0, color='gray', lw=0.8, ls='--')
ax.axhline(1, color='gray', lw=0.8, ls='--')
ax.set_xlabel('x'); ax.set_ylabel('P(Y=1|x)')
ax.set_title('Solution : Régression logistique\n(prédictions ∈ (0,1))')
ax.legend(fontsize=8)

plt.tight_layout()
plt.savefig('_static/11_intro.png', dpi=120, bbox_inches='tight')
plt.show()
```

## La fonction logistique et le modèle

### Du linéaire aux probabilités : le logit

Le modèle de régression logistique postule :

$$P(Y = 1 \mid \mathbf{x}) = \sigma(\mathbf{x}^T\boldsymbol{\beta}) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p)}}$$

La transformation inverse s'appelle le **logit** (ou log-odds) :

$$\text{logit}(p) = \log\frac{p}{1-p} = \beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p$$

Le ratio $p/(1-p)$ est la **cote** (*odds*) : si $p = 0{,}75$, la cote est $3$ (trois chances pour une). Le logit est la transformation qui rend ce rapport linéaire en $\mathbf{x}$.

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

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Relation p <-> logit
p_vals = np.linspace(0.001, 0.999, 500)
logit_vals = np.log(p_vals / (1 - p_vals))

ax = axes[0]
ax.plot(p_vals, logit_vals, 'steelblue', lw=2.5)
ax.axhline(0, color='tomato', ls='--', lw=1.2, label='logit = 0 ↔ p = 0.5')
ax.axvline(0.5, color='tomato', ls='--', lw=1.2)
ax.set_xlabel('Probabilité p'); ax.set_ylabel('logit(p) = log(p/(1-p))')
ax.set_title('Transformation logit')
ax.legend(fontsize=9)

# Odds ratio illustration
ax = axes[1]
betas = np.array([-2, -1, 0, 1, 2])
x_plot = np.linspace(-3, 3, 200)
colors = sns.color_palette("coolwarm", len(betas))
for b, c in zip(betas, colors):
    ax.plot(x_plot, expit(b * x_plot), color=c, lw=2, label=f'β₁ = {b}')
ax.axhline(0.5, color='gray', ls=':', lw=1)
ax.set_xlabel('x'); ax.set_ylabel('P(Y=1|x)')
ax.set_title('Effet du coefficient β₁\nsur la courbe sigmoïde')
ax.legend(fontsize=8, loc='upper left')

plt.tight_layout()
plt.savefig('_static/11_logit.png', dpi=120, bbox_inches='tight')
plt.show()
```

## Estimation par maximum de vraisemblance

Contrairement à la régression linéaire (qui minimise les moindres carrés), la régression logistique est estimée par **maximum de vraisemblance** (MLE). La log-vraisemblance est :

$$\ell(\boldsymbol{\beta}) = \sum_{i=1}^n \left[ y_i \log \hat{p}_i + (1 - y_i) \log(1 - \hat{p}_i) \right]$$

où $\hat{p}_i = \sigma(\mathbf{x}_i^T\boldsymbol{\beta})$. Il n'existe pas de solution analytique fermée ; on utilise des méthodes itératives (Newton-Raphson, IRLS).

```{admonition} Pourquoi pas les MCO ?
:class: note

Les moindres carrés ordinaires supposent des résidus gaussiens homoscédastiques. Pour une variable binaire, la variance est $p(1-p)$, qui dépend de $p$ lui-même (hétéroscédasticité structurelle), et les résidus ne sont pas gaussiens. Le MLE est l'approche naturelle et efficiente.
```

## Exemple fil rouge : risque de défaut de crédit

```{code-cell} python
import statsmodels.api as sm
import statsmodels.formula.api as smf
from sklearn.metrics import (confusion_matrix, classification_report,
                              roc_curve, auc, precision_recall_curve)
from sklearn.model_selection import train_test_split

# Simulation d'un jeu de données de crédit
np.random.seed(2024)
n = 800

age = np.random.normal(40, 10, n).clip(20, 70)
revenu = np.random.lognormal(10.5, 0.5, n)  # revenus annuels
dette_revenu = np.random.beta(2, 5, n) * 0.8  # ratio dette/revenu
historique = np.random.choice([0, 1], n, p=[0.3, 0.7])  # bon historique

# Log-odds vrai
log_odds = (-1.5
            + 0.03 * (age - 40)
            - 0.4 * np.log(revenu / 50000)
            + 3.0 * dette_revenu
            - 1.2 * historique)
p_defaut = expit(log_odds)
defaut = np.random.binomial(1, p_defaut)

df = pd.DataFrame({
    'defaut': defaut,
    'age': age.round(1),
    'revenu_k': (revenu / 1000).round(1),
    'dette_revenu': dette_revenu.round(3),
    'bon_historique': historique
})

print(f"Taux de défaut : {defaut.mean():.1%}")
print(f"\nAperçu du jeu de données :")
df.head()
```

```{code-cell} python
# Ajustement du modèle avec statsmodels
model = smf.logit(
    'defaut ~ age + np.log(revenu_k) + dette_revenu + bon_historique',
    data=df
).fit(disp=False)

print(model.summary())
```

### Interprétation des coefficients : odds ratios

Le coefficient $\beta_j$ représente la variation du log-odds associée à une augmentation unitaire de $x_j$. L'**odds ratio** (OR) est $e^{\beta_j}$ :

- OR > 1 : le facteur **augmente** la cote de défaut
- OR < 1 : le facteur **diminue** la cote de défaut
- OR = 1 : pas d'effet

```{code-cell} python
# Odds ratios avec intervalles de confiance à 95%
conf = model.conf_int()
conf.columns = ['IC_inf', 'IC_sup']
resultats = pd.DataFrame({
    'Coefficient': model.params,
    'OR': np.exp(model.params),
    'OR_IC_inf': np.exp(conf['IC_inf']),
    'OR_IC_sup': np.exp(conf['IC_sup']),
    'p_valeur': model.pvalues
})
print(resultats.round(4))
```

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

# Visualisation des odds ratios avec IC
fig, ax = plt.subplots(figsize=(9, 5))

params = resultats.drop('Intercept')
y_pos = range(len(params))
colors_or = ['tomato' if v > 1 else 'steelblue' for v in params['OR']]

ax.barh(y_pos, params['OR'], xerr=[
    params['OR'] - params['OR_IC_inf'],
    params['OR_IC_sup'] - params['OR']
], color=colors_or, alpha=0.75, capsize=4, height=0.5)
ax.axvline(1, color='black', lw=1.5, ls='--', label='OR = 1 (pas d\'effet)')
ax.set_yticks(y_pos)
ax.set_yticklabels(['Âge', 'log(Revenu k€)', 'Dette/Revenu', 'Bon historique'])
ax.set_xlabel('Odds Ratio (avec IC 95%)')
ax.set_title('Odds ratios — Régression logistique\nModèle de risque de défaut de crédit')
ax.legend(fontsize=9)
for i, (v, p) in enumerate(zip(params['OR'], params['p_valeur'])):
    stars = '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else 'ns'))
    ax.text(v + 0.05, i, f'{v:.2f} {stars}', va='center', fontsize=9)

plt.tight_layout()
plt.savefig('_static/11_or.png', dpi=120, bbox_inches='tight')
plt.show()
```

## Évaluation du modèle

### Matrice de confusion et métriques

```{code-cell} python
# Prédictions (seuil = 0.5 par défaut)
df_train, df_test = train_test_split(df, test_size=0.3, random_state=42, stratify=df['defaut'])

model_train = smf.logit(
    'defaut ~ age + np.log(revenu_k) + dette_revenu + bon_historique',
    data=df_train
).fit(disp=False)

p_pred = model_train.predict(df_test)
y_pred = (p_pred >= 0.5).astype(int)
y_true = df_test['defaut'].values

cm = confusion_matrix(y_true, y_pred)
print("Matrice de confusion :")
print(pd.DataFrame(cm,
    index=['Réel : Non-défaut', 'Réel : Défaut'],
    columns=['Prédit : Non-défaut', 'Prédit : Défaut']))
print()
print(classification_report(y_true, y_pred,
    target_names=['Non-défaut', 'Défaut']))
```

### Courbe ROC et AUC

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

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

# Courbe ROC
fpr, tpr, thresholds_roc = roc_curve(y_true, p_pred)
roc_auc = auc(fpr, tpr)

ax = axes[0]
ax.plot(fpr, tpr, 'steelblue', lw=2.5, label=f'ROC (AUC = {roc_auc:.3f})')
ax.plot([0, 1], [0, 1], 'gray', ls='--', lw=1.5, label='Classifieur aléatoire')
ax.fill_between(fpr, tpr, alpha=0.15, color='steelblue')

# Annoter le seuil 0.5
idx = np.argmin(np.abs(thresholds_roc - 0.5))
ax.scatter(fpr[idx], tpr[idx], color='tomato', s=100, zorder=5,
           label=f'Seuil=0.5 (FPR={fpr[idx]:.2f}, TPR={tpr[idx]:.2f})')
ax.set_xlabel('Taux de faux positifs (FPR)')
ax.set_ylabel('Taux de vrais positifs (TPR = Rappel)')
ax.set_title('Courbe ROC')
ax.legend(fontsize=8)

# Courbe Précision-Rappel
precision, recall, thresholds_pr = precision_recall_curve(y_true, p_pred)
pr_auc = auc(recall, precision)

ax = axes[1]
ax.plot(recall, precision, 'darkgreen', lw=2.5, label=f'Précision-Rappel (AUC={pr_auc:.3f})')
ax.axhline(y_true.mean(), color='gray', ls='--', lw=1.5,
           label=f'Baseline ({y_true.mean():.2f})')
idx2 = np.argmin(np.abs(thresholds_pr - 0.5))
ax.scatter(recall[idx2], precision[idx2], color='tomato', s=100, zorder=5,
           label=f'Seuil=0.5')
ax.set_xlabel('Rappel'); ax.set_ylabel('Précision')
ax.set_title('Courbe Précision-Rappel')
ax.legend(fontsize=8)

plt.tight_layout()
plt.savefig('_static/11_roc.png', dpi=120, bbox_inches='tight')
plt.show()
```

### Impact du seuil de classification

```{admonition} Choix du seuil
:class: important

Le seuil à 0,5 n'est pas toujours optimal. Dans un contexte de détection de fraude ou de maladie grave, on préfère un seuil plus bas (augmenter le rappel quitte à baisser la précision). Le seuil optimal dépend des coûts relatifs des erreurs de type I et type II.
```

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

thresholds = np.linspace(0.05, 0.95, 200)
precisions, recalls, f1s = [], [], []

for t in thresholds:
    y_t = (p_pred >= t).astype(int)
    tp = ((y_t == 1) & (y_true == 1)).sum()
    fp = ((y_t == 1) & (y_true == 0)).sum()
    fn = ((y_t == 0) & (y_true == 1)).sum()
    prec = tp / (tp + fp + 1e-9)
    rec = tp / (tp + fn + 1e-9)
    f1 = 2 * prec * rec / (prec + rec + 1e-9)
    precisions.append(prec)
    recalls.append(rec)
    f1s.append(f1)

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(thresholds, precisions, 'steelblue', lw=2, label='Précision')
ax.plot(thresholds, recalls, 'tomato', lw=2, label='Rappel')
ax.plot(thresholds, f1s, 'darkgreen', lw=2, label='F1-score')
best_t = thresholds[np.argmax(f1s)]
ax.axvline(best_t, color='black', ls='--', lw=1.5,
           label=f'Seuil optimal F1 : {best_t:.2f}')
ax.axvline(0.5, color='orange', ls=':', lw=1.5, label='Seuil = 0.5')
ax.set_xlabel('Seuil de classification')
ax.set_ylabel('Score')
ax.set_title('Précision, Rappel et F1 en fonction du seuil')
ax.legend(fontsize=9)
plt.tight_layout()
plt.savefig('_static/11_seuil.png', dpi=120, bbox_inches='tight')
plt.show()
print(f"Seuil optimal (F1) : {best_t:.2f}")
```

## Régression logistique multinomiale

Quand la variable réponse est **multi-classes** (plus de 2 catégories non ordonnées), on étend la régression logistique au modèle **softmax** :

$$P(Y = k \mid \mathbf{x}) = \frac{e^{\mathbf{x}^T\boldsymbol{\beta}_k}}{\sum_{j=1}^K e^{\mathbf{x}^T\boldsymbol{\beta}_j}}$$

On fixe une classe de référence (par exemple $k = K$) et on estime $K-1$ jeux de coefficients.

```{code-cell} python
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import make_classification

# Exemple : classification 3 classes
np.random.seed(42)
X_multi, y_multi = make_classification(
    n_samples=600, n_features=2, n_redundant=0, n_informative=2,
    n_clusters_per_class=1, n_classes=3, random_state=42
)

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_multi)

clf_multi = LogisticRegression(solver='lbfgs', max_iter=500)
clf_multi.fit(X_scaled, y_multi)
print(f"Précision multinomiale : {clf_multi.score(X_scaled, y_multi):.3f}")
print(f"Coefficients (forme) : {clf_multi.coef_.shape}  (K classes × p features)")
```

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

# Frontières de décision
fig, ax = plt.subplots(figsize=(8, 6))
xx, yy = np.meshgrid(np.linspace(-3, 3, 300), np.linspace(-3, 3, 300))
Z = clf_multi.predict(np.c_[xx.ravel(), yy.ravel()]).reshape(xx.shape)

ax.contourf(xx, yy, Z, alpha=0.25, cmap='Set1')
scatter = ax.scatter(X_scaled[:, 0], X_scaled[:, 1], c=y_multi,
                     cmap='Set1', edgecolors='k', linewidths=0.4, s=40, alpha=0.8)
ax.set_xlabel('Feature 1 (standardisée)')
ax.set_ylabel('Feature 2 (standardisée)')
ax.set_title('Régression logistique multinomiale — Frontières de décision')
plt.colorbar(scatter, ax=ax, label='Classe')
plt.tight_layout()
plt.savefig('_static/11_multi.png', dpi=120, bbox_inches='tight')
plt.show()
```

## Régression logistique ordinale

Quand les catégories sont **ordonnées** (par exemple notes A/B/C/D, niveaux de satisfaction 1-5), le modèle **proportional odds** (logit cumulatif) est adapté :

$$\log \frac{P(Y \le k \mid \mathbf{x})}{P(Y > k \mid \mathbf{x})} = \alpha_k - \mathbf{x}^T\boldsymbol{\beta}$$

```{code-cell} python
# Simulation d'une variable ordinale (note de crédit : 1=mauvais, 5=excellent)
from statsmodels.miscmodels.ordinal_model import OrderedModel

np.random.seed(42)
n_ord = 500
x1 = np.random.normal(0, 1, n_ord)
x2 = np.random.normal(0, 1, n_ord)
latent = 0.8 * x1 - 0.5 * x2 + np.random.normal(0, 1, n_ord)
cuts = [-1.5, -0.5, 0.5, 1.5]
y_ord = np.digitize(latent, cuts)  # 0, 1, 2, 3, 4

df_ord = pd.DataFrame({'y': y_ord, 'x1': x1, 'x2': x2})
print("Distribution de la variable ordinale :")
print(df_ord['y'].value_counts().sort_index())

model_ord = OrderedModel(df_ord['y'], df_ord[['x1', 'x2']], distr='logit')
result_ord = model_ord.fit(method='bfgs', disp=False)
print(result_ord.summary())
```

## Hypothèses et diagnostics

```{admonition} Hypothèses du modèle logistique
:class: note

1. **Indépendance des observations** : pas de structure de dépendance (groupes, séries temporelles).
2. **Pas de séparation parfaite** : si un prédicteur sépare parfaitement les classes, les coefficients divergent vers ±∞.
3. **Pas de multicolinéarité sévère** : vérifier avec le VIF (comme en régression linéaire).
4. **Taille d'échantillon** : règle empirique — au moins 10 événements par variable (EPV ≥ 10).
5. **Linéarité du logit** : la relation entre les variables continues et le logit doit être linéaire.
```

```{code-cell} python
# VIF pour diagnostiquer la multicolinéarité
from statsmodels.stats.outliers_influence import variance_inflation_factor

X_vif = df[['age', 'revenu_k', 'dette_revenu', 'bon_historique']].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(X_vif.shape[1])]
})
print("Facteur d'inflation de variance (VIF) :")
print(vif_data.to_string(index=False))
print("\nVIF < 5 : acceptable | VIF > 10 : problématique")
```

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

# Diagnostic : courbe logit en fonction des variables continues
fig, axes = plt.subplots(1, 2, figsize=(12, 4))

vars_cont = ['age', 'revenu_k']
labels = ['Âge', 'Revenu (k€)']
for ax, var, label in zip(axes, vars_cont, labels):
    # Logit empirique par décile
    df_diag = df[['defaut', var]].copy()
    df_diag['decile'] = pd.qcut(df_diag[var], 10)
    grouped = df_diag.groupby('decile', observed=True)['defaut'].agg(['mean', 'count'])
    grouped['logit'] = np.log(grouped['mean'].clip(0.01, 0.99) /
                               (1 - grouped['mean'].clip(0.01, 0.99)))
    grouped['x_mid'] = grouped.index.map(lambda x: x.mid)

    ax.scatter(grouped['x_mid'], grouped['logit'],
               s=grouped['count']/5, color='steelblue', alpha=0.7)
    # Tendance linéaire
    coef_diag = np.polyfit(grouped['x_mid'], grouped['logit'], 1)
    x_range = np.linspace(grouped['x_mid'].min(), grouped['x_mid'].max(), 100)
    ax.plot(x_range, np.polyval(coef_diag, x_range), 'tomato', lw=2)
    ax.set_xlabel(label)
    ax.set_ylabel('Logit empirique')
    ax.set_title(f'Linéarité du logit : {label}\n(taille ∝ n dans le décile)')

plt.tight_layout()
plt.savefig('_static/11_diag.png', dpi=120, bbox_inches='tight')
plt.show()
```

## Résumé

```{admonition} Points clés — Régression logistique
:class: note

- La régression logistique modélise $P(Y=1|\mathbf{x})$ via la fonction sigmoïde, garantissant des prédictions dans $(0,1)$.
- Le modèle est estimé par **maximum de vraisemblance**, pas par MCO.
- Les **odds ratios** ($e^\beta$) quantifient l'effet multiplicatif sur la cote.
- L'évaluation passe par la matrice de confusion, la courbe ROC/AUC et la courbe précision-rappel.
- Le **seuil de classification** doit être choisi selon le coût des erreurs, pas automatiquement à 0,5.
- L'extension **multinomiale** (softmax) gère les multi-classes ; le modèle **ordinal** préserve l'ordre des catégories.
- Diagnostics indispensables : VIF, EPV, linéarité du logit, absence de séparation parfaite.
```