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

Show code cell source

Hide code cell source

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\%\).

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

Nombre de tests | Taux FWER simulé | Taux FWER théorique
-------------------------------------------------------
              1 |            0.042 |               0.050
              5 |            0.224 |               0.226
             10 |            0.400 |               0.401
             20 |            0.622 |               0.642
             50 |            0.908 |               0.923
            100 |            0.994 |               0.994

Show code cell source

Hide code cell source

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

Show code cell source

Hide code cell source

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

Avec test répété jusqu'à n=100 (H₀ vraie) :
Taux de 'découvertes' : 31.6%
(devrait être 5% pour un test unique honnête)

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

Conservatisme de Bonferroni

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

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

Comparaison des méthodes de correction :
Test      p_brute   p_bonf   p_holm  p_fdr_bh  vrai_effet
 H_1 5.459427e-07 0.000005 0.000005  0.000005        True
 H_2 2.775236e-05 0.000278 0.000222  0.000093        True
 H_3 8.272765e-06 0.000083 0.000074  0.000041        True
 H_4 2.968676e-01 1.000000 1.000000  0.424097       False
 H_5 2.038460e-01 1.000000 1.000000  0.339743       False
 H_6 3.809623e-02 0.380962 0.266674  0.095241       False
 H_7 7.886507e-01 1.000000 1.000000  0.876279       False
 H_8 3.709269e-01 1.000000 1.000000  0.463659       False
 H_9 7.019461e-02 0.701946 0.421168  0.140389       False
H_10 8.833582e-01 1.000000 1.000000  0.883358       False

Rejets au seuil 0.05 :
  Sans correction                : 4 rejets (3 vrais positifs, 1 faux positifs)
  Bonferroni                     : 3 rejets (3 vrais positifs, 0 faux positifs)
  Holm                           : 3 rejets (3 vrais positifs, 0 faux positifs)
  FDR Benjamini-Hochberg         : 3 rejets (3 vrais positifs, 0 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).

Show code cell source

Hide code cell source

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

Show code cell source

Hide code cell source

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

Sans correction : 37 gènes significatifs (|FC|>1)
Après FDR BH   : 16 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.

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

ANOVA globale :
    Source        F    p_unc      np2
traitement 3.047068 0.019695 0.092204

Tukey HSD :
            A  B       diff       se         T  p_tukey
T1 (contrôle) T2  -4.869111 3.057277 -1.592630 0.505027
T1 (contrôle) T3 -10.333994 3.057277 -3.380130 0.008524
T1 (contrôle) T4  -6.578402 3.057277 -2.151720 0.205482
T1 (contrôle) T5  -3.947462 3.057277 -1.291169 0.697144
           T2 T3  -5.464882 3.057277 -1.787500 0.385633
           T2 T4  -1.709291 3.057277 -0.559089 0.980605
           T2 T5   0.921650 3.057277  0.301461 0.998179
           T3 T4   3.755591 3.057277  1.228411 0.734850
           T3 T5   6.386532 3.057277  2.088961 0.231633
           T4 T5   2.630941 3.057277  0.860550 0.910629

Games-Howell (robuste aux variances inégales) :

            A  B       diff       se         T     pval
T1 (contrôle) T2  -4.869111 2.844908 -1.711518 0.436823
T1 (contrôle) T3 -10.333994 2.546409 -4.058262 0.001663
T1 (contrôle) T4  -6.578402 2.467859 -2.665631 0.075131
T1 (contrôle) T5  -3.947462 3.716810 -1.062056 0.824545
           T2 T3  -5.464882 2.677213 -2.041258 0.263319
           T2 T4  -1.709291 2.602614 -0.656759 0.964451
           T2 T5   0.921650 3.807617  0.242054 0.999209
           T3 T4   3.755591 2.272508  1.652619 0.472309
           T3 T5   6.386532 3.590075  1.778941 0.400982
           T4 T5   2.630941 3.534794  0.744298 0.944396

Show code cell source

Hide code cell source

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

Show code cell source

Hide code cell source

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

Pre-enregistrement et open science

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

Show code cell source

Hide code cell source

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

            Méthode          Contrôle                  Puissance                Quand l'utiliser scipy method
         Bonferroni              FWER                     Faible Peu de tests, coût des FP élevé   bonferroni
    Holm-Bonferroni              FWER                    Modérée   FWER avec meilleure puissance         holm
 Benjamini-Hochberg               FDR                      Bonne       Analyses à grande échelle       fdr_bh
Benjamini-Yekutieli FDR (dépendances)                     Faible     Tests corrélés positivement       fdr_by
             Aucune                 — Maximale (faux positifs++)   Analyse exploratoire initiale            —

Usage avec statsmodels :
  from statsmodels.stats.multitest import multipletests
  reject, p_corr, _, _ = multipletests(p_vals, alpha=0.05, method='fdr_bh')

Comment choisir sa méthode de correction ?

• 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