ANOVA non paramétrique à 2 facteurs ou le test de Scheirer-Ray-Hare

anova non paramétrique à 2 facteurs

Jusque très récemment, je ne savais pas qu’il existait une approche permettant de réaliser une ANOVA non paramétrique à deux facteurs. C’est un lecteur du blog (Philippe V.) qui me l’a fait découvrir. Cette alternative, c’est le test de Scheirer-Ray-Hare ; je n’en avais jamais entendu parler.

J’ai fait quelques recherches pour en savoir plus, et je n’ai pas trouvé grand-chose. Je suis quand même tombée sur une page wikipédia. Elle est plutôt succincte, mais il y a quelques informations intéressantes :

  • le test est basé sur les rangs (classique pour une approche non paramétrique) ; par contre la statistique du test n’est pas décrite

  • la puissance de ce test (capacité à mettre en évidence un effet significatif lorsque cet effet existe réellement) serait bien inférieure à celle l’ANOVA paramétrique à deux facteurs. Cela serait la raison, pour laquelle le test Scheirer-Ray-Hare n’est pas populaire (en plus du fait qu’il a été décrit après les approches paramétriques et non paramétriques de l’analyse de variance (à un facteur j’imagine – c’est-à-dire le test de Kruskal Wallis).

Philippe m’avait également fait suivre cette page du site rcompanion dédiée à l’utilisation du test de Scheirer-Ray-Hare avec R, ce qui m’a permis de le mettre en oeuvre facilement.

Dans cet article, je vous propose donc d’explorer un peu ce test, en comparant les résultats obtenus avec ceux de l’ANOVA paramétrique à deux facteurs et ceux obtenus en employant une approche par permutation, le tout sur deux jeux de données que je connais bien :

  • warpbreaks : contenu dans le package dataset, et qui concerne les ruptures de fils sur un métier à tisser en fonction de deux types de laine et de trois niveaux de tension.

  • ToothGrowth : contenu dans le package dataset, et qui concerne la longueur des cellules de croissance des dents chez des cobayes qui reçoivent de la vitamine C selon trois doses et deux méthodes d’administration (jus d’orange ou comprimé) différentes.

Comme je connais ces deux jeux de données, je sais que :

  • Les interactions entre les deux facteurs étudiés sont significatives (lorsqu’on emploie l’anova paramétrique). L’évaluation de l’interaction est importante, car les analyses subséquentes peuvent être différentes si l’effet est statistiquement significatif, ou non.

  • Et que les conditions d’applications de l’anova paramétrique sont satisfaites pour ToothGrowth, mais pas pour warpbreaks (l’hypothèse d’homogénéité est rejetée) ; le test de Scheirer-Ray-Hare aura donc tout son sens.

Table des matières

Analyse des données ToothGrowth

Data

# renommer
TG <- ToothGrowth
TG$dose <- as.factor(TG$dose)

library(tidyverse)
TG %>% 
    group_by(supp, dose) %>% 
    count()
## # A tibble: 6 x 3
## # Groups:   supp, dose [6]
##   supp  dose      n
##   <fct> <fct> <int>
## 1 OJ    0.5      10
## 2 OJ    1        10
## 3 OJ    2        10
## 4 VC    0.5      10
## 5 VC    1        10
## 6 VC    2        10 

Le jeu de données est équilibré avec 10 données par conditions (1 dose, 1 mode d’administration). Nous pouvons facilement visualiser les données comme ceci :

library(RcmdrMisc)
plotMeans(TG$len,TG$dose,TG$supp) 
jeux de données ToothGrowth

Le graphique nous montre que l’évolution de la réponse en fonction de la dose n’est pas identique pour les deux modes d’administration, ce qui laisse supposer une interaction. L’évolution se faisant dans la même direction, on parlera plutôt d’interaction quantitative.

Anova paramétrique à 2 facteurs

TG.aov <- aov(len~dose*supp, 
              contrasts = list(dose=contr.sum,
                             supp=contr.sum),
              data=TG)
Anova(TG.aov, type=3)
## Anova Table (Type III tests)
## 
## Response: len
##              Sum Sq Df  F value    Pr(>F)    
## (Intercept) 21236.5  1 1610.393 < 2.2e-16 ***
## dose         2426.4  2   92.000 < 2.2e-16 ***
## supp          205.4  1   15.572 0.0002312 ***
## dose:supp     108.3  2    4.107 0.0218603 *  
## Residuals     712.1 54                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Les tests de tous les effets sont significatifs, notamment celui de l’interaction (pval =0.022).

Nous pouvons évaluer les conditions d’application de l’ANOVA paramétrique (normalité et homogénéité des résidus) de manière visuelle et en employant les tests de Shapiro-Wilk (Normalité des résidus) et Bartlett (homogénéité des résidus) :

library(performance)
check_model(TG.aov) 
Evaluation des hypothèse de normalité et d'homogénéité
# test de shapiro wilk sur les résidus 
check_normality(TG.aov)
## OK: residuals appear as normally distributed (p = 0.669).

# test de Bartlet :
check_homogeneity(TG.aov)
## OK: There is not clear evidence for different variances across groups (Bartlett Test, p = 0.226). 

Les deux hypothèses sont satisfaites.

Anova non paramétrique à 2 facteurs : Test de Scheirer-Ray-Hare.

Le test de Scheirer-Ray-Hare peut être mis en place avec la fonctioncheirerRayeHare() du package rcompanion.

### Scheirer–Ray–Hare test

library(rcompanion)
scheirerRayHare(len ~ dose * supp,
                data = TG, type=2)
                
                
## 
## DV:  len 
## Observations:  60 
## D:  0.999222 
## MS total:  305
##           Df  Sum Sq      H p.value
## dose       2 12394.4 40.669 0.00000
## supp       1  1050.0  3.445 0.06343
## dose:supp  2   515.5  1.692 0.42923
## Residuals 54  4021.1                 

Contrairement à l’anova paramétrique à deux facteurs, ici, seul l’effet de la dose est significatif au seuil de 5%. Et on peut voir que la pvalue de l’interaction est très élevée (pvalue = 0.43).

Anova à 2 facteurs par permutations

library(lmPerm)
set.seed(1234)
TG.aovp <- aovp(len~dose*supp, 
              contrasts = list(dose=contr.sum,
                             supp=contr.sum),
              data=TG)
## [1] "Settings:  unique SS "
summary(TG.aovp)
## Component 1 :
##             Df R Sum Sq R Mean Sq Iter Pr(Prob)    
## dose1        2  2426.43   1213.22 5000  < 2e-16 ***
## supp1        1   205.35    205.35 5000  < 2e-16 ***
## dose1:supp1  2   108.32     54.16 3123  0.04195 *  
## Residuals   54   712.11     13.19                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Ici, comme avec l’ANOVA paramétrique à deux facteurs, tous les effets sont significatifs au seuil de 5%, et notamment celui de l’interaction (pvalue =0.042).

Analyse des données warpbreaks

Data

View(warpbreaks)
WB <- warpbreaks
library(tidyverse)
WB %>% 
    group_by(wool, tension) %>% 
    count()
## # A tibble: 6 x 3
## # Groups:   wool, tension [6]
##   wool  tension     n
##   <fct> <fct>   <int>
## 1 A     L           9
## 2 A     M           9
## 3 A     H           9
## 4 B     L           9
## 5 B     M           9
## 6 B     H           9 

Le jeu de données est équilibré avec 9 données par conditions (1 type de laine, 1 tension). Nous pouvons facilement visualier les doonnées comme ceci :

library(RcmdrMisc)
plotMeans(WB$breaks,WB$tension,WB$wool) 
Repérsentation des données warpbreaks

On peut voir que l’évolution du nombre de ruptures des fils, en fonction du niveau de tension, n’est pas identique pour les deux types de laines, puisque les profils se croisent. On peut donc suspecter la présence d’une interaction de type qualitative entre le type de laine et la tension.

Anova paramétrique à deux facteurs

WB.aov <- aov(breaks~wool*tension, 
              contrasts = list(wool=contr.sum,
                             tension=contr.sum),
              data=WB)
Anova(WB.aov, type=3)
## Anova Table (Type III tests)
## 
## Response: breaks
##              Sum Sq Df  F value    Pr(>F)    
## (Intercept)   42785  1 357.4672 < 2.2e-16 ***
## wool            451  1   3.7653 0.0582130 .  
## tension        2034  2   8.4980 0.0006926 ***
## wool:tension   1003  2   4.1891 0.0210442 *  
## Residuals      5745 48                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Ici, l’effet de l’interaction et de la tension sont significatifs au seuil de 5% et celui de la laine au seuil de 10% (pvalue = 0.058). Nous pouvons vérifier les hypothèses de normalité et d’homogénéité des résidus.

library(performance)
check_model(WB.aov) 
Diagnostic de régression avant l'anova non paramétrique à deux facteurs
check_normality(WB.aov)
## OK: residuals appear as normally distributed (p = 0.816).
check_homogeneity(WB.aov)
## Warning: Variances differ between groups (Bartlett Test, p = 0.024). 

L’hypothèse de normalité est acceptée, mais pas l’hypothèse d’homogénéité. Cela n’est pas étonnant, car la variable réponse est de type comptage. Ces données suivent théoriquement une distribution de Poisson qui spécifie que la variance est proportionnelle à la moyenne. Ainsi, plus le nombre de ruptures est important, plus la variabilité est importante. C’est ce qu’on voit sur le graph “Homogeneity of Variance”).

Anova non paramétrique à 2 facteurs : Test de Scheirer-Ray-Hare.

### Scheirer–Ray–Hare test

library(rcompanion)

scheirerRayHare(breaks~wool*tension,
                data = WB, type=2)
## 
## DV:  breaks 
## Observations:  54 
## D:  0.9980941 
## MS total:  247.5
##              Df Sum Sq       H  p.value
## wool          1  327.6  1.3261 0.249508
## tension       2 2670.2 10.8093 0.004496
## wool:tension  2  899.8  3.6427 0.161810
## Residuals    48 9194.9 

Avec le test de Scheirer-Ray-Hare, seul l’effet de la tension est significatif au seuil de 5% ; l’effet l’interaction n’est pas mis en évidence, même au seuil de 10% (pvalue=0.16).

 

Anova à deux facteurs par permutation

library(lmPerm)
set.seed(1234)
WB.aovp <- aovp(breaks~wool*tension, 
              contrasts = list(wool=contr.sum,
                             tension=contr.sum),
              data=WB)
## [1] "Settings:  unique SS "
summary(WB.aovp)
## Component 1 :
##                Df R Sum Sq R Mean Sq Iter Pr(Prob)    
## wool1           1    450.7    450.67 2452  0.03956 *  
## tension1        2   2034.3   1017.13 5000  0.00020 ***
## wool1:tension1  2   1002.8    501.39 5000  0.01660 *  
## Residuals      48   5745.1    119.69                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Les résultats sont très proches de ceux obtenus avec l’ANOVA paramétrique à deux facteurs : les effets de l’interaction et de la tension sont significatifs au seuil de 5%. Il en est de même pour la laine dont la pvalue est de 4%.

Remarque

Certains auteurs rapportent que les approches par permutations nécessitent l’homogénéité des résidus (je pense que c’est vrai, car sans cela, les étiquettes ne sont pas interchangeables puisque la probabilité que certaines valeurs appartiennent à certains groupes n’est pas équiprobable). Cette hypothèse n’est pas satisfaite ici ;  les résultats manquent alors sans doute de robustesse. Néanmoins, à titre exploratoire, nous pouvons quand même constater que les résultats obtenus sont, là encore, très similaires à ceux de l’approche paramétrique, et en ce qui concerne l’interaction, contraires à ceux du test de Scheirer-Ray-Hare.

Conclusion

Les résultats obtenus avec les deux jeux de données semblent bien confirmer le manque de puissance du test de Scheirer-Ray-Hare, notamment pour détecter l’interaction entre les deux facteurs.

Cette lacune est majeure, car l’évaluation de la significativité de l’interaction est l’élément clé de l’ANOVA à deux facteurs ; la suite de l’analyse en dépend.

Au final, je ne pense pas me resservir beaucoup de ce test. En cas de défaut important de la normalité, il me semble préférable d’employer une approche par permutation. Et en cas de défaut d’homogénéité des variances, je pense que j’essaierai plutôt d’employer une transformation (log ou boxcox). Mais si celle-ci ne fonctionne pas, le test de Scheirer-Ray-Hare, pourrait quand même être utile (mais je pense que je vérifierais la robustesse des résultats en employant également une anova paramétrique et une approche par permutation ! )

Dites-moi en commentaire ce que vous pensez de ce test de Scheirer-Ray-Hare.

Poursuivez votre lecture

Soutenir le blog

Vous souhaitez soutenir mon travail ? Vous pouvez le faire en réalisant un don libre sur la page Tipeee du blog.

11 réponses

  1. Génial, je vais explorer ça !

    Cela signifie que l’on peut explorer les effets d’interactions, ce qui est compliqué avec les analyses non-paramétriques.

    En revanche, pourriez-vous indiquer comment on rédige les résultats du test dans un compte-rendu ou article scientifique.

    Il faut indiquer la valeur du H et la p-value j’imagine ?

    1. Bonjour Laurence,

      oui on peut explorer les interactions, mais en gardant à l’esprit que le test manque de puissance. C’est pour ça que je testerai la robustesse des résultats avec l’anova paramétrique et le test de permutation.
      Concernant la rédaction des résultats, je pense que l’on peut s’inspirer de ce que fait le package report. Vous trouverez un exemple dans cet article “Obtenir un reporting automatique des analyses statistiques

  2. Bonjour Claire et merci pour ces explications. Finalement dans cet exemple, vous préférez opter pour une ANOVA paramétrique malgré une absence d’homogénéité des variances. Dans ce cas là, n’est-il pas tout simplement préférable de partir sur un test de Kruskal-Wallis ?

    Bien cordialement

    1. Bonjour,

      Le test de Kruskall-Wallis est une anova non paramétrique à un seul facteur. Ici on est sur une anova non paramétrique à 2 facteurs. Et comme le test manque de puissance, je pense que c’est pas mal de vérifier la robustesse des résultats avec une méthode paramétrique, même si les hypothèses ne sont pas satisfaites.
      Bonne continuation

  3. Merci Claire pour ce partage et les visualisations avec le package performance.
    Rcompanion est une telle mine qu’il y a toujours à découvrir… tout comme DellaData.fr !

  4. Merci beaucoup pour cette page comme toujours extrêmement didactique et source d’inspirations.

    Vous avez pour clarifier l’aspect “nouveauté” dans la troisième édition du Sokal & Rohl (1991) en page 445-447 un exemple traité de ce test
    avec les calculs à la main (box 13.12) pour bien en comprendre construction et logique.
    L’ouvrage de Jerrold Zar (1984) présente également ce test en chapitre 13.4, p. 219 avec une approche de comparaisons multiples remaniée dans les éditions ultérieures.
    En tout cas dans le domaine de la biologie expérimentale, voilà un test bien intéressant.

  5. Claire, je vous tire mon chapeau. Vous aviez une méthodo simple pour faire passer la pilule.
    J’ai beaucoup aimé ce tuto.
    Bon vent!!

  6. Je pense avoir oublié de citer mes sources :
    Sokal, R.R. & Rohlf, F.J., 1995. Biometry. The principles and practice of statistics in biological research. Third edition. New York, Freeman, W.H. and Co. : 887 p.
    qui fait duite à l’édition de 198 en deux volumes (livre + tables)
    et
    Zar, J.H., 1984. Biostatistical analysis. Englewood Cliffs, New Jersey, Prentice Hall : 718 p.
    avec sa nouvelle édition
    Zar, J.H., 2014. Biostatistical analysis. Harlow, UK, Pearson Education Limited : 756 p.
    Cordialement

    1. Bonjour Pierre Guy,
      merci pour ces compléments d’information.
      J’ajoute aussi la publication originale :
      Scheirer, C. J., Ray, W. S., & Hare, N. (1976). The analysis of ranked data derived from completely randomized factorial designs. Biometrics, 429-434.

Laisser un commentaire

Votre adresse e-mail ne sera pas publiée.

Bonjour !

vous venez souvent ?

Identifiez-vous pour avoir accès à toutes les fontionnalités !

Aide mémoire off'R ;)

Enregistrez vous pour recevoir gratuitement mes fiches “aide mémoire” (ou cheat sheets) qui vous permettront de réaliser facilement les principales analyses biostatistiques avec le logiciel R et pour être informés des mises à jour du site.