Tutoriel : comparaison de deux moyennes avec le logiciel R

Alors vous y voilà, vous avez recueilli deux échantillons d’observations numériques continues, et vous vous demandez si les données de ces deux groupes sont différentes ou non. Autrement formulé vous vous demandez si les populations dont sont issus vos deux échantillons sont différentes ou pas.

 

Ce type d’approche appartient aux analyses inférentielles : les conclusions tirées de l’étude des échantillons seront généralisées à l’échelle des populations.

 

Pour comparer deux échantillons d’observations il est classique d’utiliser un paramètre de tendance centrale, et généralement il s’agit de la moyenne ,voir cet article pour plus de détails.

 

Deux tests statistiques, le test de Student et le test de Wilcoxon, sont généralement employés pour comparer deux moyennes. Il existe cependant des variantes de ces deux tests, pour répondre à différentes situations, comme la non indépendance des échantillons par exemple. Ces différentes versions des tests correspondent à des options dans les arguments des fonctions R. Il est nécessaire de les connaitre et de les maîtriser pour réaliser des analyses adéquates aux données et à la question posée.

 

Dans ce tutoriel, je vais donc vous montrer, en pas à pas, la démarche que j’utilise pour comparer deux moyennes avec le logiciel R, dans un maximum de situations.

 

Pour commencer, je vais faire quelques brefs rappels sur les hypothèses et les principes des tests de Student et de Wilcoxon

1. Rappels

1.1 Les Hypothèses des tests de Student et de Wilcoxon

Les hypothèses sont :

  • H0 : µ1 = µ2
  • H1 : µ1 ≠ µ2 (bilatérale) ou H1 : µ1 < µ2 (unilatérale) ou H1 :µ1 > µ2 (unilatérale)

µ1 et µ2 sont les moyennes des populations dont sont issues les échantillons de données.

 

Certains auteurs disent que les hypothèses du test de Wilcoxon concernent les médianes et pas les moyennes, puisque c’est un test basé sur les rangs. D’autres encore comme (Conover 1999), affirment que les hypothèses concernent les distributions, mais que le test peut être utilisé pour comparer des moyennes.

 

1.2 Principe

Quel que soit le test employé, le principe consiste toujours à calculer une statistique (c’est une quantité calculée selon une équation qui dépend du test employé).

 

Par exemple, la statistique d’une des forme du test de Student est :

test de Student

 

Le logiciel va ensuite comparer la statistique calculée avec les valeurs théoriques attendues sous l’hypothèse nulle ** à partir d’une loi de distribution théorique ou établie par tests de permutation. Et il va de cette façon déterminer **la probabilité de la statistique calculée sous l’hypothèse nulle. C’est la p-value.

 

La p-value est donc la probabilité, sous l’hypothèse nulle que les moyennes ne sont pas différentes, d’observer la valeur de statistique qui a été observée 😉 *

 

Le seuil de significativité, ou risque alpha est très très généralement fixé à 0.05.

 

Si la pvalue est < 0.05 alors l’hypothèse nulle (H0) est rejetée.Et dans ce cas on va conclure à la différence entre les deux moyennes des populations dont sont issus les échantillons.

 

En revanche, si la pvalue est > 0.05 alors l’hypothèse nulle n’est pas rejetée. Et on va conclure à l’absence de différence entre les deux populations

 

1.3 Le test de Student

Il s’agit d’un test paramétrique, c’est à dire qui se base sur une distribution statistique théorique (ici la distribution de Student), qui dépend de paramètres (ici le degrés de liberté = n1 + n2 -1). n1 et n2 sont les tailles des échantillons.

 

Il existe plusieurs statistique du test de Student selon que les variances des échantillons sont différentes ou pas, ou encore selon que les échantillons sont indépendants ou appariés. Pour plus de détails vous pouvez consulter à peu près n’importe quel livre de biostatistiques ou simplement la page de wikipedia du test de Student en anglais.

 

Lorsque les échantillons sont indépendants, le test de Student nécessite que les observations de chacun d’eux suivent une distribution Normale.

 

Lorsque les échantillons sont appariés c’est la différence entre les observations qui doit suivre une loi Normale.

 

Lorsque les échantillons sont indépendants, l’homogénéité des variances (cad leur égalité) doit être testée par un test F. En cas de rejet de l’hypothèse, une variante du test de Student peut être employée : le test de Welch.

 

Le test de Student est censé être plus puissant que le test de Wilcoxon. Cependant, il est sensible aux outliers et à l’assymétrie des données qui peuvent entraîner une augmentation invisible du risque alpha. Ceci va se traduire, par exemple, par une p-value anormalement
basse.
 

1.4 Le test de Wilcoxon

Il s’agit d’un test non paramétrique, qui est basé sur le calcul desommes de rangs. Sa statistique est décrite sur cette page wikipedia.

 

Le test de Wilcoxon ne nécessite aucune condition d’utilisation. C’est ce qui en fait un test particulièrement robuste. En revanche, il est censé être moins puissant que le test de Student. Cela est sans doute vrai lorsque les données suivent parfaitement une distribution normale mais comme cela est rarement le cas, l’avantage du test de Student n’est sans doute pas toujours fondée.

 

2. Comment choisir le test à utiliser?

Voici la règle que j’utilise lorsque les échantillons sont indépendants :

  • si n1 et/ou n2 < 6 : test de Wilcoxon
  • si n1 et n2 > 6
    • si les hypothèses de normalité sont accpetées (par test de Shapiro Wilk) et qu’aucune asymétrie marquée, ni outlier n’ont été mis en évidence (visuellement ou par test)
      • si les variances sont égales: alors test de Student
      • si les variances sont différentes, alors test de Welch
    • dans tous les autres cas : test de Wilcoxon.

 

Lorsque les échantillons sont appariés, j’utilise les mêmes règles à l’exception du test de l’égalité des variances puisqu’une seule variance est employée, celle de la différence.

 

En cas de doute, les deux tests peuvent aussi être réalisés de façon complémentaire et un seul sera cité dans le compte rendu d’analyse ou de la publication.

 

3. Visualisation des données

La première démarche à faire est toujours de visualiser les données observées.

J’utilise généralement des boxplots, en superposant les observations, et en ajoutant la moyenne. Cette représentation permet d’évaluer rapidement une éventuelle asymétrie des données (si les bords supérieur et inférieur des boites ne sont pas symétriques, ou si la moyenne n’est pas sur la barre à l’intérieure de la boite), ainsi que la présence éventuelle d’outliers (les points qui sont au delà des traits de parts et d’autres des boites.). Voir [cet article]
(https://en.wikipedia.org/wiki/Wilcoxon_signed-rank_test) pour plus de détails sur les boxplots.

Voici un exemple :

library(tidyverse)
iris %>%
    filter(Species %in% c("setosa", "virginica")) %>%
    ggplot(aes(x=Species, y=Sepal.Length, fill=Species, colour=Species)) +
        geom_boxplot(alpha=0.25, outlier.alpha=0) +
        geom_jitter(fill="black")+
        stat_summary(fun.y=mean, colour="white", geom="point", shape=18, size=5) 

comparaison de moyennes

 

L’argument outlier.alpha=0 dans la fonction geom_boxplot permet de ne pas représenter deux fois les points outliers (1 fois par la fonction geom_boxplot et 1 fois par la fonction geom_point).

 

Ici, les boxplots montrent que les observations de l’échantillon setosa sont symétriques, sans outlier évident. En revanche, les observations de l’échantillon Virginica montrent une distribution légèrement asymétrique, avec la présence d’un outliers. Ces données ont probablement un défaut de normalité.

 

J’ai également l’habitude de représenter les densités, en ajoutant les distributions marginales avecgeom_rug sur l’axe des x. Ca permet d’avoir une première idée sur la normalité ou non des données, et de mieux voir comment les données des deux échantillons se chevauchent.

iris %>%
    filter(Species %in% c("setosa", "virginica")) %>%
    ggplot(aes(Sepal.Length, fill=Species, colour=Species)) +
    geom_density(alpha=0.25)+
    geom_rug()

comparaisons de moyennes

 

4. Evaluations préliminaires aux tests de comparaison

4.1 Normalité

La seconde chose à faire est d’évaluer la normalité des observations.

4.1.1 Cas des échantillons indépendants

Si les échantillons sont indépendants, il s’agit d’évaluer la normalité des observations des deux échantillons. On peut employer un test statistique comme le test de de Shapiro-Wilk avec la fonction Shapiro.test. L’hypothèse de normalité est acceptée si la pvalue est > 0.05.

setosa <- iris %>%
    filter(Species=="setosa") 
shapiro.test(setosa$Sepal.Length)
## 
##  Shapiro-Wilk normality test
## 
## data:  setosa$Sepal.Length
## W = 0.9777, p-value = 0.4595

Ici, l’hypothèse de normalité est acceptée.

On peut également utiliser une méthode visuelle comme les QQplot, qui ont l’avantage de mettre également en évidence les outliers. Pour satisfaire l’hypothèse de normalité les points doivent être alignés le long de la droite.

qqnorm(setosa$Sepal.Length)
qqline(setosa$Sepal.Length)

hypothèse de normalité

C’est le cas ici, l’hypothèse de normalité est acceptée.

 

virginica <- iris %>%
    filter(Species=="virginica") 
shapiro.test(virginica$Sepal.Length)
## 
##  Shapiro-Wilk normality test
## 
## data:  virginica$Sepal.Length
## W = 0.97118, p-value = 0.2583

L’hypothèse de normalité des observations de l’échantillon Virginica est accepté par le test de Shapiro Wilks.

 

qqnorm(virginica$Sepal.Length)
qqline(virginica$Sepal.Length)

hypothèse de normalité

Le QQ plot permet de mettre en évidence un défaut de normalité ainsi que la présence d’au moins un outlier.

 

4.1.2 Cas des échantillons appariés

Dans ce cas, c’est la normalité de la différence entre les deuxobservations qui doit être évaluée. On peut créer une variable “diff”” avec la fonction mutate du package dplyr.

library(MASS)
data(immer)
head(immer,5)

##   Loc Var    Y1   Y2
## 1  UF   M  81.0 80.7
## 2  UF   S 105.4 82.3
## 3  UF   V 119.7 80.4
## 4  UF   T 109.7 87.2
## 5  UF   P  98.3 84.2

immer2 <- immer %>%
    mutate(diff=Y2-Y1)

shapiro.test(immer2$diff)

## 
##  Shapiro-Wilk normality test
## 
## data:  immer2$diff
## W = 0.93785, p-value = 0.07959

qqnorm(immer2$diff)
qqline(immer2$diff)

comparaison de moyennes

 

Ici l’acceptation de l’hypothèse de normalité ne me semble pas tout à fait évidente. Dans ce genre de situation j’utilise de façon complémentaire les tests de Student et de Wilcoxon.

 

4.2 Asymétrie

L’asymétrie peut se juger visuellement avec un boxplot, commme montré lors de l’étape de visualisation. Si on souhaite une évaluation plus objective, on peut utiliser la fonction stat.desc du package pastecs, avec l’arugment norm=TRUE. Si la quantité skew.2E est > 1 alors l’asymétrie est significative.

library(pastecs)
stat.desc(setosa$Sepal.Length, norm=TRUE)

##      nbr.val     nbr.null       nbr.na          min          max 
##  50.00000000   0.00000000   0.00000000   4.30000000   5.80000000 
##        range          sum       median         mean      SE.mean 
##   1.50000000 250.30000000   5.00000000   5.00600000   0.04984957 
## CI.mean.0.95          var      std.dev     coef.var     skewness 
##   0.10017646   0.12424898   0.35248969   0.07041344   0.11297784 
##     skew.2SE     kurtosis     kurt.2SE   normtest.W   normtest.p 
##   0.16782174  -0.45087239  -0.34058520   0.97769855   0.45951315

stat.desc(virginica$Sepal.Length, norm=TRUE)

##      nbr.val     nbr.null       nbr.na          min          max 
##  50.00000000   0.00000000   0.00000000   4.90000000   7.90000000 
##        range          sum       median         mean      SE.mean 
##   3.00000000 329.40000000   6.50000000   6.58800000   0.08992695 
## CI.mean.0.95          var      std.dev     coef.var     skewness 
##   0.18071498   0.40434286   0.63587959   0.09652089   0.11102862 
##     skew.2SE     kurtosis     kurt.2SE   normtest.W   normtest.p 
##   0.16492631  -0.20325972  -0.15354068   0.97117940   0.25831475

stat.desc(immer2$diff, norm=TRUE)

##       nbr.val      nbr.null        nbr.na           min           max 
##   30.00000000    0.00000000    0.00000000  -59.50000000   40.00000000 
##         range           sum        median          mean       SE.mean 
##   99.50000000 -477.40000000  -21.70000000  -15.91333333    4.78742302 
##  CI.mean.0.95           var       std.dev      coef.var      skewness 
##    9.79137947  687.58257471   26.22179579   -1.64778775    0.63503836 
##      skew.2SE      kurtosis      kurt.2SE    normtest.W    normtest.p 
##    0.74379207   -0.51304345   -0.30804332    0.93784528    0.07959101

 

Dans le cas des échantillons appariés, l’asymétrie de la différence est relativement élevée même si elle n’est pas significative. Cela incite à la prudence quant à l’utilisation du test de Student.

 

4.3 Outliers

L’approche classiquement employée mettre en évidence les outliers est celle du boxplot : tout point en dehors de l’intervalle 1.5 X IQR est considéré comme outlier. L’IQR est l’intervalle inter quartile, autrement dit l’étendue entre le 1er et le 3ème quartile. Pour définir les outliers supérieurs, on se place au 3ème quartile (le haut de la boite), on ajoute 1.5 X IQR, et tous les points supérieurs à cette limite sont considérés comme outliers. Pour les outliers inférieurs, on se place au 1er quartile (le bas de la boite) et on soustrait 1,5 IQR. Ces limites correspondent aux traits verticaux des boxplots.

ggplot(setosa, aes(y=Sepal.Length, x=1))+
geom_boxplot()

 

ggplot(virginica, aes(y=Sepal.Length, x=1))+
geom_boxplot()

outlier

 

ggplot(immer2, aes(y=diff, x=1))+
geom_boxplot()

outlier

 

Ici un seul outlier a été mis en évidence, il concerne l’échantillon del’espèce virginica. Pour connaitre la valeur de l’outlier on peut utiliser la fonction boxplot.stat()$out du package grDevices (chargé par défaut).

boxplot.stats(virginica$Sepal.Length)$out
## [1] 4.9

which(virginica$Sepal.Length %in% boxplot.stats(virginica$Sepal.Length)$out)
## [1] 7

La deuxième ligne de commande permet d’identifier l’indice (numéro de ligne) de la donnée outlier.

 

4.4 Egalité des variances

L’égalité des variance ne s’évalue que dans le cas des échantillons indépendants, puisque dans le cas des échantillons appariés on travaille avec la différence. Et on ne l’évalue que si les autres conditions d’application du test de Student ont été remplies (cad normalité, pas d’asymétrie marquée, pas d’outliers significatifs).

 

Le test employé pour comparer les variances est le test F. Son hypothèse nulle est l’égalité des variances. Aussi, pour accepter cette hypothèse la p-value du test doit être > 0.05. Le test F n’est pas très robuste car il est sensible au défauts de normalité. Pour l’employer on utilise
la commande var.test()

var.test(setosa$Sepal.Length, virginica$Sepal.Length)
## 
##  F test to compare two variances
## 
## data:  setosa$Sepal.Length and virginica$Sepal.Length
## F = 0.30729, num df = 49, denom df = 49, p-value = 6.366e-05
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1743776 0.5414962
## sample estimates:
## ratio of variances 
##          0.3072862

 

Ici l’hypothèse d’égalité des variance est rejetée. Le test peut aussi s’employer avec la syntaxe var.test(var~groupe), comme ceci:

iris2 <- iris %>%
    filter(Species %in% c("setosa", "virginica"))

var.test(iris2$Sepal.Length ~ iris2$Species)
## 
##  F test to compare two variances
## 
## data:  iris2$Sepal.Length by iris2$Species
## F = 0.30729, num df = 49, denom df = 49, p-value = 6.366e-05
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1743776 0.5414962
## sample estimates:
## ratio of variances 
##          0.3072862

 

5. Syntaxe des tests

Les test de Student et de Wilcoxon peuvent être employés respectivement avec les fonctions t.test, wilcox.test. Ensuite, se sont les arguments de ces fonctions qui permettent de faire de tests pour échantillons indépendants, échantillons appariés etc…

 

5.1 Les tests de Student

5.1.1 Test de Student pour échantillons indépendants et variances inégales (Test de Welch)

Par défaut la fonction t.test réalise un test pour variances inégales, c’est à dire un test de Welch (c’est écrit sur la première ligne de la sortie).

t.test(iris2$Sepal.Length~iris2$Species)
## 
##  Welch Two Sample t-test
## 
## data:  iris2$Sepal.Length by iris2$Species
## t = -15.386, df = 76.516, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.78676 -1.37724
## sample estimates:
##    mean in group setosa mean in group virginica 
##                   5.006                   6.588

 

On peut aussi utiliser la syntaxe suivante :

t.test(setosa$Sepal.Length, virginica$Sepal.Length)
## 
##  Welch Two Sample t-test
## 
## data:  setosa$Sepal.Length and virginica$Sepal.Length
## t = -15.386, df = 76.516, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.78676 -1.37724
## sample estimates:
## mean of x mean of y 
##     5.006     6.588

 

Dans la sortie logiciel, on retrouve notamment la statistique du test, la p-value, les estimations des deux moyennes, ainsi que l’intervalle de confiance à 95% de la différence des moyennes.

 

5.1.2 Test de Student pour échantillons indépendants et variances égales

Il faut ajouter var.equal=TRUE.

t.test(iris2$Sepal.Length~iris2$Species, var.equal=TRUE)
## 
##  Two Sample t-test
## 
## data:  iris2$Sepal.Length by iris2$Species
## t = -15.386, df = 98, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.786042 -1.377958
## sample estimates:
##    mean in group setosa mean in group virginica 
##                   5.006                   6.588

Remarquez le changement de nom sur la première ligne de la sortie.

 

5.1.3 Test de Student pour échantillon apparié.

Dans ce cas, il faut ajouter paired=TRUE.

t.test(immer2$Y1, immer2$Y2, paired=TRUE)
## 
##  Paired t-test
## 
## data:  immer2$Y1 and immer2$Y2
## t = 3.324, df = 29, p-value = 0.002413
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   6.121954 25.704713
## sample estimates:
## mean of the differences 
##                15.91333

Là encore l’intitulé de la première ligne à été modifié.

 

5.2 Les tests de Wilcoxon

5.2.1 Test de Wilcoxon pour échantillons indépendants

Le test réalisé par défaut est le test pour échantillons indépendants.

wilcox.test(iris2$Sepal.Length~iris2$Species)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  iris2$Sepal.Length by iris2$Species
## W = 38.5, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

 

5.2.2 Test de Wilcoxon pour échantillons appariés.

Comme précédemment, il suffit d’ajouter paired=TRUE.

wilcox.test(immer2$Y1, immer2$Y2, paired=TRUE)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  immer2$Y1 and immer2$Y2
## V = 368.5, p-value = 0.005318
## alternative hypothesis: true location shift is not equal to 0

 

5.3 Choisir une hypothèse bilatérale ou unilatérale

Par défaut les tests utilisent une hypothèse alternative bilatérale. Pour effectuer un test unilatéral, il suffit d’employer l’argument alternative="less"(correspond à m1 < m2 ), ou
alternative="greater"(correspond à m1 > m2).

wilcox.test(setosa$Sepal.Length, virginica$Sepal.Length, alternative="less")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  setosa$Sepal.Length and virginica$Sepal.Length
## W = 38.5, p-value < 2.2e-16
## alternative hypothesis: true location shift is less than 0

wilcox.test(virginica$Sepal.Length, setosa$Sepal.Length, alternative="greater")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  virginica$Sepal.Length and setosa$Sepal.Length
## W = 2461.5, p-value < 2.2e-16
## alternative hypothesis: true location shift is greater than 0

 

6. Pour aller plus loin

Certains, comme Guillaume Rousselet, défendent l’idée que la moyenne n’est un paramètre de tendance centrale adapté à la comparaison des échantillons, car elle est sensible à l’asymétrie et aux outliers. Ils proposent, à la place d’utiliser une moyenne tronquée (trimmed mean), ou encore d’utiliser une approche basée sur la comparaison des percentiles (Wilcox and Rousselet 2017).

 

Et vous, est ce que vous vous utilisez une autre règle pour choisir entre le test de Student et le test de Wilcoxon ? Si oui expliquez la moi dans un commentaire.

 

Dans tous les cas, j’espère que ce tutoriel permettra au plus grand nombre de mieux appréhender la comparaisons de deux moyennes avec le logiciel R.

 
Et si cet article vous a plu ou vous a été utile, partagez le !

Références

Conover, W J. 1999. Practical Nonparametric Statistics. John Wiley.

Wilcox, Rand R, and Guillaume A Rousselet. 2017. “A guide to robust statistical methods in neuroscience.” bioRxiv. Cold Spring Harbor Laboratory. doi:10.1101/151811.

Crédits photos : Pauline Rosenberg

Partager l'article
  •  
  •  
  •  
  •  
  •  
    3
    Partages
  • 3
  •  
  •  

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *