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

# Tests paramétriques

Les tests paramétriques constituent l'arsenal classique de l'inférence statistique : ils permettent de répondre à des questions concrètes — « ce traitement est-il efficace ? », « ces deux machines produisent-elles des pièces de même calibre ? », « les trois méthodes d'enseignement donnent-elles des résultats équivalents ? » — avec une rigueur formelle. Ce chapitre présente les tests les plus utilisés, leurs hypothèses, leur implémentation en Python et leur interprétation correcte.

```{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.anova import AnovaRM
from statsmodels.stats import power as smp
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')
```

## Rappels : cadre des tests d'hypothèses

Un test statistique repose sur une logique de réfutation : on formule une hypothèse nulle $H_0$ (statu quo, absence d'effet), on calcule la probabilité d'observer des données au moins aussi extrêmes si $H_0$ était vraie, et on décide de rejeter ou non $H_0$.

**Éléments fondamentaux :**

- **$H_0$ (hypothèse nulle)** : hypothèse par défaut, que l'on cherche à réfuter. Exemple : $\mu = \mu_0$, ou $\mu_1 = \mu_2$.
- **$H_1$ (hypothèse alternative)** : ce que l'on cherche à établir. Peut être bilatérale ($\mu \neq \mu_0$) ou unilatérale ($\mu > \mu_0$).
- **p-valeur** : $P(\text{statistique} \geq t_{\text{obs}} \mid H_0)$. Ce n'est *pas* la probabilité que $H_0$ soit vraie — c'est la probabilité d'obtenir ce résultat ou plus extrême *si $H_0$ était vraie*.
- **Niveau $\alpha$** : seuil de rejet fixé *avant* le test (souvent 0,05). Si $p < \alpha$, on rejette $H_0$.

**Les deux types d'erreurs :**

| Décision \ Réalité | $H_0$ vraie | $H_0$ fausse |
|---|---|---|
| Ne pas rejeter $H_0$ | Décision correcte | **Erreur de type II** ($\beta$) |
| Rejeter $H_0$ | **Erreur de type I** ($\alpha$) | Décision correcte |

La **puissance** du test est $1 - \beta$ : la probabilité de détecter un effet réel.

```{admonition} p-valeur : ce qu'elle n'est pas
:class: warning

La p-valeur n'est *pas* :
- La probabilité que $H_0$ soit vraie
- La probabilité que les résultats soient dus au hasard
- Une mesure de l'importance pratique de l'effet

Une p-valeur < 0,05 signifie seulement que les données sont incompatibles avec $H_0$ au seuil $\alpha = 5\%$. La *taille d'effet* mesure l'importance pratique.
```

## Test t à un échantillon

### Contexte et hypothèses

On dispose d'un échantillon $x_1, \ldots, x_n$ et on veut tester si la moyenne de la population dont il est issu est égale à une valeur de référence $\mu_0$.

$$H_0 : \mu = \mu_0 \qquad H_1 : \mu \neq \mu_0$$

La statistique de test est :
$$t = \frac{\bar{x} - \mu_0}{s / \sqrt{n}}$$

Sous $H_0$, $t \sim \mathcal{T}(n-1)$ (loi de Student à $n-1$ degrés de liberté), à condition que les observations soient indépendantes et que la population soit normale (ou $n$ suffisamment grand pour le TCL).

### Exemple : calibrage d'une machine

Une machine est censée produire des pièces de 50 mm. On mesure 30 pièces. La machine est-elle bien calibrée ?

```{code-cell} python3
# Simulation : machine légèrement déréglée
mu_reel = 50.8  # vrai écart
sigma = 2.0
n = 30
mesures = rng.normal(mu_reel, sigma, n)

# Test t à un échantillon
t_stat, p_val = stats.ttest_1samp(mesures, popmean=50.0)

print(f"Moyenne observée : {mesures.mean():.3f} mm")
print(f"Statistique t    : {t_stat:.3f}")
print(f"p-valeur         : {p_val:.4f}")
print(f"Degrés de liberté: {n - 1}")
print()
if p_val < 0.05:
    print("→ Rejet de H0 : la machine est probablement déréglée.")
else:
    print("→ Non-rejet de H0 : pas de preuve de dérèglage.")
```

### Taille d'effet : d de Cohen

La p-valeur seule est insuffisante. Le **d de Cohen** mesure l'ampleur de l'écart en unités d'écarts-types :

$$d = \frac{\bar{x} - \mu_0}{s}$$

| Valeur absolue de $d$ | Interprétation |
|---|---|
| 0,2 | Petit effet |
| 0,5 | Effet modéré |
| 0,8 | Grand effet |

```{code-cell} python3
# Taille d'effet
cohens_d = (mesures.mean() - 50.0) / mesures.std(ddof=1)
print(f"d de Cohen : {cohens_d:.3f}")

# Avec pingouin (résultat complet)
result = pg.ttest(mesures, 50.0)
print("\nRésultat complet avec pingouin :")
print(result[['T', 'dof', 'alternative', 'p_val', 'CI95', 'cohen_d', 'power']].to_string())
```

### Visualisation : distribution sous H₀ et statistique observée

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

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

# Axe gauche : données brutes
ax = axes[0]
ax.hist(mesures, bins=10, color='steelblue', alpha=0.7, edgecolor='white')
ax.axvline(50.0, color='crimson', lw=2, linestyle='--', label='Référence μ₀ = 50')
ax.axvline(mesures.mean(), color='navy', lw=2, label=f'Moyenne obs. = {mesures.mean():.2f}')
ax.set_xlabel('Longueur (mm)')
ax.set_ylabel('Effectif')
ax.set_title('Distribution des mesures')
ax.legend()

# Axe droit : distribution de Student sous H0 avec région de rejet
ax = axes[1]
df_val = n - 1
x = np.linspace(-4, 4, 400)
y = stats.t.pdf(x, df=df_val)
ax.plot(x, y, 'k-', lw=2, label=f'Student({df_val})')

# Régions de rejet bilatérales (alpha = 0.05)
t_crit = stats.t.ppf(0.975, df=df_val)
mask_left = x < -t_crit
mask_right = x > t_crit
ax.fill_between(x[mask_left], y[mask_left], alpha=0.4, color='crimson', label='Région de rejet (5%)')
ax.fill_between(x[mask_right], y[mask_right], alpha=0.4, color='crimson')

# Statistique observée
ax.axvline(t_stat, color='navy', lw=2.5, linestyle='--', label=f't obs. = {t_stat:.2f}')
ax.set_xlabel('t')
ax.set_ylabel('Densité')
ax.set_title('Distribution sous H₀ et statistique de test')
ax.legend(fontsize=9)

plt.tight_layout()
plt.savefig('_static/06_ttest1.png', bbox_inches='tight')
plt.show()
```

## Test t à deux échantillons indépendants

### Student vs Welch

Deux groupes indépendants $x_1, \ldots, x_{n_1}$ et $y_1, \ldots, y_{n_2}$. On teste $H_0 : \mu_X = \mu_Y$.

**Test de Student** (variantes égales) : $t = \frac{\bar{x} - \bar{y}}{s_p \sqrt{1/n_1 + 1/n_2}}$ avec $s_p^2$ variance poolée.

**Test de Welch** (variances inégales, recommandé par défaut) : ajuste les degrés de liberté selon l'approximation de Welch-Satterthwaite.

```{admonition} Quelle variante utiliser ?
:class: tip

En pratique, utilisez **toujours le test de Welch** (`equal_var=False` dans scipy). Il est robuste même lorsque les variances sont égales, alors que le test de Student peut être anti-conservateur lorsque les variances diffèrent. Vérifier l'égalité des variances avec un test de Levene et adapter en conséquence est une bonne pratique.
```

### Exemple : deux méthodes pédagogiques

```{code-cell} python3
# Deux groupes d'élèves (scores sur 100)
n_A, n_B = 35, 40
groupe_A = rng.normal(68, 12, n_A)  # méthode traditionnelle
groupe_B = rng.normal(74, 15, n_B)  # méthode active

# Test de Levene (égalité des variances)
lev_stat, lev_p = stats.levene(groupe_A, groupe_B)
print(f"Test de Levene : F = {lev_stat:.3f}, p = {lev_p:.4f}")
print(f"→ Variances {'probablement égales' if lev_p > 0.05 else 'probablement inégales'}")
print()

# Test de Welch (recommandé)
t_welch, p_welch = stats.ttest_ind(groupe_A, groupe_B, equal_var=False)
# Test de Student (pour comparaison)
t_student, p_student = stats.ttest_ind(groupe_A, groupe_B, equal_var=True)

print(f"Test de Welch   : t = {t_welch:.3f}, p = {p_welch:.4f}")
print(f"Test de Student : t = {t_student:.3f}, p = {p_student:.4f}")
print()

# Avec pingouin
res = pg.ttest(groupe_A, groupe_B, correction=True)
print("Résultat pingouin (Welch) :")
print(res[['T', 'dof', 'p_val', 'CI95', 'cohen_d', 'power']].to_string())
```

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

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

# Boxplot + swarmplot
ax = axes[0]
df_scores = pd.DataFrame({
    'Score': np.concatenate([groupe_A, groupe_B]),
    'Méthode': ['Traditionnelle'] * n_A + ['Active'] * n_B
})
sns.boxplot(data=df_scores, x='Méthode', y='Score', ax=ax,
            palette=['#4C72B0', '#DD8452'], width=0.5)
sns.stripplot(data=df_scores, x='Méthode', y='Score', ax=ax,
              color='black', alpha=0.4, size=4, jitter=True)
ax.set_title('Distribution des scores par méthode')
ax.set_ylabel('Score (/100)')

# Intervalles de confiance sur les moyennes
ax = axes[1]
for i, (groupe, label, color) in enumerate(zip(
    [groupe_A, groupe_B], ['Traditionnelle', 'Active'],
    ['#4C72B0', '#DD8452']
)):
    mean = groupe.mean()
    sem = stats.sem(groupe)
    ci = stats.t.interval(0.95, len(groupe)-1, loc=mean, scale=sem)
    ax.errorbar(i, mean, yerr=[[mean - ci[0]], [ci[1] - mean]],
                fmt='o', color=color, capsize=8, capthick=2,
                markersize=10, linewidth=2, label=label)

ax.set_xticks([0, 1])
ax.set_xticklabels(['Traditionnelle', 'Active'])
ax.set_ylabel('Score moyen ± IC 95%')
ax.set_title('Comparaison des moyennes avec IC 95%')
ax.set_xlim(-0.5, 1.5)
ax.legend()

plt.tight_layout()
plt.savefig('_static/06_ttest2.png', bbox_inches='tight')
plt.show()
```

## Test t apparié

### Données appariées : avant/après

Lorsque les deux mesures sont liées (même individu mesuré deux fois, jumeaux, paires matchées), les observations ne sont pas indépendantes. Le test t apparié exploite cette structure en travaillant sur les différences $d_i = x_{i,\text{après}} - x_{i,\text{avant}}$.

$$H_0 : \mu_d = 0 \qquad t = \frac{\bar{d}}{s_d / \sqrt{n}}$$

Le test t apparié est *beaucoup plus puissant* que le test à deux échantillons indépendants quand la corrélation intra-paire est élevée, car il élimine la variabilité inter-individus.

```{code-cell} python3
# Exemple médical : pression artérielle avant/après traitement
n_patients = 25
avant = rng.normal(145, 15, n_patients)  # mmHg, avant traitement
# Effet du traitement : réduction de ~10 mmHg + variabilité individuelle
effet_traitement = rng.normal(10, 5, n_patients)
apres = avant - effet_traitement

differences = apres - avant

# Test t apparié
t_app, p_app = stats.ttest_rel(avant, apres)

# Test t indépendant (incorrect ici, pour illustration)
t_ind, p_ind = stats.ttest_ind(avant, apres)

print("Test t apparié (correct) :")
print(f"  t = {t_app:.3f}, p = {p_app:.6f}")
print()
print("Test t indépendant (incorrect pour ces données) :")
print(f"  t = {t_ind:.3f}, p = {p_ind:.4f}")
print()
print(f"Différence moyenne : {differences.mean():.2f} mmHg")
print(f"IC 95% sur la différence : {np.percentile(differences, [2.5, 97.5])}")
print()

res_app = pg.ttest(avant, apres, paired=True)
print("Résultat pingouin (apparié) :")
print(res_app[['T', 'dof', 'p_val', 'CI95', 'cohen_d', 'power']].to_string())
```

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

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

# Paires avant/après
ax = axes[0]
for i in range(n_patients):
    ax.plot([0, 1], [avant[i], apres[i]], 'o-', color='steelblue', alpha=0.4, linewidth=1)
ax.plot([0, 1], [avant.mean(), apres.mean()], 'o-', color='crimson',
        linewidth=3, markersize=10, label='Moyenne')
ax.set_xticks([0, 1])
ax.set_xticklabels(['Avant', 'Après'])
ax.set_ylabel('Pression artérielle (mmHg)')
ax.set_title('Pression artérielle par patient')
ax.legend()

# Distribution des différences
ax = axes[1]
ax.hist(differences, bins=10, color='steelblue', alpha=0.7, edgecolor='white')
ax.axvline(0, color='crimson', lw=2, linestyle='--', label='H₀ : diff = 0')
ax.axvline(differences.mean(), color='navy', lw=2, label=f'Moyenne = {differences.mean():.1f}')
ax.set_xlabel('Différence (Après − Avant) en mmHg')
ax.set_ylabel('Effectif')
ax.set_title('Distribution des différences')
ax.legend()

# Comparaison puissance : apparié vs indépendant
ax = axes[2]
n_range = np.arange(5, 50)
# Approximation : effet moyen / std différences vs effet moyen / std poolée
d_apparie = abs(differences.mean()) / differences.std(ddof=1)
sigma_pool = np.sqrt((avant.var(ddof=1) + apres.var(ddof=1)) / 2)
d_indep = abs(differences.mean()) / sigma_pool

pow_app = [smp.TTestPower().power(effect_size=d_apparie, nobs=n, alpha=0.05) for n in n_range]
pow_ind = [smp.TTestIndPower().power(effect_size=d_indep, nobs1=n, alpha=0.05, ratio=1) for n in n_range]

ax.plot(n_range, pow_app, 'steelblue', lw=2, label='T apparié')
ax.plot(n_range, pow_ind, 'crimson', lw=2, linestyle='--', label='T indépendant')
ax.axhline(0.8, color='gray', linestyle=':', lw=1.5, label='Puissance 80%')
ax.set_xlabel('Taille d\'échantillon n')
ax.set_ylabel('Puissance')
ax.set_title('Puissance : apparié vs indépendant')
ax.legend(fontsize=9)
ax.set_ylim(0, 1.05)

plt.tight_layout()
plt.savefig('_static/06_apparie.png', bbox_inches='tight')
plt.show()
```

## ANOVA à un facteur

### Principe

L'ANOVA (ANalysis Of VAriance) étend le test t à plus de deux groupes. On décompose la variance totale en variance *inter-groupes* (due aux traitements) et variance *intra-groupes* (résiduelle).

$$F = \frac{SS_{\text{inter}} / (k-1)}{SS_{\text{intra}} / (N-k)} = \frac{MS_{\text{entre}}}{MS_{\text{intra}}}$$

où $k$ est le nombre de groupes et $N$ le nombre total d'observations.

**Hypothèses de l'ANOVA :**

1. Les observations sont indépendantes
2. Les résidus sont normalement distribués
3. Les variances sont homogènes (homoscédasticité) — vérifiée par Levene ou Bartlett

```{admonition} Pourquoi ne pas faire plusieurs tests t ?
:class: warning

Avec $k$ groupes, faire tous les tests t par paires donne $\binom{k}{2}$ tests. Au seuil $\alpha = 0,05$, la probabilité d'au moins une fausse alarme est $1 - (1-\alpha)^m$ où $m$ est le nombre de tests. Pour $k = 5$ groupes (10 tests), ce taux atteint ~40%. L'ANOVA contrôle globalement $\alpha$.
```

### Exemple : trois régimes alimentaires

```{code-cell} python3
# Trois régimes, effet sur la glycémie à jeun (mg/dL)
regime_A = rng.normal(95, 10, 30)   # régime standard
regime_B = rng.normal(88, 10, 30)   # régime faible glucides
regime_C = rng.normal(92, 10, 30)   # régime méditerranéen

# ANOVA avec scipy
f_stat, p_anova = stats.f_oneway(regime_A, regime_B, regime_C)
print(f"ANOVA à un facteur :")
print(f"  F = {f_stat:.3f}")
print(f"  p = {p_anova:.4f}")
print()

# Table ANOVA complète avec pingouin
df_glyc = pd.DataFrame({
    'glycemie': np.concatenate([regime_A, regime_B, regime_C]),
    'regime': ['Standard'] * 30 + ['Faible glucides'] * 30 + ['Méditerranéen'] * 30
})

aov = pg.anova(data=df_glyc, dv='glycemie', between='regime', detailed=True)
print("Table ANOVA (pingouin) :")
print(aov.to_string(index=False))
print()
print(f"η² (eta-carré, taille d'effet) : {aov.loc[0, 'np2']:.3f}")
```

### Tests post-hoc : Tukey HSD et Bonferroni

L'ANOVA indique si *au moins une* paire de moyennes diffère — les tests post-hoc précisent lesquelles.

```{code-cell} python3
# Tests post-hoc avec pingouin
print("Post-hoc Tukey HSD :")
tukey = pg.pairwise_tukey(data=df_glyc, dv='glycemie', between='regime')
print(tukey[['A', 'B', 'mean_A', 'mean_B', 'diff', 'se', 'T', 'p_tukey']].to_string(index=False))
print()

print("Post-hoc Bonferroni :")
bonf = pg.pairwise_tests(data=df_glyc, dv='glycemie', between='regime', padjust='bonf')
print(bonf[['A', 'B', 'T', 'dof', 'p_unc', 'p_corr']].to_string(index=False))
```

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

fig, axes = plt.subplots(1, 3, figsize=(16, 4.5))

# Boxplot comparatif
ax = axes[0]
palette = {'Standard': '#4C72B0', 'Faible glucides': '#55A868', 'Méditerranéen': '#C44E52'}
sns.boxplot(data=df_glyc, x='regime', y='glycemie', ax=ax, palette=palette)
sns.stripplot(data=df_glyc, x='regime', y='glycemie', ax=ax,
              color='black', alpha=0.3, size=3)
ax.set_title('Glycémie par régime')
ax.set_xlabel('Régime')
ax.set_ylabel('Glycémie à jeun (mg/dL)')
ax.tick_params(axis='x', rotation=10)

# Décomposition de la variance (barres)
ax = axes[1]
ss_inter = aov.loc[0, 'SS']
ss_intra = aov.loc[1, 'SS']
ss_total = ss_inter + ss_intra
bars = ax.bar(['SS inter\n(régimes)', 'SS intra\n(résidus)'],
              [ss_inter, ss_intra],
              color=['#C44E52', '#4C72B0'], alpha=0.8)
ax.bar_label(bars, fmt='%.1f', padding=3)
ax.set_ylabel('Somme des carrés')
ax.set_title(f'Décomposition de la variance\n(η² = {ss_inter/ss_total:.2f})')

# Distribution F sous H0 avec statistique observée
ax = axes[2]
df1, df2 = 2, len(df_glyc) - 3
x_f = np.linspace(0, 8, 400)
y_f = stats.f.pdf(x_f, df1, df2)
ax.plot(x_f, y_f, 'k-', lw=2, label=f'F({df1}, {df2})')
f_crit = stats.f.ppf(0.95, df1, df2)
mask_rejet = x_f > f_crit
ax.fill_between(x_f[mask_rejet], y_f[mask_rejet], alpha=0.4, color='crimson', label='Région de rejet (5%)')
ax.axvline(f_stat, color='navy', lw=2.5, linestyle='--', label=f'F obs. = {f_stat:.2f}')
ax.set_xlabel('F')
ax.set_ylabel('Densité')
ax.set_title('Distribution de Fisher sous H₀')
ax.legend(fontsize=9)

plt.tight_layout()
plt.savefig('_static/06_anova1.png', bbox_inches='tight')
plt.show()
```

## ANOVA à deux facteurs

L'ANOVA à deux facteurs analyse simultanément deux variables catégorielles et leur éventuelle **interaction**.

### Exemple : médicament × dose

```{code-cell} python3
# Plan factoriel 2×3 : 2 médicaments × 3 doses
rng2 = np.random.default_rng(123)
n_par_cellule = 15

medicaments = ['Méd. A', 'Méd. B']
doses = ['Faible', 'Moyenne', 'Forte']

data_rows = []
for med in medicaments:
    for dose in doses:
        # Effets principaux + interaction
        base = 80
        effet_med = 5 if med == 'Méd. B' else 0
        effet_dose = {'Faible': 0, 'Moyenne': 8, 'Forte': 14}[dose]
        # Interaction : Méd B est plus efficace à forte dose
        interaction = 6 if (med == 'Méd. B' and dose == 'Forte') else 0
        mu = base - effet_med - effet_dose - interaction
        vals = rng2.normal(mu, 8, n_par_cellule)
        for v in vals:
            data_rows.append({'medicament': med, 'dose': dose, 'reponse': v})

df_fact = pd.DataFrame(data_rows)

# ANOVA à deux facteurs avec pingouin
aov2 = pg.anova(data=df_fact, dv='reponse', between=['medicament', 'dose'], detailed=True)
print("ANOVA à deux facteurs :")
print(aov2[['Source', 'SS', 'DF', 'MS', 'F', 'p_unc', 'np2']].to_string(index=False))
```

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

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

# Interaction plot
ax = axes[0]
means = df_fact.groupby(['medicament', 'dose'])['reponse'].mean().reset_index()
dose_order = ['Faible', 'Moyenne', 'Forte']
for med, color in zip(medicaments, ['#4C72B0', '#DD8452']):
    sub = means[means['medicament'] == med].copy()
    sub['dose'] = pd.Categorical(sub['dose'], categories=dose_order, ordered=True)
    sub = sub.sort_values('dose')
    ax.plot(sub['dose'], sub['reponse'], 'o-', color=color, lw=2, markersize=8, label=med)
ax.set_xlabel('Dose')
ax.set_ylabel('Réponse moyenne')
ax.set_title('Graphique d\'interaction\n(croisement = interaction significative)')
ax.legend()

# Boxplot groupé
ax = axes[1]
sns.boxplot(data=df_fact, x='dose', y='reponse', hue='medicament',
            order=dose_order, ax=ax, palette=['#4C72B0', '#DD8452'])
ax.set_xlabel('Dose')
ax.set_ylabel('Réponse')
ax.set_title('Distribution par dose et médicament')
ax.legend(title='Médicament')

plt.tight_layout()
plt.savefig('_static/06_anova2.png', bbox_inches='tight')
plt.show()
```

## Test de Fisher : comparer deux variances

Le test F de Fisher teste $H_0 : \sigma_1^2 = \sigma_2^2$ via la statistique $F = s_1^2 / s_2^2$.

```{code-cell} python3
# Exemple : variabilité de deux procédés industriels
proc_A = rng.normal(100, 5, 30)
proc_B = rng.normal(100, 8, 30)

var_A = proc_A.var(ddof=1)
var_B = proc_B.var(ddof=1)
F_obs = var_A / var_B

# p-valeur bilatérale
df1, df2 = len(proc_A) - 1, len(proc_B) - 1
p_droite = stats.f.sf(F_obs, df1, df2)
p_gauche = stats.f.cdf(F_obs, df1, df2)
p_fisher = 2 * min(p_droite, p_gauche)

print(f"Variance procédé A : {var_A:.2f}")
print(f"Variance procédé B : {var_B:.2f}")
print(f"F = s²A / s²B = {F_obs:.3f}")
print(f"p-valeur (bilatéral) = {p_fisher:.4f}")
print()
print("Note : Levene est plus robuste que Fisher pour tester l'égalité des variances.")
lev_stat, lev_p = stats.levene(proc_A, proc_B)
print(f"Test de Levene : F = {lev_stat:.3f}, p = {lev_p:.4f}")
```

## Vérification des hypothèses

### Normalité : test de Shapiro-Wilk

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

fig, axes = plt.subplots(2, 4, figsize=(16, 7))

# Générer différentes distributions pour illustration
distributions = {
    'Normal': rng.normal(0, 1, 60),
    'Légèrement asymétrique': rng.exponential(1, 60) - 1,
    'Fortement asymétrique': rng.exponential(0.3, 60),
    'Bimodale': np.concatenate([rng.normal(-2, 0.8, 30), rng.normal(2, 0.8, 30)]),
}

for i, (name, data) in enumerate(distributions.items()):
    # Histogramme + courbe normale ajustée
    ax = axes[0, i]
    ax.hist(data, bins=15, density=True, alpha=0.6, color='steelblue', edgecolor='white')
    x_range = np.linspace(data.min(), data.max(), 200)
    ax.plot(x_range, stats.norm.pdf(x_range, data.mean(), data.std()), 'crimson', lw=2)
    stat, p = stats.shapiro(data)
    ax.set_title(f'{name}\nShapiro p = {p:.3f}', fontsize=9)
    ax.set_xlabel('Valeur')

    # QQ-plot
    ax = axes[1, i]
    (osm, osr), (slope, intercept, r) = stats.probplot(data, dist='norm')
    ax.plot(osm, osr, 'o', color='steelblue', alpha=0.6, markersize=4)
    ax.plot(osm, slope * np.array(osm) + intercept, 'crimson', lw=2)
    ax.set_title(f'QQ-plot', fontsize=9)
    ax.set_xlabel('Quantiles théoriques')
    ax.set_ylabel('Quantiles observés')

axes[0, 0].set_ylabel('Densité')
axes[1, 0].set_ylabel('Quantiles observés')
plt.suptitle('Vérification de normalité : histogrammes et QQ-plots', y=1.01, fontsize=12)
plt.tight_layout()
plt.savefig('_static/06_normalite.png', bbox_inches='tight')
plt.show()
```

### Homoscédasticité : Levene et Bartlett

```{code-cell} python3
# Trois groupes, dont un avec variance différente
g1 = rng.normal(50, 5, 30)
g2 = rng.normal(52, 5, 30)
g3 = rng.normal(54, 15, 30)  # variance beaucoup plus grande

lev_stat, lev_p = stats.levene(g1, g2, g3)
bart_stat, bart_p = stats.bartlett(g1, g2, g3)

print("Test de Levene (robuste aux non-normalités) :")
print(f"  F = {lev_stat:.3f}, p = {lev_p:.4f}")
print()
print("Test de Bartlett (sensible aux non-normalités) :")
print(f"  χ² = {bart_stat:.3f}, p = {bart_p:.4f}")
print()
print(f"Variance G1: {g1.var(ddof=1):.1f}")
print(f"Variance G2: {g2.var(ddof=1):.1f}")
print(f"Variance G3: {g3.var(ddof=1):.1f}")
```

```{admonition} Levene vs Bartlett
:class: note

Le test de **Levene** est préféré en pratique : il est robuste aux écarts à la normalité, contrairement au test de **Bartlett** qui suppose que les groupes sont normalement distribués. Si Levene est significatif, envisagez le test de Welch (deux groupes) ou le test de Brown-Forsythe (plusieurs groupes).
```

## Puissance statistique

La puissance $1 - \beta$ est la probabilité de détecter un effet réel. Elle dépend de :

- La taille de l'effet ($d$ de Cohen, $f$ de Cohen, $\eta^2$...)
- La taille de l'échantillon $n$
- Le seuil $\alpha$

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

# Courbes de puissance pour le test t (deux échantillons indépendants)
fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))

# Puissance vs n pour différentes tailles d'effet
ax = axes[0]
n_range = np.arange(5, 150)
for d, label, color in zip([0.2, 0.5, 0.8, 1.2],
                            ['Petit (d=0.2)', 'Modéré (d=0.5)', 'Grand (d=0.8)', 'Très grand (d=1.2)'],
                            ['#e41a1c', '#ff7f00', '#4daf4a', '#377eb8']):
    power_vals = [smp.TTestIndPower().power(effect_size=d, nobs1=n, alpha=0.05, ratio=1)
                  for n in n_range]
    ax.plot(n_range, power_vals, color=color, lw=2, label=label)

ax.axhline(0.80, color='black', linestyle=':', lw=1.5, label='Puissance 80%')
ax.axhline(0.90, color='gray', linestyle=':', lw=1.5, label='Puissance 90%')
ax.set_xlabel('Taille d\'échantillon (par groupe)')
ax.set_ylabel('Puissance (1 - β)')
ax.set_title('Courbes de puissance — Test t indépendant')
ax.legend(fontsize=8)
ax.set_ylim(0, 1.05)

# Calcul du n nécessaire
ax = axes[1]
alphas = [0.01, 0.05, 0.10]
d_range = np.linspace(0.2, 1.5, 100)
for alpha, color in zip(alphas, ['#e41a1c', '#4daf4a', '#377eb8']):
    n_needed = [smp.TTestIndPower().solve_power(effect_size=d, alpha=alpha, power=0.80, ratio=1)
                for d in d_range]
    ax.plot(d_range, n_needed, color=color, lw=2, label=f'α = {alpha}')

ax.set_xlabel('Taille d\'effet (d de Cohen)')
ax.set_ylabel('N nécessaire par groupe')
ax.set_title('Taille d\'échantillon nécessaire (puissance = 80%)')
ax.legend()
ax.set_ylim(0, 400)
ax.axvline(0.5, color='gray', linestyle=':', lw=1)
ax.text(0.5, 380, 'Effet modéré', ha='center', fontsize=9, color='gray')

plt.tight_layout()
plt.savefig('_static/06_puissance.png', bbox_inches='tight')
plt.show()
```

```{code-cell} python3
# Calcul de puissance et de n nécessaire
print("=== Calcul de puissance ===")
analysis = smp.TTestIndPower()

# Puissance pour n donné
n_obs = 30
d_cible = 0.5
pow_calc = analysis.power(effect_size=d_cible, nobs1=n_obs, alpha=0.05, ratio=1)
print(f"Puissance avec n={n_obs} par groupe, d={d_cible} : {pow_calc:.1%}")

# N nécessaire pour puissance cible
pow_cible = 0.80
n_nec = analysis.solve_power(effect_size=d_cible, alpha=0.05, power=pow_cible, ratio=1)
print(f"\nN nécessaire pour 80% de puissance, d={d_cible} : {n_nec:.0f} par groupe")

# ANOVA : puissance
print("\n=== ANOVA — puissance ===")
# f de Cohen pour ANOVA (f = 0.10 petit, 0.25 modéré, 0.40 grand)
f_cohen = 0.25
k_groupes = 3
anova_power = smp.FTestAnovaPower()
n_anova = anova_power.solve_power(effect_size=f_cohen, alpha=0.05, power=0.80, k_groups=k_groupes)
print(f"N par groupe pour ANOVA ({k_groupes} groupes), f={f_cohen}, puissance=80% : {n_anova:.0f}")
```

## Résumé : choisir le bon test

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

# Tableau récapitulatif
recap = pd.DataFrame({
    'Situation': [
        '1 groupe, μ connu',
        '2 groupes indépendants',
        '2 groupes indépendants (var ≠)',
        '2 groupes appariés',
        'k groupes indépendants',
        'k groupes, 2 facteurs',
        '2 variances',
    ],
    'Test': [
        't à 1 échantillon',
        't de Student',
        't de Welch (recommandé)',
        't apparié',
        'ANOVA à 1 facteur',
        'ANOVA à 2 facteurs',
        'Test F de Fisher',
    ],
    'Fonction scipy': [
        'ttest_1samp()',
        'ttest_ind(equal_var=True)',
        'ttest_ind(equal_var=False)',
        'ttest_rel()',
        'f_oneway()',
        'smf.ols() + anova_lm()',
        'f.ppf() / levene()',
    ],
    'Taille d\'effet': [
        'd de Cohen',
        'd de Cohen',
        'd de Cohen',
        'd de Cohen',
        'η² (eta-carré)',
        'η² partiel',
        'R² de variance',
    ]
})

print(recap.to_string(index=False))
```

```{admonition} Bonnes pratiques des tests paramétriques
:class: tip

1. **Vérifier les hypothèses** avant le test : normalité (Shapiro-Wilk, QQ-plot), homoscédasticité (Levene).
2. **Rapporter la taille d'effet** (d de Cohen, η²) en plus de la p-valeur.
3. **Calculer la puissance** a priori pour déterminer la taille d'échantillon, et a posteriori pour interpréter les résultats non significatifs.
4. **Utiliser Welch par défaut** pour les comparaisons à deux groupes — il est robuste.
5. **Ne pas confondre** significativité statistique et importance pratique : un effet très faible peut être significatif avec un grand $n$.
6. **Fixer $\alpha$ avant** de voir les données.
```