---
jupytext:
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.0
kernelspec:
  name: python3
  display_name: Python 3
---

# Modèles linéaires généralisés (GLM)

> Un seul cadre pour des réponses très différentes — comptages, proportions, coûts, durées.

## Introduction

La régression linéaire ordinaire repose sur deux hypothèses fortes : la réponse est continue et les résidus sont gaussiens homoscédastiques. En pratique, de nombreuses variables réponse violent ces hypothèses :

- **Comptages** (nombre d'accidents, d'appels, de sinistres) : entiers positifs, asymétriques
- **Coûts et durées** : continus mais strictement positifs et souvent très asymétriques
- **Proportions** : bornées dans $[0, 1]$

Les **modèles linéaires généralisés** (GLM, Nelder & Wedderburn, 1972) unifient régression linéaire, logistique, de Poisson et bien d'autres dans un cadre cohérent.

```{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
import statsmodels.api as sm
import statsmodels.formula.api as smf
import warnings
warnings.filterwarnings('ignore')

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

# Illustration : distributions de la famille exponentielle
fig, axes = plt.subplots(2, 3, figsize=(15, 8))

np.random.seed(42)
n_demo = 5000

# Gaussienne
ax = axes[0, 0]
x_gauss = np.random.normal(5, 2, n_demo)
ax.hist(x_gauss, bins=50, density=True, color='steelblue', alpha=0.7, edgecolor='white')
x_range = np.linspace(-2, 12, 300)
ax.plot(x_range, stats.norm.pdf(x_range, 5, 2), 'tomato', lw=2.5)
ax.set_title('Gaussienne\nRégression linéaire'); ax.set_xlabel('y')

# Binomiale
ax = axes[0, 1]
x_binom = np.random.binomial(1, 0.4, n_demo)
ax.bar([0, 1], [np.mean(x_binom==0), np.mean(x_binom==1)],
       color=['steelblue', 'tomato'], alpha=0.8, edgecolor='white')
ax.set_xticks([0, 1]); ax.set_xticklabels(['0', '1'])
ax.set_title('Bernoulli/Binomiale\nRégression logistique'); ax.set_xlabel('y')

# Poisson
ax = axes[0, 2]
x_pois = np.random.poisson(4, n_demo)
vals, counts = np.unique(x_pois, return_counts=True)
ax.bar(vals, counts/n_demo, color='steelblue', alpha=0.8, edgecolor='white')
x_range_p = np.arange(0, 15)
ax.plot(x_range_p, stats.poisson.pmf(x_range_p, 4), 'tomato', lw=2.5, marker='o', ms=4)
ax.set_title('Poisson\nRégression de Poisson'); ax.set_xlabel('y')

# Binomiale négative
ax = axes[1, 0]
x_nb = np.random.negative_binomial(3, 0.4, n_demo)
vals2, counts2 = np.unique(x_nb, return_counts=True)
ax.bar(vals2[:20], counts2[:20]/n_demo, color='steelblue', alpha=0.8, edgecolor='white')
ax.set_title('Binomiale Négative\nSur-dispersion'); ax.set_xlabel('y')

# Gamma
ax = axes[1, 1]
x_gamma = np.random.gamma(2, 5, n_demo)
ax.hist(x_gamma, bins=60, density=True, color='steelblue', alpha=0.7, edgecolor='white')
x_range_g = np.linspace(0, 50, 300)
ax.plot(x_range_g, stats.gamma.pdf(x_range_g, 2, scale=5), 'tomato', lw=2.5)
ax.set_title('Gamma\nCoûts, durées'); ax.set_xlabel('y')

# Tweedie (compound Poisson-Gamma)
ax = axes[1, 2]
# Approximation : mélange de zéros et gamma
zeros = np.zeros(int(n_demo * 0.4))
nonzeros = np.random.gamma(2, 300, n_demo - len(zeros))
x_tw = np.concatenate([zeros, nonzeros])
ax.hist(x_tw[x_tw > 0], bins=60, density=True, color='steelblue', alpha=0.7,
        edgecolor='white', label='Non-nul')
ax.bar(0, 0.4, width=30, color='tomato', alpha=0.5, label=f'Zéros ({0.4:.0%})')
ax.set_title('Tweedie\nSinistres assurance'); ax.set_xlabel('y')
ax.legend(fontsize=8)

plt.suptitle('Familles de distributions dans les GLM', fontsize=13, fontweight='bold', y=1.01)
plt.tight_layout()
plt.savefig('_static/12_familles.png', dpi=120, bbox_inches='tight')
plt.show()
```

## Structure d'un GLM

Un GLM est défini par trois composantes :

### 1. Composante aléatoire (famille)

La variable réponse $Y_i$ suit une distribution de la **famille exponentielle** :

$$f(y_i; \theta_i, \phi) = \exp\left[\frac{y_i \theta_i - b(\theta_i)}{a(\phi)} + c(y_i, \phi)\right]$$

Les principales familles sont : Gaussienne, Binomiale, Poisson, Gamma, Gaussienne inverse, Tweedie.

### 2. Composante systématique (prédicteur linéaire)

$$\eta_i = \mathbf{x}_i^T\boldsymbol{\beta} = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip}$$

### 3. Fonction de lien

La fonction de lien $g$ relie l'espérance $\mu_i = E[Y_i]$ au prédicteur linéaire :

$$g(\mu_i) = \eta_i$$

| Famille | Lien canonique | Liens alternatifs | Usage |
|---------|---------------|-------------------|-------|
| Gaussienne | Identité | — | Continu non borné |
| Binomiale | Logit | Probit, complémentaire log-log | Binaire/proportion |
| Poisson | Log | Racine carrée, identité | Comptages |
| Gamma | Inverse | Log, identité | Coûts, durées |
| Gaussienne inverse | Inverse carré | Log, identité | Extrêmement asymétrique |

## Régression de Poisson : modéliser les comptages

### Modèle et interprétation

Pour un comptage $Y_i \sim \text{Poisson}(\mu_i)$ :

$$\log(\mu_i) = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip}$$

Le coefficient $\beta_j$ représente la variation du log-espérance. L'**Incidence Rate Ratio** (IRR = $e^{\beta_j}$) est le ratio multiplicatif sur le taux attendu pour une augmentation unitaire de $x_j$.

```{code-cell} python
# Exemple : nombre d'accidents de travail
np.random.seed(2024)
n = 500

anciennete = np.random.uniform(0, 20, n)       # années d'expérience
formation = np.random.binomial(1, 0.6, n)       # formation reçue ou non
risque_secteur = np.random.choice([1, 2, 3], n, p=[0.4, 0.35, 0.25])  # 1=faible, 3=élevé

log_mu = (1.2
          - 0.06 * anciennete
          - 0.7 * formation
          + 0.5 * (risque_secteur - 1)
          + np.random.normal(0, 0.1, n))
mu = np.exp(log_mu)
accidents = np.random.poisson(mu)

df_pois = pd.DataFrame({
    'accidents': accidents,
    'anciennete': anciennete.round(1),
    'formation': formation,
    'risque_secteur': risque_secteur
})

print(f"Distribution des accidents :")
print(df_pois['accidents'].describe())
print(f"\nVariance / Moyenne = {df_pois['accidents'].var():.2f} / {df_pois['accidents'].mean():.2f}"
      f" = {df_pois['accidents'].var()/df_pois['accidents'].mean():.2f}")
print("(ratio proche de 1 → Poisson bien ajusté)")
```

```{code-cell} python
# Modèle Poisson avec statsmodels
model_pois = smf.glm(
    'accidents ~ anciennete + formation + C(risque_secteur)',
    data=df_pois,
    family=sm.families.Poisson()
).fit()

print(model_pois.summary())
```

```{code-cell} python
# IRR avec intervalles de confiance
conf_pois = model_pois.conf_int()
irr = pd.DataFrame({
    'IRR': np.exp(model_pois.params),
    'IRR_IC_inf': np.exp(conf_pois[0]),
    'IRR_IC_sup': np.exp(conf_pois[1]),
    'p_valeur': model_pois.pvalues
})
print("Incidence Rate Ratios :")
print(irr.round(4))
```

### Sur-dispersion et binomiale négative

L'hypothèse fondamentale de Poisson est $\text{Var}(Y) = E[Y]$ (équidispersion). En pratique, les comptages réels ont souvent $\text{Var}(Y) > E[Y]$ : c'est la **sur-dispersion**.

```{code-cell} python
# Simuler un jeu de données sur-dispersé
np.random.seed(42)
n_od = 400

x_od = np.random.normal(0, 1, n_od)
# Sur-dispersion via mélange (paramètre de dispersion r = 3)
mu_od = np.exp(1.0 + 0.8 * x_od)
r_od = 3  # paramètre de dispersion
p_od = r_od / (r_od + mu_od)
y_od = np.random.negative_binomial(r_od, p_od, n_od)

df_od = pd.DataFrame({'y': y_od, 'x': x_od})
print(f"Variance / Moyenne = {y_od.var():.1f} / {y_od.mean():.1f} = {y_od.var()/y_od.mean():.1f}")
print("(ratio >> 1 → sur-dispersion !!)")

# Comparer Poisson vs Binomiale Négative
model_pois_od = smf.glm('y ~ x', data=df_od, family=sm.families.Poisson()).fit()
model_nb = smf.glm('y ~ x', data=df_od,
                   family=sm.families.NegativeBinomial()).fit()

print(f"\nAIC Poisson         : {model_pois_od.aic:.1f}")
print(f"AIC Binomiale Neg.  : {model_nb.aic:.1f}")
print(f"(AIC plus bas = meilleur ajustement)")
```

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

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

# Comparaison Poisson vs NB
ax = axes[0]
x_pred = np.linspace(-3, 3, 200)
df_pred = pd.DataFrame({'x': x_pred})
mu_pois_pred = model_pois_od.predict(df_pred)
mu_nb_pred = model_nb.predict(df_pred)

ax.scatter(x_od, y_od, alpha=0.25, color='steelblue', s=20, label='Données')
ax.plot(x_pred, mu_pois_pred, 'tomato', lw=2.5, label='Poisson')
ax.plot(x_pred, mu_nb_pred, 'darkgreen', lw=2.5, ls='--', label='Binomiale Négative')
ax.set_xlabel('x'); ax.set_ylabel('Comptage y')
ax.set_title('Poisson vs Binomiale Négative\n(données sur-dispersées)')
ax.legend(fontsize=9)

# Résidus Pearson Poisson
ax = axes[1]
pearson_res = (y_od - model_pois_od.fittedvalues) / np.sqrt(model_pois_od.fittedvalues)
ax.scatter(model_pois_od.fittedvalues, pearson_res, alpha=0.4, color='steelblue', s=20)
ax.axhline(0, color='tomato', ls='--', lw=1.5)
ax.axhline(2, color='orange', ls=':', lw=1.2, label='±2')
ax.axhline(-2, color='orange', ls=':', lw=1.2)
ax.set_xlabel('Valeurs ajustées (μ̂)')
ax.set_ylabel('Résidus de Pearson')
ax.set_title('Résidus de Pearson — Modèle Poisson\n(sur-dispersion visible)')
ax.legend(fontsize=9)

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

## Régression Gamma : coûts et durées

La loi Gamma est adaptée aux variables **continues, strictement positives et asymétriques** (coûts médicaux, durées de vie, primes d'assurance). Elle peut prendre des formes exponentielles, peaked ou très asymétriques selon ses paramètres.

```{code-cell} python
# Exemple : modéliser des coûts hospitaliers
np.random.seed(2024)
n_gamma = 600

age_pat = np.random.uniform(20, 80, n_gamma)
comorbidites = np.random.poisson(1.5, n_gamma)
urgence = np.random.binomial(1, 0.3, n_gamma)

# Coût log-linéaire avec distribution Gamma
log_mu_cout = (6.5
               + 0.02 * (age_pat - 50)
               + 0.3 * comorbidites
               + 0.8 * urgence
               + np.random.normal(0, 0.1, n_gamma))
mu_cout = np.exp(log_mu_cout)

# Gamma(shape=2) : coûts avec variance = mu²/2
shape_cout = 2.0
cout = np.random.gamma(shape_cout, mu_cout / shape_cout)

df_gamma = pd.DataFrame({
    'cout': cout.round(2),
    'age': age_pat.round(1),
    'comorbidites': comorbidites,
    'urgence': urgence
})

print(f"Coût médian : {np.median(cout):.0f} €")
print(f"Coût moyen  : {cout.mean():.0f} €")
print(f"Skewness    : {stats.skew(cout):.2f}")
```

```{code-cell} python
# Modèle Gamma avec lien log
model_gamma = smf.glm(
    'cout ~ age + comorbidites + urgence',
    data=df_gamma,
    family=sm.families.Gamma(link=sm.families.links.Log())
).fit()

print(model_gamma.summary())

# Comparaison avec OLS sur log(cout)
model_ols_log = smf.ols('np.log(cout) ~ age + comorbidites + urgence',
                         data=df_gamma).fit()
print(f"\nComparaison AIC :")
print(f"  Gamma GLM : {model_gamma.aic:.1f}")
print(f"  OLS log(y): {model_ols_log.aic:.1f}")
```

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

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

# Distribution des coûts
ax = axes[0, 0]
ax.hist(cout, bins=80, density=True, color='steelblue', alpha=0.7, edgecolor='white')
ax.set_xlabel('Coût (€)'); ax.set_title('Distribution des coûts\n(asymétrie positive typique)')

# Résidus de Pearson Gamma
ax = axes[0, 1]
pearson_g = (df_gamma['cout'] - model_gamma.fittedvalues) / model_gamma.fittedvalues
ax.scatter(model_gamma.fittedvalues, pearson_g, alpha=0.3, color='steelblue', s=15)
ax.axhline(0, color='tomato', ls='--', lw=1.5)
ax.axhline(2, color='orange', ls=':', lw=1)
ax.axhline(-2, color='orange', ls=':', lw=1)
ax.set_xlabel('Valeurs ajustées (μ̂)'); ax.set_ylabel('Résidus de Pearson')
ax.set_title('Résidus de Pearson — Modèle Gamma')

# QQ-plot des résidus de deviance
ax = axes[1, 0]
dev_res_g = model_gamma.resid_deviance
qq = stats.probplot(dev_res_g, dist='norm')
ax.scatter(qq[0][0], qq[0][1], color='steelblue', s=20, alpha=0.5)
ax.plot(qq[0][0], qq[1][0] * qq[0][0] + qq[1][1], 'tomato', lw=2)
ax.set_xlabel('Quantiles théoriques'); ax.set_ylabel('Résidus de déviance')
ax.set_title('Q-Q plot des résidus de déviance\nModèle Gamma')

# Prédictions vs observés
ax = axes[1, 1]
ax.scatter(model_gamma.fittedvalues, df_gamma['cout'],
           alpha=0.3, color='steelblue', s=15)
max_val = max(model_gamma.fittedvalues.max(), df_gamma['cout'].max())
ax.plot([0, max_val], [0, max_val], 'tomato', ls='--', lw=2, label='y = ŷ')
ax.set_xlabel('Coût prédit (€)'); ax.set_ylabel('Coût observé (€)')
ax.set_title('Prédictions vs Observé — Gamma GLM')
ax.legend(fontsize=9)

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

## Déviance et AIC pour les GLM

En régression linéaire, le $R^2$ mesure la part de variance expliquée. Pour les GLM, on utilise la **déviance** :

$$D(\boldsymbol{\beta}) = 2[\ell(\boldsymbol{\hat{\beta}}_{\text{saturé}}) - \ell(\boldsymbol{\hat{\beta}})]$$

- **Déviance nulle** : modèle sans prédicteur (intercept seul)
- **Déviance résiduelle** : modèle ajusté
- **$R^2$ de McFadden** : $1 - D_{\text{résid}}/D_{\text{nulle}}$

Le critère d'information d'Akaike (**AIC**) permet la comparaison de modèles :

$$\text{AIC} = -2\ell(\hat{\boldsymbol{\beta}}) + 2p$$

```{code-cell} python
# Comparaison de modèles par AIC
models_compare = {
    'Intercepte seul': smf.glm('accidents ~ 1', data=df_pois, family=sm.families.Poisson()).fit(),
    'Ancienneté': smf.glm('accidents ~ anciennete', data=df_pois, family=sm.families.Poisson()).fit(),
    'Ancienneté + Formation': smf.glm('accidents ~ anciennete + formation', data=df_pois,
                                       family=sm.families.Poisson()).fit(),
    'Modèle complet': model_pois
}

comparaison = pd.DataFrame({
    'AIC': {k: v.aic for k, v in models_compare.items()},
    'Déviance résid.': {k: v.deviance for k, v in models_compare.items()},
    'df résid.': {k: v.df_resid for k, v in models_compare.items()}
})
print("Comparaison de modèles GLM Poisson :")
print(comparaison.round(2))
```

## Quasi-vraisemblance et correction de sur-dispersion

Quand la sur-dispersion est modérée, on peut ajuster les erreurs standard sans changer la famille via la **quasi-vraisemblance** :

```{code-cell} python
# Quasi-Poisson : même estimation, erreurs standard corrigées
model_qpois = smf.glm(
    'accidents ~ anciennete + formation + C(risque_secteur)',
    data=df_pois,
    family=sm.families.Poisson()
).fit(cov_type='HC3')  # erreurs robustes

# Comparaison des erreurs standard
se_pois = model_pois.bse
se_qpois = model_qpois.bse

print("Comparaison des erreurs standard :")
print(pd.DataFrame({
    'SE Poisson': se_pois,
    'SE Quasi-Poisson (robuste)': se_qpois,
    'Ratio': se_qpois / se_pois
}).round(4))
```

## Diagnostics GLM complets

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

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

# Données : modèle Poisson sur accidents
fitted = model_pois.fittedvalues
y_pois_obs = df_pois['accidents'].values

# 1. Résidus de Pearson vs ajustés
ax = axes[0, 0]
pears_res = (y_pois_obs - fitted) / np.sqrt(fitted)
ax.scatter(fitted, pears_res, alpha=0.4, color='steelblue', s=20)
ax.axhline(0, color='tomato', ls='--', lw=1.5)
ax.axhline(2, color='orange', ls=':', lw=1)
ax.axhline(-2, color='orange', ls=':', lw=1)
ax.set_xlabel('Valeurs ajustées'); ax.set_ylabel('Résidus de Pearson')
ax.set_title('Résidus de Pearson vs Ajustés\nPoisson GLM')

# 2. Résidus de déviance vs ajustés
ax = axes[0, 1]
dev_res = model_pois.resid_deviance
ax.scatter(fitted, dev_res, alpha=0.4, color='steelblue', s=20)
ax.axhline(0, color='tomato', ls='--', lw=1.5)
ax.axhline(2, color='orange', ls=':', lw=1)
ax.axhline(-2, color='orange', ls=':', lw=1)
ax.set_xlabel('Valeurs ajustées'); ax.set_ylabel('Résidus de déviance')
ax.set_title('Résidus de déviance vs Ajustés\nPoisson GLM')

# 3. QQ-plot résidus Poisson
ax = axes[1, 0]
qq_p = stats.probplot(dev_res, dist='norm')
ax.scatter(qq_p[0][0], qq_p[0][1], color='steelblue', s=20, alpha=0.5)
ax.plot(qq_p[0][0], qq_p[1][0] * qq_p[0][0] + qq_p[1][1], 'tomato', lw=2)
ax.set_xlabel('Quantiles théoriques (Normale)'); ax.set_ylabel('Résidus de déviance')
ax.set_title('Q-Q Plot des résidus de déviance')

# 4. Comparaison distributions prédites
ax = axes[1, 1]
vals_obs, counts_obs = np.unique(y_pois_obs, return_counts=True)
ax.bar(vals_obs - 0.2, counts_obs/len(y_pois_obs), width=0.4,
       color='steelblue', alpha=0.8, label='Observé', edgecolor='white')
# Prédit sous Poisson (simulé)
y_sim = np.random.poisson(fitted)
vals_sim, counts_sim = np.unique(y_sim, return_counts=True)
ax.bar(vals_sim + 0.2, counts_sim/len(y_sim), width=0.4,
       color='tomato', alpha=0.8, label='Simulé (modèle)', edgecolor='white')
ax.set_xlabel('Nombre d\'accidents')
ax.set_ylabel('Fréquence relative')
ax.set_title('Distributions observée vs simulée\nPosterior predictive check')
ax.legend(fontsize=9)
ax.set_xlim(-0.5, 12.5)

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

## Régression de Tweedie

La distribution de Tweedie est une famille générale paramétrisée par un indice de puissance $p$ :

| Paramètre p | Distribution |
|-------------|-------------|
| $p = 0$ | Gaussienne |
| $p = 1$ | Poisson |
| $1 < p < 2$ | Tweedie compound (mélange Poisson-Gamma) |
| $p = 2$ | Gamma |
| $p = 3$ | Gaussienne inverse |

```{code-cell} python
# Modèle Tweedie pour sinistres assurance
np.random.seed(42)
n_tw = 700

prime = np.random.lognormal(6, 0.4, n_tw)         # prime annuelle
anciennete_contrat = np.random.uniform(0, 10, n_tw)
age_conducteur = np.random.uniform(18, 75, n_tw)

log_mu_tw = (4.5
             + 0.3 * np.log(prime / 500)
             - 0.02 * anciennete_contrat
             + 0.01 * np.abs(age_conducteur - 40))
mu_tw = np.exp(log_mu_tw)

# Tweedie p=1.5 : mélange Poisson-Gamma
n_sinistres = np.random.poisson(mu_tw / 200)
sinistre_total = np.array([
    np.random.gamma(2, 100, max(ns, 0)).sum() if ns > 0 else 0.0
    for ns in n_sinistres
])

df_tw = pd.DataFrame({
    'sinistre': sinistre_total,
    'log_prime': np.log(prime),
    'anciennete': anciennete_contrat,
    'age': age_conducteur
})

print(f"Proportion zéros : {(sinistre_total == 0).mean():.1%}")
print(f"Montant médian (non-nuls) : {sinistre_total[sinistre_total > 0].mean():.0f} €")

model_tw = smf.glm(
    'sinistre ~ log_prime + anciennete + age',
    data=df_tw,
    family=sm.families.Tweedie(var_power=1.5, link=sm.families.links.Log())
).fit()

print(f"\nAIC Tweedie : {model_tw.aic:.1f}")
print(model_tw.params.round(4))
```

## Résumé comparatif

```{admonition} Choisir la bonne famille GLM
:class: note

| Réponse | Famille | Lien | Paramètre |
|---------|---------|------|-----------|
| Continu non borné | Gaussienne | Identité | σ² |
| Binaire (0/1) | Binomiale | Logit | — |
| Comptage (Var ≈ µ) | Poisson | Log | — |
| Comptage (Var >> µ) | Binomiale Négative | Log | r (dispersion) |
| Continu positif asymétrique | Gamma | Log | φ (dispersion) |
| Mélange zéros + continu | Tweedie | Log | p (puissance) |

**Règle pratique** : commencer par visualiser la distribution de la réponse, calculer ratio Var/Moyenne, puis choisir la famille. Valider avec les résidus de Pearson et un Q-Q plot.
```

```{code-cell} python
# Tableau récapitulatif avec test de sur-dispersion pour Poisson
from scipy.stats import chi2

deviance = model_pois.deviance
df_res = model_pois.df_resid
ratio_dispers = deviance / df_res

print(f"Test de sur-dispersion (Poisson) :")
print(f"  Déviance résiduelle : {deviance:.1f}")
print(f"  Degrés de liberté   : {df_res}")
print(f"  Ratio (D/df)        : {ratio_dispers:.3f}")
print(f"  p-valeur            : {chi2.sf(deviance, df_res):.4f}")
if ratio_dispers > 1.5:
    print("  → Possible sur-dispersion, envisager Binomiale Négative ou Quasi-Poisson")
else:
    print("  → Pas de sur-dispersion significative, Poisson convenable")
```