Comparez plusieurs groupes avec Kruskal-Wallis et Friedman

Une entreprise hésite entre trois prestataires de livraison et mesure les délais de cinq colis chez chacun. Une équipe de data science compare quatre modèles sur six jeux de données. Dans les deux cas, il y a plus de deux groupes à comparer, et l’on voudrait savoir, d’un seul test, si au moins un groupe se distingue.

En statistique paramétrique, c’est le rôle de l’analyse de la variance (ANOVA). Ce chapitre présente ses deux équivalents de rang : le test de Kruskal-Wallis pour des groupes indépendants, et le test de Friedman pour des mesures organisées en blocs. Vous verrez aussi comment trouver quels groupes diffèrent, sans multiplier les faux positifs.

Plusieurs groupes indépendants : Kruskal-Wallis

Le fil rouge : trois prestataires

Délais de livraison (en heures) de cinq colis par prestataire :

Prestataire A	32	41	28	35	39
Prestataire B	45	52	38	49	56
Prestataire C	30	27	36	33	25

Pourquoi ne pas faire trois tests de Mann-Whitney, A contre B, A contre C, B contre C ? Parce que chaque test a 5 % de chances de conclure à tort. Avec trois tests, le risque de se tromper au moins une fois monte jusqu’à $1 - 0{,}95^3 \approx 14\ \%$ (si les tests étaient indépendants). Il faut d’abord un test global.

La statistique H

Le principe prolonge celui de Mann-Whitney : on classe les $N$ observations ensemble, puis on compare les rangs moyens des groupes. Sous $H_0$ (tous les groupes viennent de la même loi), chaque rang moyen devrait être proche du rang moyen global $(N+1)/2$. La statistique de Kruskal-Wallis mesure l’écart à cette situation :

$$H = \frac{12}{N(N+1)} \sum_{j=1}^{k} \frac{R_j^2}{n_j} - 3(N+1)$$

• $k$ est le nombre de groupes et $n_j$ l’effectif du groupe $j$ ;
• $N = n_1 + \dots + n_k$ est l’effectif total ;
• $R_j$ est la somme des rangs du groupe $j$ dans le classement global.

Sous $H_0$, et si chaque groupe compte au moins cinq observations environ, $H$ suit approximativement une loi du $\chi^2$ à $k - 1$ degrés de liberté. On rejette $H_0$ quand $H$ est grand.

Exemple

En classant les 15 délais, les rangs sont :

• A : 5, 11, 3, 7, 10, soit $R_A = 36$ (rang moyen 7,2) ;
• B : 12, 14, 9, 13, 15, soit $R_B = 63$ (rang moyen 12,6) ;
• C : 4, 2, 8, 6, 1, soit $R_C = 21$ (rang moyen 4,2).

Le rang moyen global vaut $(15+1)/2 = 8$. Alors :

$$H = \frac{12}{15 \times 16} \left( \frac{36^2 + 63^2 + 21^2}{5} \right) - 3 \times 16 = 0{,}05 \times 1\,141{,}2 - 48 = 9{,}06$$

Avec $k - 1 = 2$ degrés de liberté, $p = P(\chi^2_2 \ge 9{,}06) = e^{-9{,}06/2} \approx 0{,}011$. On rejette $H_0$ : au moins un prestataire a des délais différents.

Sous H₀, les trois rangs moyens devraient tourner autour de 8. Celui de B est bien au-dessus et celui de C bien en dessous : H additionne ces écarts au carré.

À retenir

Kruskal-Wallis est l’ANOVA à un facteur des rangs. Il compare les rangs moyens de $k$ groupes indépendants ; sous $H_0$, $H$ suit approximativement un $\chi^2_{k-1}$. Pour $k = 2$, il équivaut à Mann-Whitney bilatéral (en approximation normale).

H = 9,06 dépasse le seuil de 5,99 ; l'aire restante à droite, p ≈ 0,011, est la p-valeur.

Une mesure de taille d’effet simple est $\varepsilon^2 = H/(N-1)$, entre 0 et 1 : ici $9{,}06/14 \approx 0{,}65$, un effet très fort.

Trouver quels groupes diffèrent

Le test global dit « au moins un groupe diffère », sans dire lequel. On enchaîne alors des comparaisons deux à deux (post hoc), en corrigeant le seuil pour le nombre de comparaisons. Le test de Dunn compare les rangs moyens de deux groupes, en réutilisant le classement global :

$$z_{ij} = \frac{\bar R_i - \bar R_j}{\sqrt{\dfrac{N(N+1)}{12}\left(\dfrac{1}{n_i} + \dfrac{1}{n_j}\right)}}$$

où $\bar R_i$ est le rang moyen du groupe $i$. On compare $z_{ij}$ à la loi normale, puis on corrige les p-valeurs, par exemple avec la méthode de Bonferroni (multiplier chaque p-valeur par le nombre de comparaisons) ou celle, moins conservatrice, de Holm.

Comparaison	$z$	$p$ brute	$p$ Bonferroni (× 3)
B contre C	2,97	0,003	0,009
A contre B	−1,91	0,056	0,17
A contre C	1,06	0,29	0,87

Seul l’écart entre B (le plus lent) et C (le plus rapide) est significatif après correction.

Attention

Comme Mann-Whitney, Kruskal-Wallis compare des médianes seulement si les lois des groupes ont la même forme et ne diffèrent que par un décalage. Sinon, il détecte plus largement une tendance d’un groupe à produire des valeurs plus grandes que les autres.

Des mesures en blocs : Friedman

Le fil rouge : quatre modèles, six jeux de données

Quatre modèles de classification (régression logistique, arbre de décision, forêt aléatoire, boosting) sont évalués sur six jeux de données. Voici leurs taux de bonnes réponses :

Jeu	LogReg	Arbre	Forêt	Boosting
1	0,81	0,76	0,85	0,86
2	0,72	0,70	0,78	0,77
3	0,90	0,84	0,91	0,93
4	0,65	0,66	0,71	0,73
5	0,88	0,82	0,87	0,90
6	0,77	0,71	0,80	0,79

Les modèles ne sont pas des groupes indépendants : chaque jeu de données est un bloc, sur lequel les quatre modèles sont mesurés. Certains jeux sont faciles (autour de 0,90), d’autres difficiles (autour de 0,70). Si l’on ignorait les blocs, cette variabilité entre jeux noierait les différences entre modèles. C’est le même raisonnement que pour les données appariées du chapitre 2, étendu à plus de deux traitements.

La statistique Q

Le test de Friedman classe les traitements à l’intérieur de chaque bloc (ici, 1 pour le meilleur modèle du jeu, 4 pour le moins bon), puis additionne les rangs par traitement. Avec $n$ blocs et $k$ traitements :

$$Q = \frac{12}{nk(k+1)} \sum_{j=1}^{k} R_j^2 - 3n(k+1)$$

où $R_j$ est la somme des rangs du traitement $j$ sur les $n$ blocs. Sous $H_0$ (dans chaque bloc, tous les classements sont équiprobables), $Q$ suit approximativement un $\chi^2_{k-1}$.

Exemple

Rangs dans chaque jeu (1 = meilleur) :

Jeu	LogReg	Arbre	Forêt	Boosting
1	3	4	2	1
2	3	4	1	2
3	3	4	2	1
4	4	3	2	1
5	2	4	3	1
6	3	4	1	2
Somme $R_j$	18	23	11	8
Rang moyen	3,00	3,83	1,83	1,33

$$Q = \frac{12}{6 \times 4 \times 5} (18^2 + 23^2 + 11^2 + 8^2) - 3 \times 6 \times 5 = 0{,}1 \times 1\,038 - 90 = 13{,}8$$

Avec 3 degrés de liberté, $p \approx 0{,}003$. Les modèles ne se valent pas.

Sur les valeurs brutes, l'écart entre jeux faciles et difficiles est plus grand que l'écart entre modèles. Classer à l'intérieur de chaque jeu efface cet effet de bloc et fait ressortir un ordre stable : boosting et forêt devant, arbre derrière.

La quantité $W = Q/(n(k-1))$, ici $13{,}8/18 \approx 0{,}77$, s’appelle le coefficient de concordance de Kendall : elle mesure à quel point les jeux de données s’accordent sur le classement des modèles (0 : aucun accord, 1 : même classement partout).

Le diagramme de différence critique

Pour localiser les différences entre modèles, la démarche classique en apprentissage automatique (popularisée par un article de Demšar en 2006) utilise le test post hoc de Nemenyi : deux modèles diffèrent significativement si leurs rangs moyens s’écartent d’au moins la différence critique

$$CD = q_\alpha \sqrt{\frac{k(k+1)}{6n}}$$

où $q_\alpha$ est une valeur tabulée (2,569 pour $k = 4$ et $\alpha = 5\ \%$). Ici, $CD = 2{,}569 \times \sqrt{20/36} \approx 1{,}91$. Le boosting (1,33) et la forêt (1,83) se distinguent de l’arbre (3,83), mais aucune autre paire ne dépasse l’écart de 1,91 : on ne peut pas départager boosting et forêt, ni boosting et régression logistique (écart de 1,67). Six jeux de données, c’est peu.

Deux modèles reliés par une barre ne se distinguent pas (écart de rang moyen inférieur à CD). Boosting et forêt se séparent de l'arbre, mais la régression logistique reste reliée aux deux groupes.

En entretien

« Pour comparer plusieurs modèles, je les évalue sur plusieurs jeux de données, je les classe sur chaque jeu, puis j’applique le test de Friedman et, s’il rejette, un post hoc de Nemenyi avec un diagramme de différence critique. Cela évite de conclure sur un seul jeu de données ou avec des tests multiples non corrigés. »

En pratique : Python et R

from scipy import stats

A = [32, 41, 28, 35, 39]
B = [45, 52, 38, 49, 56]
C = [30, 27, 36, 33, 25]
print(stats.kruskal(A, B, C))      # H = 9.06, p ≈ 0.011

logreg   = [0.81, 0.72, 0.90, 0.65, 0.88, 0.77]
arbre    = [0.76, 0.70, 0.84, 0.66, 0.82, 0.71]
foret    = [0.85, 0.78, 0.91, 0.71, 0.87, 0.80]
boosting = [0.86, 0.77, 0.93, 0.73, 0.90, 0.79]
# un argument par traitement, les blocs sont les positions dans les listes
print(stats.friedmanchisquare(logreg, arbre, foret, boosting))   # Q = 13.8, p ≈ 0.003

Les tests post hoc (Dunn, Nemenyi) ne sont pas dans scipy ; le paquet scikit-posthocs les fournit (posthoc_dunn, posthoc_nemenyi_friedman). En R :

delais <- c(32, 41, 28, 35, 39,  45, 52, 38, 49, 56,  30, 27, 36, 33, 25)
presta <- factor(rep(c("A", "B", "C"), each = 5))
kruskal.test(delais ~ presta)          # H = 9.06, df = 2, p = 0.0108

acc <- matrix(c(0.81, 0.76, 0.85, 0.86,
                0.72, 0.70, 0.78, 0.77,
                0.90, 0.84, 0.91, 0.93,
                0.65, 0.66, 0.71, 0.73,
                0.88, 0.82, 0.87, 0.90,
                0.77, 0.71, 0.80, 0.79),
              nrow = 6, byrow = TRUE,
              dimnames = list(NULL, c("LogReg", "Arbre", "Foret", "Boosting")))
friedman.test(acc)                     # une ligne par bloc : Q = 13.8, df = 3, p = 0.0032

Piège

Le sens des rangs de Friedman n’a pas d’importance pour le test (classer du meilleur au pire ou l’inverse donne le même $Q$), mais il en a pour l’interprétation des rangs moyens. Précisez toujours si « rang 1 » veut dire « meilleur » ou « moins bon ».

Récapitulatif des tests de rang

Plan d’expérience	Test paramétrique	Test de rang
Un échantillon, ou deux mesures appariées	t à un échantillon ou apparié	signe, Wilcoxon signé
Deux groupes indépendants	t de Student ou de Welch	Mann-Whitney
$k$ groupes indépendants	ANOVA à un facteur	Kruskal-Wallis
$k$ traitements mesurés sur $n$ blocs	ANOVA à mesures répétées	Friedman

En résumé

• Pour plus de deux groupes, on commence par un test global pour éviter l’inflation des faux positifs.
• Kruskal-Wallis compare les rangs moyens de $k$ groupes indépendants : $H \approx \chi^2_{k-1}$ sous $H_0$.
• Friedman classe les traitements à l’intérieur de chaque bloc : $Q \approx \chi^2_{k-1}$ sous $H_0$ ; il neutralise la variabilité entre blocs.
• Après un rejet, les comparaisons post hoc (Dunn pour Kruskal-Wallis, Nemenyi pour Friedman) localisent les différences, avec une correction pour tests multiples.
• Pour comparer des modèles sur plusieurs jeux de données : Friedman puis Nemenyi, résumés par un diagramme de différence critique.
• scipy.stats.kruskal, scipy.stats.friedmanchisquare ; en R, kruskal.test et friedman.test.

Exercices

Exercice · Un test de Friedman à la main

facile

Quatre juges classent trois recettes de cookies (1 = préférée). Classements : juge 1 : (1, 2, 3) ; juge 2 : (1, 3, 2) ; juge 3 : (1, 2, 3) ; juge 4 : (2, 1, 3), pour les recettes (X, Y, Z). Calculez $Q$ et concluez sachant que $P(\chi^2_2 \ge q) = e^{-q/2}$.

Voir la solution

Sommes des rangs : $R_X = 1 + 1 + 1 + 2 = 5$, $R_Y = 2 + 3 + 2 + 1 = 8$, $R_Z = 3 + 2 + 3 + 3 = 11$ (total $24 = 4 \times 6$, c’est cohérent). Avec $n = 4$ et $k = 3$ :

$$Q = \frac{12}{4 \times 3 \times 4}(25 + 64 + 121) - 3 \times 4 \times 4 = 0{,}25 \times 210 - 48 = 4{,}5$$

$p \approx e^{-2{,}25} \approx 0{,}105$. On ne rejette pas $H_0$ au seuil de 5 %. La recette X semble préférée, mais quatre juges ne suffisent pas pour l’affirmer. (Avec aussi peu de blocs, l’approximation du $\chi^2$ est d’ailleurs grossière ; des tables exactes existent.)

En code

Python

import numpy as np
from scipy import stats

# Une ligne par juge, une colonne par recette (X, Y, Z) ; 1 = préférée
rangs = np.array([[1, 2, 3], [1, 3, 2], [1, 2, 3], [2, 1, 3]])
print(rangs.sum(axis=0))                   # [ 5  8 11]

res = stats.friedmanchisquare(*rangs.T)    # un argument par recette
print(res.statistic, res.pvalue)           # 4.5 0.105

R

# Une ligne par juge (bloc), une colonne par recette (X, Y, Z) ; 1 = préférée
rangs <- matrix(c(1, 2, 3,  1, 3, 2,  1, 2, 3,  2, 1, 3), nrow = 4, byrow = TRUE,
                dimnames = list(NULL, c("X", "Y", "Z")))
colSums(rangs)                  # 5 8 11
friedman.test(rangs)            # Friedman chi-squared = 4.5, df = 2, p-value = 0.1054

Exercice · Petits groupes et approximation

moyen

Trois groupes de trois mesures : G1 = (1,2 ; 3,4 ; 2,2), G2 = (4,1 ; 5,0 ; 3,9), G3 = (0,8 ; 2,9 ; 1,5).

1. Calculez $H$.
2. L’approximation du $\chi^2_2$ donne $p \approx 0{,}061$. Par énumération des $1\,680$ répartitions possibles des rangs, on trouve 84 répartitions avec $H \ge$ la valeur observée. Quelle est la p-valeur exacte ? Qu’en concluez-vous ?

Voir la solution

1. Classement des 9 valeurs : 0,8 (1), 1,2 (2), 1,5 (3), 2,2 (4), 2,9 (5), 3,4 (6), 3,9 (7), 4,1 (8), 5,0 (9). Sommes de rangs : G1 = 2 + 6 + 4 = 12, G2 = 8 + 9 + 7 = 24, G3 = 1 + 5 + 3 = 9. Alors :

$$H = \frac{12}{9 \times 10} \times \frac{144 + 576 + 81}{3} - 3 \times 10 = \frac{12}{90} \times 267 - 30 = 5{,}6$$

2. $p_{\text{exact}} = 84/1\,680 = 0{,}05$. L’approximation (0,061) est ici un peu conservatrice. Avec de si petits groupes, il faut utiliser la loi exacte (tables ou énumération) ; l’approximation du $\chi^2$ n’est fiable qu’à partir d’environ cinq observations par groupe.

En code

Python

import itertools
import numpy as np
from scipy import stats

G1, G2, G3 = [1.2, 3.4, 2.2], [4.1, 5.0, 3.9], [0.8, 2.9, 1.5]
res = stats.kruskal(G1, G2, G3)
print(res.statistic, res.pvalue)           # 5.6 0.061 (approximation du chi2)

def H(r1, r2, r3, N=9):
    return 12 / (N * (N + 1)) * sum(sum(r) ** 2 / 3 for r in (r1, r2, r3)) - 3 * (N + 1)

# Loi exacte : toutes les répartitions des rangs 1..9 en trois groupes de 3
rangs = set(range(1, 10))
compte = total = 0
for r1 in itertools.combinations(rangs, 3):
    for r2 in itertools.combinations(rangs - set(r1), 3):
        r3 = rangs - set(r1) - set(r2)
        total += 1
        compte += H(r1, r2, r3) >= 5.6 - 1e-9
print(total, compte, compte / total)       # 1680 84 0.05

R

G1 <- c(1.2, 3.4, 2.2); G2 <- c(4.1, 5.0, 3.9); G3 <- c(0.8, 2.9, 1.5)
kruskal.test(list(G1, G2, G3))   # chi-squared = 5.6, df = 2, p-value = 0.06081

H <- function(r1, r2, r3, N = 9)
  12 / (N * (N + 1)) * (sum(r1)^2 + sum(r2)^2 + sum(r3)^2) / 3 - 3 * (N + 1)

# Loi exacte : toutes les répartitions des rangs 1..9 en trois groupes de 3
compte <- 0; total <- 0
for (i in seq_len(choose(9, 3))) {
  r1 <- combn(9, 3)[, i]; reste <- setdiff(1:9, r1)
  for (j in seq_len(choose(6, 3))) {
    r2 <- reste[combn(6, 3)[, j]]; r3 <- setdiff(reste, r2)
    total <- total + 1
    compte <- compte + (H(r1, r2, r3) >= 5.6 - 1e-9)
  }
}
c(total, compte, compte / total)   # 1680 84 0.05

Exercice · Le bon test pour le bon plan

difficile

Pour chaque situation, choisissez le test adapté et justifiez.

1. Trente clients notent de 1 à 5 trois versions d’une application : chaque client teste les trois.
2. On compare la durée de séjour à l’hôpital (très asymétrique) de patients répartis au hasard entre quatre protocoles.
3. On compare deux modèles sur les 10 plis d’une même validation croisée.
4. Un collègue a comparé quatre modèles sur un seul jeu de données en faisant les six tests de Mann-Whitney deux à deux entre leurs scores par pli, sans correction. Qu’en pensez-vous ?

Voir la solution

1. Friedman : chaque client est un bloc, les trois versions sont les traitements, et les notes ordinales se prêtent aux rangs.
2. Kruskal-Wallis : quatre groupes indépendants, variable asymétrique. En cas de rejet, test de Dunn avec correction de Holm ou de Bonferroni.
3. Wilcoxon signé sur les 10 différences de score (chapitre 2), puisque les deux modèles sont évalués sur les mêmes plis. Avec la réserve que les plis ne sont pas indépendants, ce qui rend la p-valeur optimiste.
4. Deux erreurs. D’abord, les scores par pli sont appariés (mêmes plis pour tous les modèles), donc Mann-Whitney, qui suppose des groupes indépendants, est inadapté : il faudrait Friedman, avec les plis comme blocs. Ensuite, six tests non corrigés font grimper le risque de faux positif bien au-delà de 5 %. Enfin, un seul jeu de données ne permet pas de généraliser la conclusion à d’autres problèmes.

En code

Python

import itertools
import numpy as np
from scipy import stats

# Question 4 : quatre modèles équivalents (H0 vraie), 10 plis chacun,
# six tests de Mann-Whitney deux à deux sans correction
rng = np.random.default_rng(0)
n_sim, faux_positif = 2000, 0
for _ in range(n_sim):
    scores = rng.normal(size=(4, 10))
    p = [stats.mannwhitneyu(scores[i], scores[j]).pvalue
         for i, j in itertools.combinations(range(4), 2)]
    faux_positif += min(p) < 0.05
print(faux_positif / n_sim)   # environ 0.19 : près de quatre fois le risque annoncé de 5 %

R

# Question 4 : quatre modèles équivalents (H0 vraie), 10 plis chacun,
# six tests de Mann-Whitney deux à deux sans correction
set.seed(0)
paires <- combn(4, 2)
faux_positif <- replicate(2000, {
  scores <- matrix(rnorm(40), nrow = 4)
  p <- apply(paires, 2, function(k) wilcox.test(scores[k[1], ], scores[k[2], ])$p.value)
  min(p) < 0.05
})
mean(faux_positif)   # environ 0.20 : près de quatre fois le risque annoncé de 5 %