Comparez plusieurs groupes avec l'analyse de la variance

Au chapitre 7, vous avez comparé deux moyennes. Mais comment comparer trois modes de transport, quatre types de viande dans des saucisses, ou cinq recettes d’un même produit ? Multiplier les tests de Student deux à deux est tentant, et c’est une erreur : plus on fait de tests, plus on finit par trouver une différence par hasard.

L’analyse de la variance (ANOVA) répond à la question globale en un seul test : les moyennes des groupes sont-elles toutes égales ? Malgré son nom, elle compare bien des moyennes, en analysant comment se répartit la variabilité. C’est aussi un cas particulier du modèle linéaire des deux chapitres précédents. Notre fil rouge : des temps de trajet domicile-école.

Le fil rouge : trois façons de venir à l’INSA

On relève le temps de trajet (en minutes) d’étudiants habitant à une distance comparable du campus, selon leur mode de transport :

Mode	Temps de trajet (min)	Effectif	Moyenne
Bus	28, 34, 25, 31, 37	5	31,0
Vélo	18, 22, 15, 20, 24, 19	6	19,7
Voiture	24, 30, 21, 27, 23	5	25,0

La moyenne générale vaut $\bar y = 398/16 \approx 24{,}9$ minutes. Des boîtes à moustaches côte à côte (boxplot(temps ~ mode) en R) montrent des groupes bien décalés. Mais avec cinq ou six étudiants par groupe, ces écarts pourraient-ils venir du hasard ?

Les groupes sont nettement décalés, mais chacun est aussi dispersé. L'ANOVA confronte ces deux sources de variabilité.

Le vocabulaire : le mode de transport est un facteur (une variable explicative qualitative), et ses $I = 3$ valeurs sont ses modalités (ou niveaux).

Le modèle d’ANOVA à un facteur

On note $Y_{ij}$ le temps du $j$-ème étudiant du groupe $i$, avec $i = 1, \dots, I$ et $j = 1, \dots, n_i$ ; l’effectif total est $n = \sum_i n_i$. Le modèle s’écrit

$$Y_{ij} = \mu_i + \varepsilon_{ij}, \qquad \varepsilon_{ij} \text{ i.i.d. } \mathcal{N}(0, \sigma^2)$$

où $\mu_i$ est la moyenne du groupe $i$. Les hypothèses sont celles de la régression : erreurs indépendantes, gaussiennes et de même variance dans tous les groupes.

On l’écrit souvent sous la forme des effets :

$$Y_{ij} = \mu + \alpha_i + \varepsilon_{ij}$$

où $\mu$ est un niveau de référence et $\alpha_i$ l’effet de la modalité $i$. Mais ce modèle a $I + 1$ paramètres pour $I$ moyennes : il n’est pas identifiable (on peut ajouter 5 à $\mu$ et retirer 5 à tous les $\alpha_i$ sans rien changer). Il faut une contrainte :

• la contrainte « naturelle » $\sum_i n_i \alpha_i = 0$ : $\mu$ est alors la moyenne générale, et $\alpha_i$ l’écart du groupe $i$ à cette moyenne ;
• la contrainte de R par défaut, $\alpha_1 = 0$ : le premier groupe (par ordre alphabétique) sert de référence, $\mu$ est sa moyenne, et chaque $\alpha_i$ est l’écart à ce groupe.

Attention

Les coefficients affichés par summary(lm(temps ~ mode)) dépendent de la contrainte. Ici, R prend « bus » comme référence : la constante vaut 31,0 (moyenne du bus), le coefficient modevelo vaut $-11{,}33$ et modevoiture vaut $-6{,}00$. Ce ne sont pas les moyennes des groupes, mais leurs écarts au bus. Les prévisions, elles, ne dépendent pas de la contrainte.

Décomposer la variabilité

L’idée centrale : si les modes de transport ont des effets différents, la variabilité entre les groupes sera grande devant la variabilité à l’intérieur des groupes. On décompose la variabilité totale en deux morceaux :

$$\underbrace{\sum_{i}\sum_{j} (y_{ij} - \bar y)^2}_{SST} = \underbrace{\sum_{i} n_i (\bar y_i - \bar y)^2}_{SSM} + \underbrace{\sum_{i}\sum_{j} (y_{ij} - \bar y_i)^2}_{SSR}$$

• $SST$ (totale) mesure la dispersion de toutes les observations autour de la moyenne générale $\bar y$ ;
• $SSM$ (du modèle, ou inter-groupes) mesure la dispersion des moyennes de groupes $\bar y_i$ autour de $\bar y$, pondérée par les effectifs ;
• $SSR$ (résiduelle, ou intra-groupes) mesure la dispersion de chaque observation autour de la moyenne de son groupe.

Les degrés de liberté se décomposent de la même façon : $n - 1 = (I - 1) + (n - I)$.

Exemple

• $SSM = 5(31 - 24{,}875)^2 + 6(19{,}67 - 24{,}875)^2 + 5(25 - 24{,}875)^2 \approx 187{,}6 + 162{,}8 + 0{,}1 \approx 350{,}4$.
• $SSR$ : on somme les carrés des écarts dans chaque groupe, soit $90 + 49{,}33 + 50 \approx 189{,}3$.
• $SST \approx 539{,}75$, et on vérifie que $350{,}4 + 189{,}3 = 539{,}75$.

Le mode de transport explique donc $R^2 = SSM/SST \approx 65\,\%$ de la variabilité des temps de trajet.

Chaque écart total (en haut) se coupe en un écart entre la moyenne du groupe et la moyenne générale (au milieu) et un écart à l'intérieur du groupe (en bas). En sommant les carrés : 539,75 = 350,4 + 189,3.

Le test de Fisher

On teste $H_0 : \mu_1 = \mu_2 = \dots = \mu_I$ (le facteur n’a pas d’effet) contre $H_1$ : au moins deux moyennes diffèrent. On compare les carrés moyens, c’est-à-dire les sommes de carrés divisées par leurs degrés de liberté :

$$F_0 = \frac{SSM / (I - 1)}{SSR / (n - I)} \sim \mathcal{F}_{I-1,\ n-I} \text{ sous } H_0$$

Sous $H_0$, numérateur et dénominateur estiment tous deux $\sigma^2$, et $F_0$ reste proche de 1. Si les moyennes diffèrent, le numérateur gonfle. On rejette $H_0$ si $F_0 > f_{I-1, n-I; 1-\alpha}$.

À moyennes de groupe identiques, tout dépend de la dispersion interne : si elle grandit, les écarts entre moyennes ne se distinguent plus du bruit et F s'approche de 1.

On présente le tout dans la table d’ANOVA, exactement comme l’affiche anova() en R :

Source	ddl	Somme des carrés	Carré moyen	$F$	p-valeur
Mode (inter)	2	350,4	175,2	12,03	0,0011
Résidus (intra)	13	189,3	14,56

Avec $f_{2,13;0,95} \approx 3{,}81$, on rejette $H_0$ : le temps de trajet moyen dépend du mode de transport. Le carré moyen résiduel $CMR = 14{,}56$ estime $\sigma^2$, soit $\hat\sigma \approx 3{,}8$ minutes.

Sous H₀, F dépasse 3,81 une fois sur vingt ; la valeur observée 12,03 est bien plus loin, d'où la p-valeur de 0,0011.

trajets <- data.frame(
  temps = c(28, 34, 25, 31, 37, 18, 22, 15, 20, 24, 19, 24, 30, 21, 27, 23),
  mode  = factor(rep(c("bus", "velo", "voiture"), times = c(5, 6, 5)))
)
modele <- lm(temps ~ mode, data = trajets)
anova(modele)      # table d'ANOVA : F = 12.03, p = 0.0011
summary(modele)    # coefficients avec « bus » comme référence
par(mfrow = c(2, 2)); plot(modele)   # graphiques de diagnostic des résidus

from scipy import stats

bus = [28, 34, 25, 31, 37]
velo = [18, 22, 15, 20, 24, 19]
voiture = [24, 30, 21, 27, 23]
print(stats.f_oneway(bus, velo, voiture))   # statistic 12.03, pvalue 0.0011

À retenir

L’ANOVA à un facteur teste l’égalité de $I$ moyennes en comparant la variabilité entre groupes ($SSM$) à la variabilité dans les groupes ($SSR$) : $F_0 = \frac{SSM/(I-1)}{SSR/(n-I)} \sim \mathcal{F}_{I-1,n-I}$ sous $H_0$. Elle suppose des groupes indépendants, gaussiens et de même variance.

Avec deux groupes seulement, l’ANOVA redonne exactement le test de Student à variances égales du chapitre 7, avec $F_0 = t_0^2$. Et c’est bien une régression : lm(temps ~ mode) crée automatiquement des variables indicatrices (0 ou 1) pour chaque modalité sauf la référence.

Après le rejet : quels groupes diffèrent ?

Le test global dit qu’au moins deux moyennes diffèrent, pas lesquelles. On compare alors les groupes deux à deux par

$$T_{ij} = \frac{\bar Y_i - \bar Y_j}{\sqrt{CMR\left(\dfrac{1}{n_i} + \dfrac{1}{n_j}\right)}} \sim \mathcal{T}_{n-I} \text{ sous } H_0 : \mu_i = \mu_j$$

C’est un test de Student, mais qui utilise l’estimation de $\sigma^2$ fondée sur tous les groupes (le $CMR$), d’où $n - I$ degrés de liberté.

Mais avec $I$ groupes, il y a $m = I(I-1)/2$ comparaisons. Si chacune a un risque de 5 %, le risque de trouver au moins une fausse différence grimpe vite. La correction de Bonferroni fait chaque test au niveau $\alpha/m$ (ou, ce qui revient au même, multiplie chaque p-valeur par $m$). Ainsi, le risque global reste inférieur à $\alpha$.

Exemple

Trois comparaisons, donc $m = 3$ :

Comparaison	Écart (min)	$t$	p-valeur brute	p-valeur Bonferroni
Bus / vélo	11,33	4,90	0,0003	0,0009
Bus / voiture	6,00	2,49	0,027	0,082
Vélo / voiture	$-5{,}33$	$-2{,}31$	0,038	0,114

Après correction, seule la différence bus / vélo reste significative au niveau 5 %. Sans correction, on aurait conclu que les trois modes diffèrent tous, avec un risque réel d’erreur bien supérieur à 5 %.

Comparer chaque p-valeur brute à 5 % / 3 revient à multiplier les p-valeurs par 3 : seule la comparaison bus / vélo passe sous la ligne de Bonferroni.

library(emmeans)
emmeans(modele, pairwise ~ mode, adjust = "bonferroni")
# ou, sans paquet supplémentaire :
pairwise.t.test(trajets$temps, trajets$mode, p.adjust.method = "bonferroni")

En TP, nous avons suivi cette démarche sur une étude de consommateurs comparant des hot-dogs selon leur composition (bœuf, mélange de viandes, volaille ; une quatrième modalité, le soja, était retirée de l’analyse) : ANOVA sur les calories puis sur le sodium, vérification des hypothèses sur les résidus, puis comparaisons deux à deux avec emmeans et la correction de Bonferroni pour identifier le type le moins calorique.

Piège

Faire tous les tests de Student deux à deux sans correction, c’est le problème des tests multiples. Avec 20 comparaisons au niveau 5 %, on s’attend à une fausse découverte en moyenne. Bonferroni est simple mais conservateur ; les méthodes de Tukey (comparaisons de moyennes) ou de Benjamini-Hochberg (contrôle du taux de fausses découvertes) sont moins pénalisantes.

Vérifier les hypothèses et aller plus loin

Les hypothèses

Comme en régression, on vérifie les hypothèses sur les résidus $y_{ij} - \bar y_i$ :

• même variance dans tous les groupes : des boîtes à moustaches d’étendues comparables, ou un test de Bartlett (bartlett.test, sensible à la normalité) ou de Levene ;
• normalité : droite de Henry des résidus ;
• indépendance : garantie par le plan d’expérience (des individus différents, tirés au hasard).

L’ANOVA est assez robuste à un léger écart à la normalité, surtout avec des groupes de tailles égales. Si les variances diffèrent beaucoup, on peut utiliser la version de Welch (oneway.test en R).

L’ANOVA à deux facteurs

On peut croiser deux facteurs. Le modèle additif (sans interaction) s’écrit

$$Y_{ijk} = \mu + \alpha_i + \beta_j + \varepsilon_{ijk}$$

où $\alpha_i$ est l’effet de la modalité $i$ du premier facteur et $\beta_j$ celui de la modalité $j$ du second. La table d’ANOVA contient alors une ligne par facteur, chacune avec son test de Fisher.

L’intérêt est double. D’abord, étudier deux effets à la fois. Ensuite, retirer du bruit : en TP, des juges notaient l’acidité de trois cafés ; certains juges notent systématiquement plus sévèrement que d’autres. Intégrer le facteur « juge » (on parle de facteur bloc) retire cette variabilité du résidu, ce qui rend le test sur l’effet « café » plus puissant. En R, il suffit d’écrire lm(note ~ cafe + juge).

Si l’effet d’un facteur dépend du niveau de l’autre (un café paraît plus acide pour certains juges seulement), on ajoute un terme d’interaction $\gamma_{ij}$, écrit cafe * juge en R. Il faut alors plusieurs observations par combinaison de modalités pour pouvoir l’estimer.

En entretien

« L’ANOVA compare plusieurs moyennes en un seul test, en rapportant la variabilité entre groupes à la variabilité dans les groupes. Si elle est significative, je fais des comparaisons deux à deux avec une correction pour tests multiples, comme Bonferroni. »

En résumé

• ANOVA à un facteur : $Y_{ij} = \mu + \alpha_i + \varepsilon_{ij}$, avec une contrainte d’identifiabilité ($\sum n_i\alpha_i = 0$, ou premier groupe en référence comme dans R).
• $SST = SSM + SSR$, avec $n - 1 = (I - 1) + (n - I)$ degrés de liberté ; $R^2 = SSM/SST$.
• Test de Fisher : $F_0 = \frac{SSM/(I-1)}{SSR/(n-I)}$, à comparer à $\mathcal{F}_{I-1,n-I}$ ; la table d’ANOVA résume le tout.
• Après un rejet, comparaisons deux à deux avec le $CMR$ et une correction de Bonferroni ($\alpha/m$).
• Hypothèses : indépendance, normalité, égalité des variances ; à deux facteurs, un facteur bloc réduit le bruit et l’interaction capte les effets croisés.

Exercices

Exercice · Compléter une table d'ANOVA

facile

On compare 4 méthodes d’apprentissage sur 24 élèves (6 par méthode). On obtient $SST = 300$ et $SSR = 180$. Complétez la table d’ANOVA et concluez au niveau 5 % ($f_{3,20;0,95} \approx 3{,}10$).

Voir la solution

$SSM = 300 - 180 = 120$. Degrés de liberté : $I - 1 = 3$ et $n - I = 20$.

Source	ddl	SC	CM	$F$
Méthode	3	120	40	4,44
Résidus	20	180	9

$F_0 = 40/9 \approx 4{,}44 > 3{,}10$ : on rejette $H_0$, la méthode a un effet sur les résultats moyens. La méthode explique $R^2 = 120/300 = 40\,\%$ de la variabilité.

En code

Python

from scipy import stats

I, n, SST, SSR = 4, 24, 300, 180
SSM = SST - SSR                                    # 120
ddl_m, ddl_r = I - 1, n - I                        # 3 et 20
CMM, CMR = SSM / ddl_m, SSR / ddl_r                # 40 et 9
F0 = CMM / CMR
print(F0, stats.f.ppf(0.95, ddl_m, ddl_r))         # 4.44 > 3.10 : on rejette H0
print(stats.f.sf(F0, ddl_m, ddl_r), SSM / SST)     # p = 0.015, R² = 0.40

R

I <- 4; n <- 24; SST <- 300; SSR <- 180
SSM <- SST - SSR                            # 120
ddl_m <- I - 1; ddl_r <- n - I              # 3 et 20
CMM <- SSM / ddl_m; CMR <- SSR / ddl_r      # 40 et 9
F0 <- CMM / CMR
F0; qf(0.95, ddl_m, ddl_r)                  # 4.44 > 3.10 : on rejette H0
1 - pf(F0, ddl_m, ddl_r); SSM / SST         # p = 0.015, R² = 0.40

Exercice · Interpréter des coefficients de R

moyen

Une ANOVA du rendement d’une culture selon l’engrais (modalités A, B, C) donne avec summary(lm(rendement ~ engrais)) : constante 52,0 ; engraisB 4,5 ; engraisC $-2{,}0$. Les groupes ont tous 8 parcelles.

1. Quelles sont les moyennes des trois groupes ?
2. Quelle est la moyenne générale, et quels seraient les effets $\alpha_i$ avec la contrainte naturelle $\sum n_i \alpha_i = 0$ ?

Voir la solution

1. R utilise A comme référence : $\bar y_A = 52{,}0$, $\bar y_B = 52{,}0 + 4{,}5 = 56{,}5$ et $\bar y_C = 52{,}0 - 2{,}0 = 50{,}0$.
2. Les effectifs étant égaux, la moyenne générale vaut $(52{,}0 + 56{,}5 + 50{,}0)/3 = 52{,}83$. Avec la contrainte naturelle, $\mu = 52{,}83$ et les effets sont les écarts à cette moyenne : $\alpha_A \approx -0{,}83$, $\alpha_B \approx 3{,}67$ et $\alpha_C \approx -2{,}83$, dont la somme (pondérée par les effectifs égaux) est bien nulle.

En code

Python

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

# Coefficients lus dans la sortie (A en référence)
constante, effet_B, effet_C = 52.0, 4.5, -2.0
moyennes = np.array([constante, constante + effet_B, constante + effet_C])
print(moyennes)                                    # [52.  56.5 50. ]
mu = moyennes.mean()                               # effectifs égaux
print(mu, moyennes - mu)                           # 52.83  [-0.83  3.67 -2.83]

# Vérification : 8 parcelles par engrais, construites avec ces moyennes
ecarts = np.array([-3, -2, -1, 0, 0, 1, 2, 3])     # écarts centrés
df = pd.DataFrame({"engrais": np.repeat(["A", "B", "C"], 8),
                   "rendement": np.concatenate([m + ecarts for m in moyennes])})
print(smf.ols("rendement ~ engrais", df).fit().params)            # 52.0  4.5  -2.0
print(smf.ols("rendement ~ C(engrais, Sum)", df).fit().params)    # 52.83  -0.83  3.67

R

# Coefficients lus dans la sortie (A en référence)
constante <- 52.0; effet_B <- 4.5; effet_C <- -2.0
moyennes <- c(A = constante, B = constante + effet_B, C = constante + effet_C)
moyennes                                    # 52.0 56.5 50.0
mu <- mean(moyennes)                        # effectifs égaux
mu; moyennes - mu                           # 52.83  -0.83 3.67 -2.83

# Vérification : 8 parcelles par engrais, construites avec ces moyennes
ecarts <- c(-3, -2, -1, 0, 0, 1, 2, 3)      # écarts centrés
df <- data.frame(engrais   = factor(rep(c("A", "B", "C"), each = 8)),
                 rendement = c(outer(ecarts, moyennes, "+")))
coef(lm(rendement ~ engrais, data = df))    # 52.0  4.5  -2.0 (A en référence)
# Contrainte somme nulle : constante = mu, puis alpha_A et alpha_B
coef(lm(rendement ~ engrais, data = df, contrasts = list(engrais = "contr.sum")))
## 52.83  -0.83  3.67 (alpha_C = -(alpha_A + alpha_B) = -2.83)

Exercice · Combien de comparaisons ?

difficile

Une ANOVA sur 5 groupes est significative. On veut comparer tous les groupes deux à deux avec un risque global d’au plus 5 %.

1. Combien de comparaisons faut-il faire, et à quel niveau individuel avec Bonferroni ?
2. Si l’on faisait ces comparaisons à 5 % chacune et qu’elles étaient indépendantes, quelle serait la probabilité d’au moins une fausse découverte lorsque tous les groupes ont la même moyenne ?

Voir la solution

1. $m = 5 \times 4/2 = 10$ comparaisons, chacune au niveau $0{,}05/10 = 0{,}005$.
2. La probabilité de n’avoir aucune fausse découverte serait $0{,}95^{10} \approx 0{,}599$, donc celle d’en avoir au moins une vaut environ $1 - 0{,}599 = 0{,}40$. Sans correction, on aurait quatre chances sur dix de « découvrir » une différence inexistante. Les comparaisons ne sont pas vraiment indépendantes (elles partagent des groupes), mais l’ordre de grandeur montre bien le danger.

En code

Python

from itertools import combinations
import numpy as np
from scipy import stats

I, alpha = 5, 0.05
m = I * (I - 1) // 2
print(m, alpha / m, 1 - (1 - alpha)**m)            # 10  0.005  0.401

# Simulation : 5 groupes de même moyenne, 10 tests de Student à 5 % chacun
rng = np.random.default_rng(0)
R, sans_corr, bonf = 5000, 0, 0
for _ in range(R):
    groupes = rng.normal(0, 1, size=(I, 10))
    p = [stats.ttest_ind(groupes[i], groupes[j]).pvalue
         for i, j in combinations(range(I), 2)]
    sans_corr += min(p) < alpha                     # au moins une fausse découverte
    bonf += min(p) < alpha / m
print(sans_corr / R, bonf / R)                      # 0.29  0.043 : moins que 0.40 (tests liés), Bonferroni sous 5 %

R

I <- 5; alpha <- 0.05
m <- choose(I, 2)
c(m, alpha / m, 1 - (1 - alpha)^m)          # 10  0.005  0.401

# Simulation : 5 groupes de même moyenne, 10 tests de Student à 5 % chacun
set.seed(0)
groupe <- factor(rep(1:I, each = 10))
p_min <- replicate(5000, {
  y <- rnorm(50)
  p <- pairwise.t.test(y, groupe, pool.sd = FALSE, var.equal = TRUE,
                       p.adjust.method = "none")$p.value
  min(p, na.rm = TRUE)
})
mean(p_min < alpha); mean(p_min < alpha / m)   # 0.29  0.037 : moins que 0.40 (tests liés), Bonferroni sous 5 %