UNIVERSITÉ PARIS OUEST NANTERRE LA DÉFENSE

U.F.R. SEGMI Année universitaire 2012 – 2013

Master d’économie Cours de M. Desgraupes

MATHS/STATS

Document 6 : Exemple d’ANOVA

1 Analyse de variance 1
1.1 Énoncé du problème . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Calcul à la main . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.1 Calcul des inerties . . . . . . . . . . . . . . . . . . . . . . 3
1.2.2 Statistique de Fisher . . . . . . . . . . . . . . . . . . . . . 4
1.3 Calcul par régression . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4 Table d’analyse de la variance . . . . . . . . . . . . . . . . . . . . 7
1.5 Test HSD de Tukey . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.6 Test de Kruskal-Wallis . . . . . . . . . . . . . . . . . . . . . . . . 10

1 Analyse de variance
1.1 Énoncé du problème
On considère les données suivantes représentant des taux de créatinine rele-
vés dans quatre groupes constitués de n = 5 cobayes. Le premier groupe est un
groupe témoin, les trois autres ont subi l’injection d’un produit pharmaceutique
à trois doses différentes.
Les résultats se présentent sous la forme d’un tableau de dimension 5 × 4
comme ceci :
A B C D
11.7 11.9 13.4 12.8
12.1 12.5 12.8 12.8
11.6 12.7 12.3 13.3
12.8 12.6 11.5 13.2
12.2 12.1 13.2 13.1
Pour stocker les données dans R, on va définir un vecteur X contenant toutes
les données (parcourues en colonne) :

> X <- c(11.7, 12.1, 11.6, 12.8, 12.2, 11.9, 12.5, 12.7, 12.6, 12.1,
+ 13.4, 12.8, 12.3, 11.5, 13.2, 12.8, 12.8, 13.3, 13.2, 13.1)

1
 Il sera utile aussi de les utiliser sous forme d’une matrice x. On transforme
le vecteur X en une matrice comme ceci :
> n <- 5
> k <- 4
> x <- matrix(X, nrow=n, ncol=k)
> colnames(x) <- LETTERS[1:4]
On calcule les moyennes dans chaque groupe, c’est-à-dire dans chaque co-
lonne de la matrice, au moyen de la fonction colMeans :
> moy <- colMeans(x)

A B C D
12.08 12.36 12.64 13.04

On cherche à déterminer si le résultats obtenus laissent apparaı̂tre une dif-

férence significative entre les quatre groupes.
On peut commencer par visualiser les données de chaque groupe au moyen
de boı̂tes à moustaches.
13.0
12.5
12.0
11.5

1 2 3 4

2
1.2 Calcul à la main
1.2.1 Calcul des inerties
On commence par calculer la somme des inerties intra-groupes.
> intra <- apply(x, 2, function(y) sum((y-mean(y))^2))
Leur somme, notée SCintra, vaut 3.924 :
> SCintra <- sum(intra)

[1] 3.924

L’inertie inter-groupes est l’inertie entre le barycentre de tous les points et

les barycentres de chacun des 4 groupes.
Le barycentre de tous les points est :

> g <- mean(X)

[1] 12.53

Donc, l’inertie inter-groupes est obtenue comme ceci :

> SCinter <- n*sum( (moy - g)^2 )

[1] 2.518

L’inertie totale est la somme des carrés des écarts avec le barycentre général :
> SCtotal <- sum( (X-g)^2 )

[1] 6.442

> SCtotal

[1] 6.442

On vérifie qu’elle est la somme des inerties intra et de l’inertie inter :

> SCintra + SCinter

[1] 6.442

Une autre manière d’interpréter les inerties inter et intra est de considérer
que l’inertie inter-groupe mesure des effets et que l’inertie intra-groupe est une
mesure des erreurs.

3
1.2.2 Statistique de Fisher
Le test de comparaison de toutes les moyennes vise à tester si la moyenne
peut être considérée comme étant la même dans chaque groupe. Ce test est
valide dans les conditions suivantes :
– les k groupes sont indépendants et tirés au hasard de leurs populations
respectives ;
– les populations ont une distribution normale ;
– la variance est la même pour chaque groupe.
Lorsque les groupes ont exactement le même nombre d’individus, on dit que
le test est robuste et reste valable si on s’écarte un peu des conditions énoncées.
En revanche, s’ils n’ont pas même effectif, il faut s’assurer que les conditions
ci-dessus sont remplies.
Le test de Fisher pose l’hypothèse H0 suivante :
H0 : les moyennes de tous les groupes sont égales entre elles
Le nombre de degrés de liberté pour l’inertie inter-groupes est égal au nombre
de groupes diminué de 1 car on a estimé un paramètre (à savoir le barycentre
de l’ensemble de tous les individus) :
> dlinter <- 4-1

[1] 3

Le nombre de degrés de liberté pour l’inertie intra-groupes est égal à n − 1

pour chaque groupe car on a estimé le barycentre de chacun, ce qui fait au total
4 × (n − 1) = 16 :
> dlintra <- 4*(n-1)

[1] 16

La statistique de Fisher est la quantité (SCinter/dlinter)/(SCintra/dlintra).

On trouve :
> f <- (SCinter/dlinter)/(SCintra/dlintra)

[1] 3.422358

Si l’hypothèse H0 est vraie, cette variable suit une loi de Fisher à (3, 16)
degrés de liberté. La valeur critique au seuil 5% est :
> qf(0.95,dlinter,dlintra)

[1] 3.238872

Comme la valeur calculée 3.422358 est supérieure à la valeur critique 3.238872,

on rejette l’hypothèse avec un risque d’erreur de 5% de se tromper.
On pourrait aussi calculer la p-valeur correspondant à f :

4
> pval <- 1-pf(f,dlinter,dlintra)

[1] 0.04275917

On voit qu’elle est inférieure à 5% = 0.05 et cela confirme qu’on peut rejeter
l’hypothèse H0 .
Le coefficient de détermination R2 est la part d’inertie expliquée (inertie
inter-groupes) par rapport à l’inertie totale :
> R2 <- SCinter/SCtotal

[1] 0.3908724

Il est effectivement assez faible et indique que la seule appartenance à un

groupe particulier n’explique que faiblement les différences constatées. La part
expliquée de la variance totale est de 39.09%.

1.3 Calcul par régression

Rappelons qu’on a réuni, dans la section précédente, tous les groupes en un
unique vecteur appelé X.
Parallèlement, on fabrique un facteur (variable qualitative) appelé G qui
indique l’appartenance de chaque individu à un des 4 groupes :
> G <- factor(rep(1:4,each=5))

[1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4
Levels: 1 2 3 4

On effectue une régression pour expliquer la variable quantitative X au

moyen de la variable qualitative G. La régression est faite au moyen de la fonc-
tion lm :
> reg <- lm(X~G)

Call:
lm(formula = X ~ G)

Coefficients:
(Intercept) G2 G3 G4
12.08 0.28 0.56 0.96

La fonction summary permet d’afficher les principaux résultats de la régres-

sion avec les tests de significativité sur les coefficients :
> summary(reg)

5
Call:
lm(formula = X ~ G)

Residuals:
Min 1Q Median 3Q Max
-1.140 -0.280 0.090 0.245 0.760

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.0800 0.2215 54.544 <2e-16 ***
G2 0.2800 0.3132 0.894 0.3846
G3 0.5600 0.3132 1.788 0.0927 .
G4 0.9600 0.3132 3.065 0.0074 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4952 on 16 degrees of freedom

Multiple R-squared: 0.3909, Adjusted R-squared: 0.2767
F-statistic: 3.422 on 3 and 16 DF, p-value: 0.04276

On voit que la moyenne du premier groupe (qui joue le rôle de groupe témoin)
est très significative : elle est marquée de trois astérisques. Les coefficients trou-
vés pour les autres groupes sont la différence entre la moyenne de ces groupes et
celle du premier groupe. La moyenne du quatrième groupe est aussi considérée
comme significative (avec deux astérisques).
> coef <- coefficients(reg)

(Intercept) G2 G3 G4
12.08 0.28 0.56 0.96

Par exemple, la moyenne dans le deuxième groupe vaut 12.36 comme on l’a
vu dans la première section. On vérifie effectivement que 12.08 + 0.28 = 12.36.
La manière dont sont calculés les coefficients s’appelle un codage par contraste.
Le contraste utilisé ici – qui consiste à prendre un des groupes comme référence
et à calculer les coefficients des autres groupes par différence avec le coefficient
du groupe témoin – est un type de codage qui s’appelle contraste de traitement.
On peut vérifer que c’est bien le codage utilisé comme ceci :
> reg$contrast

$G
[1] "contr.treatment"

La dernière ligne du résumé contient les informations relatives au test de

Fisher. On y retrouve les valeurs calculées à la section précédente : la valeur de
la statistique égale à 3.422, les degrés de liberté 3 et 16 ainsi que la p-valeur.
L’avant-dernière ligne comporte le coefficient R2 égal à 0.391.

6
 Pour mieux comprendre la représentation des variables qualitatives dans un
modèle linéaire, on peut afficher la matrice du modèle :
> model.matrix(reg)

(Intercept) G2 G3 G4
1 1 0 0 0
2 1 0 0 0
3 1 0 0 0
4 1 0 0 0
5 1 0 0 0
6 1 1 0 0
7 1 1 0 0
8 1 1 0 0
9 1 1 0 0
10 1 1 0 0
11 1 0 1 0
12 1 0 1 0
13 1 0 1 0
14 1 0 1 0
15 1 0 1 0
16 1 0 0 1
17 1 0 0 1
18 1 0 0 1
19 1 0 0 1
20 1 0 0 1
attr(,"assign")
[1] 0 1 1 1
attr(,"contrasts")
attr(,"contrasts")$G
[1] "contr.treatment"

1.4 Table d’analyse de la variance

Il est traditionnel de rassembler les résultats des calculs d’inertie inter et
intra dans une table appelée table d’analyse de variance. Dans le cas présent,
cette table se présente sous la forme suivante :
d.l. Inertie Variance Stat f p-valeur
G 3 2.518 0.839 3.422 0.043
Résidus 16 3.924 0.245

La première colonne (appelée d.l.) contient les nombres de degrés de liberté.

La deuxième colonne comporte les inerties inter et intra respectivement. La
troisième colonne est obtenue en faisant le quotient de la deuxième colonne par
la première : ce sont les variances, c’est-à-dire les moyennes des inerties. La
quatrième colonne comporte la statistique qui est le quotient des deux valeurs

7
trouvées dans la troisième colonne. Enfin la dernière colonne indique la p-valeur
correspondant à la statistique f trouvée dans la quatrième colonne.
Il existe une fonction anova qui construit la table d’analyse de variance à
partir de l’objet produit par la régression linéaire :
> anv <- anova(reg)

Analysis of Variance Table

Response: X
Df Sum Sq Mean Sq F value Pr(>F)
G 3 2.518 0.83933 3.4224 0.04276 *
Residuals 16 3.924 0.24525
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

On retrouve évidemment toutes les valeurs calculées antérieurement et le fait

que la p-valeur est inférieure à 5%, ce qui conduit à rejeter l’hypothèse H0 du
test de Fisher.
On peut aussi utiliser la fonction aov qui a une syntaxe un peu différente. Elle
ne prend pas en argument le résultat de la régression mais directement la formule
qui a permis de définir le modèle de régression, à savoir l’expression X~G utilisée
au début de la section 1.3 dans la fonction lm. La fonction aov se charge elle-
même d’effectuer la régression et de bâtir un résumé d’analyse de variance. Par
exemple, dans le cas de notre exemple, il faudrait écrire l’instruction suivante :
> av <- aov(X~G)

Call:
aov(formula = X ~ G)

Terms:
G Residuals
Sum of Squares 2.518 3.924
Deg. of Freedom 3 16

Residual standard error: 0.4952272

Estimated effects may be unbalanced

L’objet renvoyé par cette fonction est de classe aov et comporte les mêmes
composantes que l’objet reg de classe lm.

1.5 Test HSD de Tukey

Le test de comparaisons multiples de Tukey, dit test HSD (abréviation de
Honestly Significant Difference) a pour but de distinguer parmi les échantillons
s’il y en a qui diffèrent significativement des autres. Dans le cas où le test

8
de Fisher a conduit à rejeter l’hypothèse que les échantillons de diffèrent pas
significativement, on cherche à savoir quels échantillons se distinguent des autres.
Pour cela, le test HSD envisage les échantilons deux par deux et calcule dans
chaque cas la statistique suivante :

|Mk − Mkj |
Qij = p i
M intra/n

où Mintra est la moyenne des inerties intra-groupes (donc le quotient de l’inertie
intra par le nombre de degrés de liberté intra) et n est le nombre d’observations
dans chaque échantillon.
Lorsque les échantillons n’ont pas la même taille on prend pour n la moyenne
harmonique des effectifs dans chaque échantillon.
Si il y a k échantillons, cela fait k(k − 1)/2 paires.
La quantité au dénominateur de Qij ne dépend pas des indices i et j.
Il existe des valeurs critiques pour les quantités Q en fonction du nombre k
d’échantillons et du nombre de degrés de liberté intra. On obtient ces valeurs
critiques avec la fonction qtukey. Par exemple, si on est à 5%, la valeur critique
pour 4 échantillons et 16 degrés de liberté est approximativement 4.05 :
> Qc <- qtukey(.95, 4, 16)

[1] 4.046093

Calculons par exemple la quantité Q14 correspondant à la différence de

moyennes entre les échantillons 1 et 4. On a
> Q <- (abs(moy[1]-moy[4]))/sqrt(SCintra/(dlintra*n))

A
4.334627

La quantité obtenue est à comparer à la valeur critique 4.05 et montre,

puisqu’elle est supérieure, que la différence entre les échantillons 1 et 4 est si-
gnificative.
La quantité au dénominateur de la statistique Q vaut ici
> sqrt(SCintra/(dlintra*n))

[1] 0.2214723

Si on la multiplie par la valeur critique 4.05, on obtient le rayon de l’intervalle

de confiance autour des différences de moyennes, à savoir
> dr <- Qc * sqrt(SCintra/(dlintra*n))

[1] 0.8960977

9
 Par exemple, la différence de moyennes entre les échantillons 1 et 4 est 0.96.
L’intervalle de confiance autour de cette valeur est donc

[0.96 − 0.8961, 0.96 + 0.8961] = [0.0639, 1.8561].

La fonction TukeyHSD peut être utilisée pour effectuer le test HSD de com-
paraisons multiples. Elle prend en argument l’objet de classe aov renvoyé par la
fonction aov. Par exemple :
> tuk <- TukeyHSD(av)

Tukey multiple comparisons of means

95% family-wise confidence level

Fit: aov(formula = X ~ G)

$G
diff lwr upr p adj
2-1 0.28 -0.61609772 1.176098 0.8080588
3-1 0.56 -0.33609772 1.456098 0.3145281
4-1 0.96 0.06390228 1.856098 0.0336561
3-2 0.28 -0.61609772 1.176098 0.8080588
4-2 0.68 -0.21609772 1.576098 0.1736154
4-3 0.40 -0.49609772 1.296098 0.5896220

La colonne diff donne les différences entre les moyennes observées. Les co-
lonnes lwr et upr donnent les bornes inférieure et supérieure de l’intervalle et
la colonne p adj donne la p-valeur après ajustement pour les comparaisons
multiples. Si la valeur 0 n’est pas dans un intervalle, c’est que la différence
correspondante peut être considérée comme significative.
Il apparaı̂t que la comparaison entre les échantillons 1 et 4 présente une
différence significative : la p-value 0.03366 est inférieure à 5%, ce qui implique
qu’on rejette l’hypothèse que les moyennes sont égales.
La figure suivante donne une représentation graphique de ces intervalles ob-
tenue au moyen de la fonction plot.
On voit clairement que la valeur 0 n’est pas comprise dans l’intervalle cor-
respondant à 1-4.

1.6 Test de Kruskal-Wallis

Lorsque les hypothèses du test de Fisher énoncées dans la section 1.2.2 ne
sont pas remplies, on a recours à un test non-paramétrique appelé le test de
Kruskal-Wallis. Ce test repose sur la façon dont sont ordonnées les valeurs consti-
tuant les groupes et est parfois qualifié d’anova sur les rangs. Il ne s’agit pas à
proprement parler d’une analyse de variance mais la démarche est très similaire.
Supposons qu’on dispose des données suivantes correspondant à trois groupes
de tailles différentes dans le contexte de l’exemple introduit dans la section 1.1 :

10
> plot(tuk)

95% family−wise confidence level

2−1
3−1
4−1
3−2
4−2
4−3

−0.5 0.0 0.5 1.0 1.5

Differences in mean levels of G

A B C
11.7 12.1 12.7
12.1 13.15 13.2
11.6 12.5 13.1
12.8 11.9 13.05
12.2 11.2 13.3
11.6 12.9
12.0 13.0
12.4
12.3
On peut stocker les données dans R sous forme d’une liste :
> k <- 3
> L <- list(
+ A=c(11.7, 12.1, 11.6, 12.8, 12.2, 11.6, 12.0, 12.4, 12.3),
+ B=c(12.1, 13.15, 12.5, 11.9, 11.2, 12.9, 13.0),
+ C=c(12.7, 13.2, 13.1, 13.05, 13.3)
+ )
On calcule les tailles et les moyennes de chaque groupe comme ceci :

11
> len <- sapply(L,length)
A B C
9 7 5
> N <- sum(len)
[1] 21
> moy <- sapply(L,mean)
A B C
12.07778 12.39286 13.07000
Le problème est de déterminer si ces moyennes diffèrent ou pas. L’hypothèse
H0 est la même que pour une anova traditionnelle :
H0 : les moyennes sont égales entre elles
Le principe est de remplacer les valeurs du tableau précédent par leur rang.
On fusionne donc toutes les valeurs et on les remplace par leur rang, c’est-à-dire
la position qu’elles occupent lorsqu’on les trie dans l’ordre croissant. Lorsqu’il
y a des valeurs ex-aequo, on remplace leur rang par la moyenne des rangs cor-
respondants. On exécute les instructions suivantes :
> Y <- c(L[[1]],L[[2]],L[[3]])
[1] 11.70 12.10 11.60 12.80 12.20 11.60 12.00 12.40 12.30 12.10 13.15 12.50
[13] 11.90 11.20 12.90 13.00 12.70 13.20 13.10 13.05 13.30
> sort(Y)
[1] 11.20 11.60 11.60 11.70 11.90 12.00 12.10 12.10 12.20 12.30 12.40 12.50
[13] 12.70 12.80 12.90 13.00 13.05 13.10 13.15 13.20 13.30
> rk <- rank(Y)
[1] 4.0 7.5 2.5 14.0 9.0 2.5 6.0 11.0 10.0 7.5 19.0 12.0 5.0 1.0 15.0
[16] 16.0 13.0 20.0 18.0 17.0 21.0
Par exemple, ici on voit que la valeur 11,6 est présente deux fois, aux rangs 2 et
3 (lorsque Y est trié). On remplace donc ces rangs par leur moyenne 2,5. Avec
R, la fonction rank effectue automatiquement ces corrections pour les ex-aequos.
On obtient maintenant le tableau des rangs :
A B C
4 7.5 13
7.5 19 20
2.5 12 18
14 5 17
9 1 21
2.5 15
6 16
11
10

12
 On va calculer l’équivalent d’une intertie inter-groupes sur les rangs. On la
note SCrg . On peut l’obtenir au moyen de deux formules différentes. La première
utilise les moyennes :
Xk
SCrg = nj (Ȳj − Ȳ )2
j=1

où Ȳj désigne la moyenne des rangs dans le groupe j et Ȳ la moyenne totale des
rangs. L’autre formule est :
k
X Sj2 S2
SCrg = −
j=1
nj N

où Sj désigne la somme des rangs dans le groupe j et S la somme totale des
rangs.
Cette inertie inter-groupes permet de définir la statistique de Kruskal-Wallis
comme ceci :
SCrg
H= (1)
N (N + 1)/12
Sous l’hypothèse H0 , cette statistique suit une loi du χ2 à k − 1 degrés de
liberté, où k est le nombre de groupes.
On peut mener les calculs avec R de la manière suivante en commençant à
calculer les sommes des rangs dans chaque groupe :
> idx <- 0
> Srk <- vector(length=k);
> for (i in 1:k) {
+ Srk[i] <- sum(rk[(idx+1):(idx+len[i])])
+ idx <- idx+len[i]
+ }
On obtient alors la statistique de Kruskal-Wallis :
> SCrg <- sum( Srk^2/len ) - sum(rk)^2/length(rk)

[1] 348.8825

> H <- 12*SCrg/(N*(N+1))

[1] 9.061884

La borne critique du χ2 à 2 degrés de liberté est :

> qchisq(0.975,k-1)

[1] 7.377759

Puisque la valeur calculée 9.061884 est supérieure à la valeur critique 7.377759,

on rejette l’hypothèse et on considère que les moyennes dans ces trois groupes
diffèrent significativement les unes des autres, au risque 5% de se tromper.

13
Remarque : il existe une fonction kruskal.test définie dans le package stats
de R qui effectue ce test automatiquement. Il faut définir un facteur permettant
de distinguer les groupes et passer en arguments le vecteur de toutes les va-
leurs ainsi que ce facteur. On retrouve approximativement les valeurs calculées
précédemment :

> g <- rep(1:3,len)

[1] 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 3 3 3 3 3

> kw <- kruskal.test(Y,g)

Kruskal-Wallis rank sum test

data: Y and g
Kruskal-Wallis chi-squared = 9.0737, df = 2, p-value = 0.01071

La p-valeur 0.01071 est inférieure à 5% et on rejette l’hypothèse H0 .

14