plot of chunk
Lise_Vaudor_headband

A quoi sert-il?

Le test U de Mann-Whitney (aussi appelé test de la somme des rangs de Wilcoxon ou plus simplement test de Wilcoxon) sert à tester l'hypothèse selon laquelle la distribution des données est la même pour deux groupes.

La p-value associée à ce test va ainsi répondre à la question suivante:

Si les données pour les deux groupes étaient issues d'une même population, quelle serait la probabilité que l'on observe par hasard une différence de rangs(ou localisations) entre les deux groupes au moins aussi extrême que celle observée sur les données?

Le test U de Mann-Whitney est souvent utilisé comme solution alternative à l'utilisation d'un test de Student (t-test) dans le cas où les données ne sont pas distribuées selon une loi normale et/ou dans le cas où les données sont peu nombreuses. Il s'agit en effet d'un test non-paramétrique, i.e. un test qui ne repose pas sur une hypothèse de distribution des données.

Un exemple

Considérons les données suivantes:

xA=c(1.5,6.3,2.4,4.1,1.2,5.3,15.2,10.6)
xB=c(2.5,3.3,1.3,2.1,5.7,1.1) 

nA=length(xA)
nB=length(xB)
n=nA+nB
groupes=c(rep("A",nA),rep("B",nB))
x=c(xA,xB)
boxplot(x~groupes, col=c("pink","lightsteelblue3"))

plot of chunk
distribs_simples

On a ici des données qui ne sont ni franchement nombreuses, ni franchement normales.

Comment faire un test U de Mann-Whitney sous R?

montest=wilcox.test(x~groupes)
montest

## 
##  Wilcoxon rank sum test
## 
## data:  x by groupes
## W = 34, p-value = 0.2284
## alternative hypothesis: true location shift is not equal to 0

Ici, la p-value est plus grande que 0.05: on ne parvient pas à rejeter l'hypothèse selon laquelle les données des deux groupes proviennent d'une même distribution au seuil α = 0.05.

Sur quelle métrique le test de Mann-Whitney se base-t-il?

Le test U de Mann-Whitney calcule une p-value qui se base non pas sur les données brutes, mais sur leurs rangs. Les rangs en question correspondent à l'emplacement de chaque donnée au sein de la série lorque cette série est ordonnée par ordre croissant.

rangs=rank(x)
dat=data.frame(groupes,x,rangs)
print(dat)

##    groupes    x rangs
## 1        A  1.5     4
## 2        A  6.3    12
## 3        A  2.4     6
## 4        A  4.1     9
## 5        A  1.2     2
## 6        A  5.3    10
## 7        A 15.2    14
## 8        A 10.6    13
## 9        B  2.5     7
## 10       B  3.3     8
## 11       B  1.3     3
## 12       B  2.1     5
## 13       B  5.7    11
## 14       B  1.1     1

Les distributions des rangs sont beaucoup moins affectées par l'existence de valeurs extrêmes que ne le sont les distributions des valeurs brutes:

boxplot(rangs~groupes, col=c("pink","lightsteelblue3"))

plot of chunk
distribs_rangs

Représentons les rangs des individus ainsi que l'appartenance de ces individus aux deux groupes. Les individus du groupe A sont représentés en rose, les individus du groupe B en bleu. Les pointillés représentent la moyenne des rangs pour chacun des groupes

plot of chunk
rangs

Le test de Mann-Whitney peut se baser sur une métrique W (en lien avec la dénomination "test de Wilcoxon") qui correspond à

    
    W=W_A=\sum_{i=1}^{n}\mathbf{1}_A R_i-n_A(n_A+1)/2\\
    

On peut aussi définir WB :

    
    W_B=\sum_{i=1}^{n}\mathbf{1}_B R_i-n_B(n_B+1)/2
    

Cette notation

    
    \sum_{i=1}^{n}\mathbf{1}_A R_i
     

correspond simplement à la somme des rangs des individus du groupe A...

Les nombres nA(nA + 1)/2 et nB(nB + 1)/2 correspondent en fait à ce que serait la somme des rangs si tous les individus du groupe A (respectivement B) avaient les rangs 1,2,3,....,nA (respectivement 1,2,3,...,nB).

W est ainsi d'autant plus grand que les individus du groupe A présentent des rangs élevés au sein de l'ensemble des individus.

Alternativement, on pourrait considérer la métrique U (en lien avec la dénomination "test U de Mann-Whitney") qui correspond à:

    
    U=min(W_A,W_B)
    

Cette métrique est quant à elle d'autant plus faible que l'un ou l'autre des groupes présente des rangs bas.

Calculons la métrique W pour nos données:

SR_A=sum(rank(x)[which(groupes=="A")])
SR_A0=nA*(nA+1)/2
W_A=SR_A-SR_A0
W=W_A
W

## [1] 34

On retrouve bien la valeur de W qui était indiquée en sortie du test de Mann-Whitney réalisé ci-dessus.

On pourrait aussi calculer WB

SR_B=sum(rank(x)[which(groupes=="B")])
SR_B0=nB*(nB+1)/2
W_B=SR_B-SR_B0

En l'occurence, par construction on a toujours

    W_A+W_B=n_An_B

On peut le vérifier ici:

W+W_B

## [1] 48

nA*nB

## [1] 48

A quoi la p-value correspond-elle?

Supposons que le groupe A ne corresponde pas à des rangs particulièrement élevés dans la série. Dans ce cas, les rangs observés pour A sont comme "attribués au hasard". Quelle est alors la distribution de W?

Pour répondre à cette question nous allons simuler un grand nombre d'échantillons au hasard

  • correspondant à des rangs entre 1 et nA + nB = 14
  • tels que nA = 8 des individus soient étiquetés comme appartenant au groupe A
  • tels que nB = 6 des individus soient étiquetés comme appartenant au groupe B

Nous calculerons pour chaque simulation la valeur de W.

r=1:n  # les rangs prennent les valeurs 1,2,3,...,n

niter=10000
vecteur_W=rep(NA,niter)
for (i in 1:niter){
  # On réattribue au hasard les rangs aux données (i.e. aux différents groupes)
  r_sim=sample(r,n,replace=F)  
  # On calcule la valeur de W pour ces rangs simulés
  SR_A=sum(r_sim[which(groupes=="A")])
  SR_A0=nA*(nA+1)/2
  Wsim=SR_A-SR_A0
  # On stocke la valeur de W pour chaque itération itération dans vecteur_W
  vecteur_W[i]=Wsim
}

hist(vecteur_W,freq=F,
     col="lemonchiffon",
     main="distrib. de W sous hypothèse nulle")
abline(v=W, lty=3,lwd=3, col="red")

plot of chunk
hist_W

La distribution empirique de W nous permet de calculer la p-value correspondant au test de l'hypothèse "A n'a pas une localisation particulièrement élevée".

wilcox.test(x~groupes,alternative="less")

## 
##  Wilcoxon rank sum test
## 
## data:  x by groupes
## W = 34, p-value = 0.9094
## alternative hypothesis: true location shift is less than 0

print(length(which(vecteur_W>=montest$statistic))/niter)

## [1] 0.1145

Ici on veut en fait faire un test bilatéral : notre hypothèse nulle "la distribution des données est la même pour les deux groupes" correspond en fait à deux sous-hypothèses:

  • "A n'a pas une localisation particulièrement élevée" ET
  • "A n'a pas une localisation particulièrement basse" (i.e. "B n'a pas une localisation particulièrement élevée")

On s'interroge ainsi non seulement sur la distribution de WA mais aussi sur la distribution de WB.

Or ces deux distributions sont en fait les mêmes... Vous pouvez le vérifier comme suit:

layout(matrix(1:2,nrow=1))
vecteur_W_A=vecteur_W
vecteur_W_B=nA*nB-vecteur_W_A

hist(vecteur_W_A,freq=F,
     col="lemonchiffon",
     main="distrib. de W_A")
hist(vecteur_W_B,freq=F,
     col="lemonchiffon",
     main="distrib. de W_B")

plot of chunk
hist_W_A_W_B

On peut ainsi retrouver la p-value fournie par le test:

print(2*length(which(vecteur_W>=montest$statistic))/niter)

## [1] 0.229

print(montest$p.value)

## [1] 0.2284

Comment R calcule-t-il la p-value?

Ci-dessus, j'ai montré comment on peut recalculer "à la main" la p-value fournie par le test. En fait, pour des valeurs de nA et nB assez élevées, la distribution de W est approximativement normale (R n'a pas besoin d'effectuer des simulations, comme nous l'avons fait, pour calculer la p-value).

Pour des valeurs de nA et nB relativement faibles, les p-values sont tabulées (c'est-à-dire que pour des valeurs nA et nB faibles, R a en mémoire les quantiles de la distribution de W). Cependant, de la même manière que dans nos simulations ci-dessus, ces p-values ont été calculées en supposant qu'il n'y avait pas d'ex-aequos (i.e. pour calculer ces quantiles, il a été supposé que les rangs étaient 1,2,3,....,n, et non, par exemple, 2,2,2,3,4,5,...,n). C'est la raison pour laquelle R affiche des Warnings dans le cas où

  • les données sont peu nombreuses ET
  • il y a des ex-aequos.

xA=c(1.5,6.3,6.3,2.7)
xB=c(2.5,3.3,1.3,2.1,5.7,1.1) 
nA=length(xA)
nB=length(xB)
n=nA+nB
groupes=c(rep("A",nA),rep("B",nB))
x=c(xA,xB)

wilcox.test(x~groupes)

## Warning: cannot compute exact p-value with ties

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  x by groupes
## W = 18, p-value = 0.2395
## alternative hypothesis: true location shift is not equal to 0

Pour des échantillons plus grands, on n'a pas ce même message d'erreur car la p-value est calculée en se basant sur l'approximation normale:

xA=rlnorm(50,3,2)
xB=rlnorm(50,4,2)
xA[1:10]=xA[11:20]
nA=length(xA)
nB=length(xB)
n=nA+nB
groupes=c(rep("A",nA),rep("B",nB))
x=c(xA,xB)

wilcox.test(x~groupes)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  x by groupes
## W = 1244, p-value = 0.9698
## alternative hypothesis: true location shift is not equal to 0