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

# Comparaisons multiples et reproductibilité

Tester une seule hypothèse à 5% de risque d'erreur, c'est raisonnable. Mais que se passe-t-il quand on effectue 20, 100, ou 10 000 tests simultanément — comme en génomique, en neuro-imagerie, ou lors d'analyses exploratoires de données ? La probabilité d'obtenir au moins un faux positif explose. Ce chapitre présente les méthodes de correction pour comparaisons multiples, illustre le problème du *p-hacking*, et aborde la crise de la reproductibilité.

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

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import seaborn as sns
from scipy import stats
from statsmodels.stats.multitest import multipletests
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')
```

## Le problème des comparaisons multiples

### Inflation du taux d'erreur de type I

Si on effectue $m$ tests indépendants, chacun au seuil $\alpha$, la probabilité d'obtenir *au moins un* faux positif est :

$$P(\text{au moins un faux positif}) = 1 - (1 - \alpha)^m$$

Pour $\alpha = 0{,}05$ et $m = 20$ tests : $1 - 0{,}95^{20} \approx 64\%$.

```{code-cell} python3
# Démonstration par simulation : taux de faux positifs réels
def simulation_faux_positifs(m_tests, alpha=0.05, n_simul=500, n_obs=30):
    """Simule m_tests sous H0 vraie et compte les faux positifs (vectorisé)."""
    # Génère toutes les données en une seule fois : (n_simul, m_tests, 2, n_obs)
    data = rng.normal(0, 1, (n_simul, m_tests, 2, n_obs))
    g1, g2 = data[:, :, 0, :], data[:, :, 1, :]
    # Welch t-test vectorisé
    m1, m2 = g1.mean(axis=2), g2.mean(axis=2)
    s1, s2 = g1.var(axis=2, ddof=1), g2.var(axis=2, ddof=1)
    se = np.sqrt(s1 / n_obs + s2 / n_obs)
    t = (m1 - m2) / se
    df = (s1 / n_obs + s2 / n_obs)**2 / (
        (s1 / n_obs)**2 / (n_obs - 1) + (s2 / n_obs)**2 / (n_obs - 1)
    )
    p_valeurs = 2 * stats.t.sf(np.abs(t), df=df)   # (n_simul, m_tests)
    return (p_valeurs < alpha).sum(axis=1)           # nb faux positifs par simulation

m_values = [1, 5, 10, 20, 50, 100]
taux_fwer = []
for m in m_values:
    fp = simulation_faux_positifs(m, n_simul=500)
    taux_fwer.append((fp >= 1).mean())

# Valeurs théoriques
taux_theorique = [1 - (1 - 0.05)**m for m in m_values]

print("Nombre de tests | Taux FWER simulé | Taux FWER théorique")
print("-" * 55)
for m, sim, theo in zip(m_values, taux_fwer, taux_theorique):
    print(f"{m:>15} | {sim:>16.3f} | {theo:>19.3f}")
```

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

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

# Courbe FWER vs nombre de tests
ax = axes[0]
m_range = np.arange(1, 101)
taux_th = 1 - (1 - 0.05)**m_range
ax.plot(m_range, taux_th, 'crimson', lw=2.5, label='α = 0.05')
ax.plot(m_range, 1 - (1 - 0.01)**m_range, 'steelblue', lw=2, linestyle='--', label='α = 0.01')
ax.scatter(m_values, taux_fwer, color='black', zorder=5, s=60, label='Simulation')
ax.axhline(0.05, color='gray', linestyle=':', lw=1.5)
ax.text(85, 0.07, 'Nominal 5%', color='gray', fontsize=9)
ax.set_xlabel('Nombre de tests m')
ax.set_ylabel('P(au moins un faux positif)')
ax.set_title('Inflation du taux d\'erreur\n(FWER — Family-Wise Error Rate)')
ax.legend()
ax.set_ylim(0, 1.05)

# Distribution des p-valeurs sous H0
ax = axes[1]
n_tests_demo = 200
p_vals_h0 = []
for _ in range(n_tests_demo):
    g1 = rng.normal(0, 1, 30)
    g2 = rng.normal(0, 1, 30)
    _, p = stats.ttest_ind(g1, g2)
    p_vals_h0.append(p)

ax.hist(p_vals_h0, bins=20, density=True, color='steelblue', alpha=0.7, edgecolor='white')
ax.axhline(1.0, color='crimson', lw=2, linestyle='--', label='Uniforme(0,1) sous H₀')
ax.axvline(0.05, color='orange', lw=2, linestyle='--', label=f'α = 0.05 ({sum(p<0.05 for p in p_vals_h0)} faux pos.)')
ax.set_xlabel('p-valeur')
ax.set_ylabel('Densité')
ax.set_title(f'Distribution des p-valeurs sous H₀\n({n_tests_demo} tests, aucun effet réel)')
ax.legend(fontsize=9)

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

## Le p-hacking : démonstration

Le *p-hacking* consiste à chercher des analyses ou des sous-groupes jusqu'à obtenir $p < 0{,}05$, puis à ne rapporter que ce résultat. C'est une forme de biais de publication et une menace majeure pour la reproductibilité.

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

# Simulation du p-hacking : on collecte des données progressivement
# et on s'arrête dès que p < 0.05

def _welch_p(g1, g2):
    """Welch t-test vectorisé sur les n premiers éléments cumulés."""
    # g1, g2 : (n_simul, n_max)  — calcule la p-valeur pour chaque taille cumulée
    # Retourne p array (n_simul, n_max - 4) pour n=5..n_max
    n_simul, n_max = g1.shape
    ps = np.ones((n_simul, n_max))
    for n in range(5, n_max + 1):
        a, b = g1[:, :n], g2[:, :n]
        m1, m2 = a.mean(1), b.mean(1)
        s1, s2 = a.var(1, ddof=1), b.var(1, ddof=1)
        se = np.sqrt(s1 / n + s2 / n)
        t = (m1 - m2) / se
        df = (s1/n + s2/n)**2 / ((s1/n)**2/(n-1) + (s2/n)**2/(n-1))
        ps[:, n-1] = 2 * stats.t.sf(np.abs(t), df=df)
    return ps

def phacking_simulation(n_max=100, alpha=0.05, n_simul=500):
    """Simule le p-hacking vectorisé. Retourne le taux de faux positifs."""
    g1 = rng.normal(0, 1, (n_simul, n_max))
    g2 = rng.normal(0, 1, (n_simul, n_max))
    ps = _welch_p(g1, g2)   # (n_simul, n_max)
    # Significatif si au moins un p < alpha en cours de route
    return (ps < alpha).any(axis=1).mean()

taux_phacking = phacking_simulation(n_max=100, n_simul=500)
print(f"Avec test répété jusqu'à n=100 (H₀ vraie) :")
print(f"Taux de 'découvertes' : {taux_phacking:.1%}")
print(f"(devrait être 5% pour un test unique honnête)")
print()

# Trajectoires de p-valeurs (vectorisées)
n_trajectoires = 15
n_max = 80
fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))

ax = axes[0]
g1_traj = rng.normal(0, 1, (n_trajectoires, n_max))
g2_traj = rng.normal(0, 1, (n_trajectoires, n_max))
ps_traj = _welch_p(g1_traj, g2_traj)   # (n_trajectoires, n_max)
ns_range = np.arange(1, n_max + 1)
n_sig = 0
for i in range(n_trajectoires):
    sig = ps_traj[i].min() < 0.05
    color = 'crimson' if sig else 'steelblue'
    alpha_traj = 0.8 if sig else 0.3
    ax.plot(ns_range[4:], ps_traj[i, 4:], color=color, alpha=alpha_traj, lw=1.5)
    if sig:
        n_sig += 1

ax.axhline(0.05, color='black', lw=2, linestyle='--', label='α = 0.05')
ax.set_xlabel('Taille de l\'échantillon')
ax.set_ylabel('p-valeur')
ax.set_title(f'Trajectoires de p-valeurs sous H₀\n(rouge = passé sous 0.05 au moins une fois)')
ax.legend()
ax.set_yscale('log')

# Distribution finale des p-valeurs "sélectionnées" (vectorisé)
ax = axes[1]
n_rep, n_tests_bias, n_obs_bias = 500, 20, 30
a_bias = rng.normal(0, 1, (n_rep, n_tests_bias, n_obs_bias))
b_bias = rng.normal(0, 1, (n_rep, n_tests_bias, n_obs_bias))
m1b, m2b = a_bias.mean(2), b_bias.mean(2)
s1b, s2b = a_bias.var(2, ddof=1), b_bias.var(2, ddof=1)
se_b = np.sqrt(s1b / n_obs_bias + s2b / n_obs_bias)
t_b = (m1b - m2b) / se_b
df_b = (s1b/n_obs_bias + s2b/n_obs_bias)**2 / (
    (s1b/n_obs_bias)**2/(n_obs_bias-1) + (s2b/n_obs_bias)**2/(n_obs_bias-1))
ps_bias = 2 * stats.t.sf(np.abs(t_b), df=df_b)   # (n_rep, n_tests_bias)
p_min_vals = ps_bias.min(axis=1)
p_single_vals = ps_bias[:, 0]

ax.hist(p_single_vals, bins=25, density=True, alpha=0.6, color='steelblue',
        edgecolor='white', label='Test unique honnête')
ax.hist(p_min_vals, bins=25, density=True, alpha=0.6, color='crimson',
        edgecolor='white', label='Minimum de 20 tests (p-hacking)')
ax.axhline(1.0, color='black', lw=1.5, linestyle=':', label='Uniforme(0,1) théorique')
ax.axvline(0.05, color='orange', lw=2, linestyle='--')
ax.set_xlabel('p-valeur')
ax.set_ylabel('Densité')
ax.set_title('Biais du p-hacking :\nl\'exploitation des tests multiples')
ax.legend(fontsize=9)

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

## Correction de Bonferroni

La correction la plus simple : diviser $\alpha$ par le nombre de tests $m$.

$$\alpha_{\text{Bonferroni}} = \frac{\alpha}{m} \qquad \text{ou équivalemment} \quad p_{\text{corr}} = \min(m \cdot p_i, 1)$$

La correction de Bonferroni contrôle le **FWER** (Family-Wise Error Rate) — la probabilité d'au moins un faux positif. Elle est valable même si les tests sont corrélés (contrôle conservateur).

```{admonition} Conservatisme de Bonferroni
:class: warning

Bonferroni est très conservateur lorsque $m$ est grand et que les tests sont corrélés positivement. Elle peut manquer des vrais positifs (faible puissance). Pour des analyses exploratoires ou en génomique, préférer les méthodes FDR.
```

## Holm-Bonferroni : plus puissant que Bonferroni

La procédure de Holm contrôle aussi le FWER mais est uniformément plus puissante que Bonferroni :

1. Trier les p-valeurs : $p_{(1)} \leq p_{(2)} \leq \cdots \leq p_{(m)}$
2. Pour $i = 1, 2, \ldots, m$ : rejeter $H_{(i)}$ si $p_{(j)} \leq \alpha / (m - j + 1)$ pour tout $j \leq i$
3. Dès qu'un test n'est pas rejeté, arrêter

```{code-cell} python3
# Exemple : 10 tests, dont 3 avec vrais effets
rng_ex = np.random.default_rng(7)

# Simuler 10 p-valeurs : 7 sous H0, 3 avec vrais effets
n_tests_tot = 10
n_vrais_effets = 3
p_vals_sim = []
vrai_effet = []

for i in range(n_tests_tot):
    if i < n_vrais_effets:
        # Vrai effet : d = 0.7
        g1 = rng_ex.normal(0, 1, 50)
        g2 = rng_ex.normal(0.7, 1, 50)
        vrai_effet.append(True)
    else:
        # Sous H0
        g1 = rng_ex.normal(0, 1, 50)
        g2 = rng_ex.normal(0, 1, 50)
        vrai_effet.append(False)
    _, p = stats.ttest_ind(g1, g2)
    p_vals_sim.append(p)

df_tests = pd.DataFrame({
    'Test': [f'H_{i+1}' for i in range(n_tests_tot)],
    'p_brute': p_vals_sim,
    'vrai_effet': vrai_effet
})

# Appliquer les corrections
for method, col_name in [('bonferroni', 'p_bonf'), ('holm', 'p_holm'),
                          ('fdr_bh', 'p_fdr_bh'), ('fdr_by', 'p_fdr_by')]:
    reject, pvals_corr, _, _ = multipletests(p_vals_sim, alpha=0.05, method=method)
    df_tests[col_name] = pvals_corr
    df_tests[f'rej_{method}'] = reject

print("Comparaison des méthodes de correction :")
print(df_tests[['Test', 'p_brute', 'p_bonf', 'p_holm', 'p_fdr_bh', 'vrai_effet']].to_string(index=False))
print()
print("Rejets au seuil 0.05 :")
for method, col in [('Sans correction', 'p_brute'), ('Bonferroni', 'p_bonf'),
                     ('Holm', 'p_holm'), ('FDR Benjamini-Hochberg', 'p_fdr_bh')]:
    n_rejets = (df_tests[col] < 0.05).sum()
    n_vrais = sum(r and e for r, e in zip(df_tests[col] < 0.05, df_tests['vrai_effet']))
    n_faux = n_rejets - n_vrais
    print(f"  {method:<30} : {n_rejets} rejets ({n_vrais} vrais positifs, {n_faux} faux positifs)")
```

## FDR — Benjamini-Hochberg

Le **Taux de Faux Découvertes** (False Discovery Rate, FDR) est plus adapté aux analyses à grande échelle. Plutôt que d'interdire tout faux positif, on contrôle la *proportion* de fausses découvertes parmi les rejets.

La procédure de Benjamini-Hochberg :

1. Trier les p-valeurs : $p_{(1)} \leq p_{(2)} \leq \cdots \leq p_{(m)}$
2. Trouver le plus grand $k$ tel que $p_{(k)} \leq \frac{k}{m} \cdot q^*$
3. Rejeter $H_{(1)}, \ldots, H_{(k)}$

où $q^*$ est le niveau FDR cible (souvent 0,05 ou 0,10).

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

# Visualisation de la procédure BH
rng_bh = np.random.default_rng(15)
m_bh = 30
n_vrais_bh = 8

# p-valeurs simulées
p_bh = []
vrais_bh = []
for i in range(m_bh):
    if i < n_vrais_bh:
        g1, g2 = rng_bh.normal(0,1,40), rng_bh.normal(0.6,1,40)
        vrais_bh.append(True)
    else:
        g1, g2 = rng_bh.normal(0,1,40), rng_bh.normal(0,1,40)
        vrais_bh.append(False)
    _, p = stats.ttest_ind(g1, g2)
    p_bh.append(p)

# Trier
order = np.argsort(p_bh)
p_sorted = np.array(p_bh)[order]
vrai_sorted = np.array(vrais_bh)[order]

q_star = 0.05
seuils_bh = np.arange(1, m_bh + 1) / m_bh * q_star

# Trouver k maximal
k_max = 0
for k in range(m_bh):
    if p_sorted[k] <= seuils_bh[k]:
        k_max = k + 1

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

# Graphique de la procédure BH
ax = axes[0]
colors_pts = ['green' if v else 'red' for v in vrai_sorted]
ax.scatter(range(1, m_bh + 1), p_sorted, c=colors_pts, s=60, zorder=5, alpha=0.8)
ax.plot(range(1, m_bh + 1), seuils_bh, 'navy', lw=2, linestyle='--', label=f'Seuil BH (q*={q_star})')
ax.axhline(0.05 / m_bh, color='orange', lw=2, linestyle=':', label=f'Bonferroni ({0.05/m_bh:.4f})')
if k_max > 0:
    ax.axvline(k_max, color='crimson', lw=2, linestyle='--', label=f'Seuil BH : k={k_max}')
    ax.fill_between(range(1, k_max + 1), p_sorted[:k_max], seuils_bh[:k_max],
                    alpha=0.1, color='steelblue')

from matplotlib.lines import Line2D
legend_elements = [
    Line2D([0], [0], marker='o', color='w', markerfacecolor='green', markersize=10, label='Vrai effet'),
    Line2D([0], [0], marker='o', color='w', markerfacecolor='red', markersize=10, label='H₀ vraie'),
    Line2D([0], [0], color='navy', lw=2, linestyle='--', label='Seuil BH'),
    Line2D([0], [0], color='orange', lw=2, linestyle=':', label='Bonferroni'),
]
ax.legend(handles=legend_elements, fontsize=9)
ax.set_xlabel('Rang de la p-valeur')
ax.set_ylabel('p-valeur')
ax.set_title(f'Procédure Benjamini-Hochberg\n({k_max} rejets au niveau FDR = {q_star})')

# Comparaison des corrections : nombres de rejets et FDR
ax = axes[1]
methods = ['Sans correction', 'Bonferroni', 'Holm', 'BH (FDR)', 'BY (FDR dép.)']
scipy_methods = [None, 'bonferroni', 'holm', 'fdr_bh', 'fdr_by']

n_rejets_list = []
n_vp_list = []
n_fp_list = []

for mname, meth in zip(methods, scipy_methods):
    if meth is None:
        rej = np.array(p_bh) < 0.05
    else:
        rej, _, _, _ = multipletests(p_bh, alpha=0.05, method=meth)
    n_r = rej.sum()
    n_vp = sum(r and v for r, v in zip(rej, vrais_bh))
    n_fp = n_r - n_vp
    n_rejets_list.append(n_r)
    n_vp_list.append(n_vp)
    n_fp_list.append(n_fp)

x = np.arange(len(methods))
width = 0.35
bars1 = ax.bar(x - width/2, n_vp_list, width, label='Vrais positifs', color='#55A868', alpha=0.8)
bars2 = ax.bar(x + width/2, n_fp_list, width, label='Faux positifs', color='#C44E52', alpha=0.8)
ax.set_xticks(x)
ax.set_xticklabels([m.replace(' ', '\n') for m in methods], fontsize=9)
ax.set_ylabel('Nombre de rejets')
ax.set_title('Comparaison des méthodes de correction\n(vrais vs faux positifs)')
ax.legend()
ax.bar_label(bars1, padding=2)
ax.bar_label(bars2, padding=2)

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

## Volcano plot : visualiser les résultats à grande échelle

Le *volcano plot* est la visualisation standard en bioinformatique et en pharmacologie pour représenter simultanément la taille d'effet et la significativité statistique.

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

# Simulation d'une étude d'expression génique : 500 gènes
n_genes = 500
n_vrais_genes = 40  # 40 gènes différentiellement exprimés

rng_genes = np.random.default_rng(42)

# Effets : la plupart sont nuls, quelques-uns ont un effet fort
effects = np.zeros(n_genes)
# Gènes sur-exprimés
effects[:20] = rng_genes.normal(1.5, 0.5, 20)
# Gènes sous-exprimés
effects[20:40] = rng_genes.normal(-1.5, 0.5, 20)

p_valeurs_genes = []
log2fc = []

for i in range(n_genes):
    n_echant = 10
    expr_ctrl = rng_genes.normal(0, 1, n_echant)
    expr_cas = rng_genes.normal(effects[i], 1, n_echant)
    _, p = stats.ttest_ind(expr_ctrl, expr_cas)
    p_valeurs_genes.append(p)
    log2fc.append(expr_cas.mean() - expr_ctrl.mean())

p_valeurs_genes = np.array(p_valeurs_genes)
log2fc = np.array(log2fc)

# Correction BH
reject_bh, p_corr_bh, _, _ = multipletests(p_valeurs_genes, alpha=0.05, method='fdr_bh')

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

ax = axes[0]
neg_log10_p = -np.log10(p_valeurs_genes + 1e-20)
colors_vol = ['gray'] * n_genes
for i in range(n_genes):
    if p_valeurs_genes[i] < 0.05 and abs(log2fc[i]) > 1:
        colors_vol[i] = '#C44E52' if log2fc[i] > 0 else '#4C72B0'

ax.scatter(log2fc, neg_log10_p, c=colors_vol, alpha=0.6, s=20)
ax.axhline(-np.log10(0.05), color='black', lw=1.5, linestyle='--', label='p = 0.05')
ax.axvline(-1, color='gray', lw=1, linestyle=':')
ax.axvline(1, color='gray', lw=1, linestyle=':')
ax.set_xlabel('log₂ Fold Change')
ax.set_ylabel('-log₁₀(p-valeur)')
ax.set_title('Volcano plot (sans correction)\n(rouge = sur-exprimé, bleu = sous-exprimé)')
legend_elts = [
    mpatches.Patch(color='#C44E52', label='Sur-exprimé (|FC|>1, p<0.05)'),
    mpatches.Patch(color='#4C72B0', label='Sous-exprimé'),
    mpatches.Patch(color='gray', label='Non significatif'),
]
ax.legend(handles=legend_elts, fontsize=8)

# Après correction FDR
ax = axes[1]
neg_log10_bh = -np.log10(p_corr_bh + 1e-20)
colors_bh = ['gray'] * n_genes
for i in range(n_genes):
    if reject_bh[i] and abs(log2fc[i]) > 1:
        colors_bh[i] = '#C44E52' if log2fc[i] > 0 else '#4C72B0'

ax.scatter(log2fc, neg_log10_bh, c=colors_bh, alpha=0.6, s=20)
ax.axhline(-np.log10(0.05), color='black', lw=1.5, linestyle='--', label='q = 0.05')
ax.axvline(-1, color='gray', lw=1, linestyle=':')
ax.axvline(1, color='gray', lw=1, linestyle=':')
ax.set_xlabel('log₂ Fold Change')
ax.set_ylabel('-log₁₀(q-valeur BH)')
ax.set_title(f'Volcano plot (après correction FDR)\n({reject_bh.sum()} rejets)')
ax.legend(handles=legend_elts, fontsize=8)

plt.tight_layout()
plt.savefig('_static/08_volcano.png', bbox_inches='tight')
plt.show()

n_sig_brut = ((p_valeurs_genes < 0.05) & (np.abs(log2fc) > 1)).sum()
n_sig_fdr = (reject_bh & (np.abs(log2fc) > 1)).sum()
print(f"Sans correction : {n_sig_brut} gènes significatifs (|FC|>1)")
print(f"Après FDR BH   : {n_sig_fdr} gènes significatifs (|FC|>1)")
```

## Tests post-hoc ANOVA

Après une ANOVA significative, on identifie les paires différentes avec des tests post-hoc spécialisés, qui contrôlent les comparaisons multiples par construction.

```{code-cell} python3
# Exemple : 5 traitements médicaux, mesure sur 25 patients chacun
rng_ph = np.random.default_rng(123)
n_pat = 25
traitements = {
    'T1 (contrôle)': rng_ph.normal(50, 10, n_pat),
    'T2': rng_ph.normal(55, 10, n_pat),
    'T3': rng_ph.normal(62, 10, n_pat),
    'T4': rng_ph.normal(58, 10, n_pat),
    'T5': rng_ph.normal(55, 14, n_pat),  # variance plus grande
}

df_ph = pd.DataFrame({
    'score': np.concatenate(list(traitements.values())),
    'traitement': np.repeat(list(traitements.keys()), n_pat)
})

# ANOVA globale
aov_ph = pg.anova(data=df_ph, dv='score', between='traitement')
print("ANOVA globale :")
print(aov_ph[['Source', 'F', 'p_unc', 'np2']].to_string(index=False))
print()

# Tukey HSD
print("Tukey HSD :")
tukey_ph = pg.pairwise_tukey(data=df_ph, dv='score', between='traitement')
print(tukey_ph[['A', 'B', 'diff', 'se', 'T', 'p_tukey']].to_string(index=False))
print()

# Games-Howell (variances inégales)
print("Games-Howell (robuste aux variances inégales) :")
gh = pg.pairwise_gameshowell(data=df_ph, dv='score', between='traitement')
print(gh[['A', 'B', 'diff', 'se', 'T', 'pval']].to_string(index=False))
```

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

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

# Boxplot avec annotations
ax = axes[0]
order_trt = list(traitements.keys())
sns.boxplot(data=df_ph, x='traitement', y='score', order=order_trt, ax=ax, palette='husl')
sns.stripplot(data=df_ph, x='traitement', y='score', order=order_trt, ax=ax,
              color='black', alpha=0.3, size=3)
ax.set_xticklabels([t.replace(' ', '\n') for t in order_trt], fontsize=9)
ax.set_title('Comparaison des traitements\n(ANOVA + post-hoc)')
ax.set_ylabel('Score')

# Heatmap des p-valeurs post-hoc (Tukey)
ax = axes[1]
trt_names = list(traitements.keys())
n_trt = len(trt_names)
pmat = np.ones((n_trt, n_trt))
short_names = [f'T{i+1}' for i in range(n_trt)]

for _, row in tukey_ph.iterrows():
    i = list(traitements.keys()).index(row['A'])
    j = list(traitements.keys()).index(row['B'])
    pmat[i, j] = row['p_tukey']
    pmat[j, i] = row['p_tukey']

np.fill_diagonal(pmat, np.nan)

mask = np.eye(n_trt, dtype=bool)
sns.heatmap(pmat, annot=True, fmt='.3f', cmap='RdYlGn_r',
            xticklabels=short_names, yticklabels=short_names,
            vmin=0, vmax=0.10, ax=ax, mask=mask,
            cbar_kws={'label': 'p-valeur Tukey HSD'})
ax.set_title('Matrice des p-valeurs post-hoc\n(Tukey HSD)')

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

## Reproductibilité et crise de la réplication

### La crise de la réplication

En 2015, le *Reproducibility Project* a tenté de répliquer 100 études en psychologie : seulement ~36% ont produit des résultats significatifs similaires. Les causes identifiées :

- **p-hacking** et *researcher degrees of freedom* (choix d'analyse après avoir vu les données)
- **Biais de publication** : seuls les résultats positifs sont publiés
- **Taille d'échantillon insuffisante** : études sous-puissantes, qui, quand elles "réussissent", surestiment les effets
- **Focalisation sur p < 0,05** plutôt que sur les tailles d'effet

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

# Illustration : biais du "winner's curse"
# Les études significatives avec petit n surestiment les effets

def simulate_winner_curse(true_effect, n_obs, n_simul=500, alpha=0.05):
    """
    Simule le biais du 'winner's curse' : seules les études significatives
    sont publiées. Quelle est l'estimation de l'effet chez les publiées ?
    """
    all_effects = []
    published_effects = []
    for _ in range(n_simul):
        g1 = rng.normal(0, 1, n_obs)
        g2 = rng.normal(true_effect, 1, n_obs)
        _, p = stats.ttest_ind(g1, g2)
        est_effect = g2.mean() - g1.mean()
        all_effects.append(est_effect)
        if p < alpha:
            published_effects.append(est_effect)
    return np.array(all_effects), np.array(published_effects)

true_d = 0.3  # petit effet réel
all_eff, pub_eff = simulate_winner_curse(true_d, n_obs=25, n_simul=500)

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

ax = axes[0]
ax.hist(all_eff, bins=40, density=True, alpha=0.6, color='steelblue',
        edgecolor='white', label=f'Toutes les études (m={len(all_eff)})')
ax.hist(pub_eff, bins=30, density=True, alpha=0.7, color='crimson',
        edgecolor='white', label=f'Études "publiées" p<0.05 (m={len(pub_eff)})')
ax.axvline(true_d, color='black', lw=2.5, linestyle='--', label=f'Vrai effet d={true_d}')
ax.axvline(np.mean(pub_eff), color='crimson', lw=2, linestyle=':', label=f'Effet moyen "publié" = {np.mean(pub_eff):.2f}')
ax.set_xlabel('Effet estimé (différence de moyennes)')
ax.set_ylabel('Densité')
ax.set_title(f'"Winner\'s curse" : surestimation des effets\n(n={25} par groupe, d réel = {true_d})')
ax.legend(fontsize=8)

# Comparaison p-valeur vs intervalle de confiance
ax = axes[1]
ax.axis('off')
conseils = [
    ["Pratique déconseillée", "Alternative recommandée"],
    ["Reporter seulement p-valeur", "Reporter IC 95% + taille d'effet"],
    ["Décision binaire p<0.05", "Interpréter la magnitude et la précision"],
    ["Tests multiples sans correction", "Correction Bonferroni/FDR selon le contexte"],
    ["Collecter des données jusqu'à p<0.05", "Fixer n avant, méthodes séquentielles"],
    ["Omettre les résultats non significatifs", "Pre-enregistrement, open data"],
    ["Test post-hoc sans ajustement", "Tukey HSD, Games-Howell, FDR"],
]
table = ax.table(cellText=conseils[1:], colLabels=conseils[0],
                 loc='center', cellLoc='left')
table.auto_set_font_size(False)
table.set_fontsize(9)
table.scale(1.2, 1.9)
for j in range(2):
    table[0, j].set_facecolor('#2c3e50')
    table[0, j].set_text_props(color='white', fontweight='bold')
# Colonne gauche en rouge clair, droite en vert clair
for i in range(1, len(conseils)):
    table[i, 0].set_facecolor('#fde8e8')
    table[i, 1].set_facecolor('#e8f5e9')
ax.set_title('Bonnes pratiques vs pratiques à éviter', pad=20)

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

```{admonition} Pre-enregistrement et open science
:class: note

Le **pre-enregistrement** consiste à déposer publiquement, avant la collecte des données, les hypothèses précises, la taille d'échantillon et la méthode d'analyse. Cela distingue les analyses *confirmatives* (hypothèses précises, puissance calculée) des analyses *exploratoires* (génération d'hypothèses, résultats préliminaires). Des plateformes comme *OSF* (Open Science Framework), *AsPredicted*, ou *ClinicalTrials.gov* permettent ce dépôt. Le partage des données et du code (open data, open code) est une condition supplémentaire de reproductibilité.
```

## Résumé des méthodes de correction

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

# Tableau récapitulatif des méthodes
recap = pd.DataFrame({
    'Méthode': ['Bonferroni', 'Holm-Bonferroni', 'Benjamini-Hochberg', 'Benjamini-Yekutieli', 'Aucune'],
    'Contrôle': ['FWER', 'FWER', 'FDR', 'FDR (dépendances)', '—'],
    'Puissance': ['Faible', 'Modérée', 'Bonne', 'Faible', 'Maximale (faux positifs++)'],
    'Quand l\'utiliser': [
        'Peu de tests, coût des FP élevé',
        'FWER avec meilleure puissance',
        'Analyses à grande échelle',
        'Tests corrélés positivement',
        'Analyse exploratoire initiale',
    ],
    'scipy method': ['bonferroni', 'holm', 'fdr_bh', 'fdr_by', '—'],
})
print(recap.to_string(index=False))
print()
print("Usage avec statsmodels :")
print("  from statsmodels.stats.multitest import multipletests")
print("  reject, p_corr, _, _ = multipletests(p_vals, alpha=0.05, method='fdr_bh')")
```

```{admonition} Comment choisir sa méthode de correction ?
:class: tip

- **Tests confirmateurs** (hypothèses précises, coût d'un faux positif élevé) → **Bonferroni** ou **Holm**
- **Analyses exploratoires** (génomique, imagerie, psychologie) → **Benjamini-Hochberg (FDR)**
- **Tests fortement corrélés** (SNPs en déséquilibre de liaison, pixels voisins) → **Benjamini-Yekutieli**
- **Post-hoc ANOVA**, toutes variances égales → **Tukey HSD**
- **Post-hoc ANOVA**, variances inégales → **Games-Howell**
- **Post-hoc ANOVA**, comparaison vs contrôle uniquement → **Test de Dunnett**
```