Tests de permutation avec le logiciel R
Note importante aux lecteurs : Depuis la première publication de cet article, deux lecteurs (Corto et Fanfoué) ont apporté des commentaires importants concernant les tests de permutation. Je vous recommande vivement de prendre le temps de lire et de prendre en compte leurs remarques (disponibles à la fin de cet article ou dans la section des commentaires) avant de mettre en œuvre les méthodes décrites ici. Leurs mises en garde apportent une dimension critique à la compréhension et à l’application des tests de permutation.
Les tests de permutation sont des approches robustes, basées sur l’aléatoire et le ré-échantilonnage. Elles permettent de réaliser des tests “classiques”, notamment des procédures paramétriques, sans que la validité des résultats ne reposent sur des distributions théoriques !
Concrètement, ça veut dire par exemple, qu’il est possible de réaliser une régression linéaire multiple, et d’avoir confiance dans ses résultats, même si les résidus ne suivent pas une loi normale, qu’ils ne sont pas distribués de façon homogène, ou encore que des outliers sont présents dans les données.
Table des matières
Rappels sur les tests d'hypothèses classiques
Prenons l’exemple du test t de Student.
Hypothèses et statistique du test
Les hypothèses nulle et alternative du test t de Student sont:
\[H_0 : m_A = m_B\]
\[H_1 : m_A \neq m_B\]
La statistique du test est :
\[T=\frac{m_A – m_B }{s_p\sqrt{1/n_A + 1/n_B}} \sim t_{n_A+n_B-2}\]
Avec :
\[s_p =\sqrt \frac{(n_A-1)s_A^2 + (n_B-1)s_B^2}{n_A+n_B-2}\]
Déroulement du test
Comme expliqué dans un de mes précédents articles les tests d’hypothèses classiques se déroulent en 3 étapes:
- Validation des hypothèses nécessaires à l’utilisation du test considéré. Pour le test t de Student, il s’agit de vérifier que les données des deux groupes suivent une loi normale, et que les variances sont homogènes. Ces hypothèses sont nécessaires pour que la statistique T, calculée à l’étape 2, suive la distribution de Student (c’est la théorie qui dit ça).
- Calcule de la statistique du test, à partir des données observées.
- Comparaison de la valeur observée à la distribution théorique de cette statistique, sous l’hypothèse nulle qui stipule que les moyennes ne sont pas différentes. Il s’agit d’une distribution de Student centrée sur 0, et dont le nombre de degrés de liberté dépend du nombre d’observations dans les groupes A et B.
- Si la valeur observée de la statistique se situe à droite du percentile 97.5 (ou à gauche du percentile 2.5) de cette distribution théorique centrée sur 0, alors le test est significatif, et on concluera que les deux moyennes sont différentes.
- sinon, le test n’est pas significatif, et on concluera que les deux moyennes ne sont pas différentes.
Les tests de permutation
Avec un test de permutation on ne va plus se baser sur la distribution théorique de la statistique pour évaluer la significativité, mais sur une distribution empirique. Elle est dite empirique parce qu’elle est construite directement partir des données.
Déroulement
Imaginons une expérimentation dans laquelle 10 souris sont soumises à un régime alimentaire A, et 10 autres à un régime alimentaire B. Et que l’on souhaite comparer les moyennes des masses des souris des deux groupes, en employant un test t de Student par exemple. Les données pourraient être celles-ci :
## 1 souris_ 1 A 65.28
## 2 souris_ 2 A 61.95
## 3 souris_ 3 A 85.05
## 4 souris_ 4 A 92.30
## 5 souris_ 5 A 107.09
## 6 souris_ 6 A 112.05
## 7 souris_ 7 A 88.86
## 8 souris_ 8 A 84.23
## 9 souris_ 9 A 72.95
## 10 souris_ 10 A 110.52
## 11 souris_ 11 B 85.19
## 12 souris_ 12 B 82.31
## 13 souris_ 13 B 66.59
## 14 souris_ 14 B 73.48
## 15 souris_ 15 B 71.03
## 16 souris_ 16 B 73.46
## 17 souris_ 17 B 71.76
## 18 souris_ 18 B 79.99
## 19 souris_ 19 B 53.46
## 20 souris_ 20 B 55.57 
Les étapes
Le test de permutation va consister à :
- Calculer les moyennes des masses observées, m_A et m_B, puis à calculer la statistique T du test de Student (comme décrit plus haut). Cette première statistique calculée sur les valeurs observées est généralement appelé t0 (c’est seulement une convention). Ici t0 = 1.59
- Placer virtuellement les masses des 20 souris dans un même sac, puis à attribuer aléatoirement 10 masses au régime alimentaire A, et les 10 masses restantes au régime alimentaire B. Par exemple, les masses des souris 13, 6, 12, 19, 9, 11, 7, 4, et 10, pourraient être attribués au régime A, et les autres au régime B.
- Calculer les nouvelles m_A e et m_B puis la nouvelle statistique T.
- Recommencer les étapes 2 et 3 pour toutes les différentes façons possibles de constituer les groupes. Ce nombre correspond au nombre de combinaisons mathématiques. Il s’obtient avec la fonction `choose(20,10)`. Dans cet exemple il y a 184756 façons possibles de constituer les groupes A et B, autrement dit à la fin de l’étape 4, 184756 statistiques T auront été estimées. On les nomme généralement t.
- L’ensemble des statistiques t ainsi estimées sont alors ordonnées ; elles forment la distribution empirique.
Ces tests sont appelés test de permutations car les différentes combinaisons peuvent aussi être vues comme des permutations des étiquettes des régimes (A ou B),entre les lignes
Calcul des p-values
La p-value bilatérale du test est estimée en calculant la proportion de ces statistiques qui se situent au dela de la valeur absolue de statistique t0 caculées sur les données observées. Elle correspond à l’aire sous la courbe au dela des valeurs t0 et -t0
Ici la p-value bilatérale est égale à 0.02, on conclut donc que les deux moyennes sont différentes.
Pour un test unilatéral, il suffit de considérer un seul coté, c’est à dire la proportion de statistique t au dela de t0 OU au delà de -t0. La p-value unilatérale à droite est ici égale à 0.01.
Principe
Les permutations ont pour conséquence de casser les relations entre la variable dépendante (ici la masse des souris) et la ou les variables indépendantes (ici le groupe). Lorsque cette relation est détruite, on est dans la situation de l’hypothèse nulle, qui postule que le groupe n’a pas d’effet. Les tests de permutation permettent donc de construire la distribution empirique de la statistique sous l’hypothèse nulle.
Remarques
Lorsque la distribution empirique sous l’hypothèse nulle est basée sur toutes les permutations possibles, le test de permutation est dit “exact”. En effet, lorsque l’hypothèse nulle est rejetée, elle l’est à un niveau alpha exact.
Le nombre de permutations possibles est généralement très élevé. Dans ce cas là, seul un échantillon des permutations possibles est constitué, en employant des approches de type Monte Carlo. On parle alors de test de permutation approximatif.
Les permutations reposent sur un générateur de nombres pseudo aléatoires. Ce nombre change chaque fois que le test est lancé. Si un test de permutation est réalisé plusieurs fois sur les mêmes données, les résultats obtenus seront légèrement différents. Pour éviter ce problème et obtenir des résultats stritement identiques,il est nécessaire de fixer ce nombre à l’aide de la fonction set.seed. Par exemple set.seed(1234) ou set.seed(678). N’importe quel nombre peut être utilisé.
Mise en application sous R
Il existe plusieurs packages sous R pour réaliser des tests de permutations. Les plus utilisés sont coin, lmPerm, et pgirmess.
Le package coin permet surtout de réaliser des tests statistiques classiques comme le test t, le test de Wilcoxon ou encore le test de Kruskal Wallis. Le package lmPerm est davantage dédié aux modèles linéaires. Le package pgirmess, quant à lui, permet d’ étendre les tests de permutations aux modèles linéaires, aux modèles linéaires généralisés, ainsi qu’aux modèles linéaires mixtes.
Le package coin
Liste des fonction disponibles du package coin
Exemple d'utilisation du package coin
Nous allons, ici, réaliser l’equivalent du test t de Student par un test de permutation. Cela nécessite la fonction oneway_test.
library(coin)
oneway_test(Masse ~ Regime, data=mydf, distribution="exact")
##
## Exact Two-Sample Fisher-Pitman Permutation Test
##
## data: Masse by Regime (A, B)
## Z = 2.2412, p-value = 0.02125
## alternative hypothesis: true mu is not equal to 0 Ces résultats peuvent être comparer avec ceux du test t :
t.test(Masse ~ Regime, data=mydf, var.equal=TRUE)
##
## Two Sample t-test
##
## data: Masse by Regime
## t = 2.5434, df = 18, p-value = 0.02038
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 2.912846 30.575154
## sample estimates:
## mean in group A mean in group B
## 88.028 71.284 Ici, les résultats sont plutôt proches. Néanmoins quand on regarde les distributions des masses, on s’aperçoit que 2 outliers sont présents dans le groupe du régime B. Cela pourrait biaiser les résultats. Le test de permutation est donc un moyen de les valider.
De façon générale, pour savoir comment utiliser les fonctions d’un package, vous pouvez consulter mon article 7 façons d’obtenir de l’aide avec le logiciel R.
Le package lmPerm
Le package lmPerm est dédié aux modèles linéaires. Il permet, en autres, de réaliser des ANOVA, des ANCOVA, des régressions linéaires simples ou encore de la régression linéaires multiples, sous la forme de tests de permutation.
Les fonction du package lmPerm
Ce package comportent deux fonctions principales :
lmp: pour réaliser des régressions linéaires simples ou multiplesaovp: pour réaliser des ANOVA à un ou deux facteurs, ou encore
des ANCOVA
Exemple d'utilisation du package lmPerm
ANOVA à un facteur (one-way)
L’ANOVA nécessite que les données de chaque groupe suivent une loi normale et que les variances soient homogènes. Ce n’est clairement pas le cas ici. Dans cette situation, il peut être intéressant d’utiliser un test de permutation.
library(lmPerm)
mod <- aovp(response~group, perm="Prob", data=mydf2)
## [1] "Settings: unique SS "
summary(mod)
## Component 1 :
## Df R Sum Sq R Mean Sq Iter Pr(Prob)
## group 2 29713 14856.5 5000 0.0348 *
## Residuals 27 136109 5041.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 A comparer avec les résultats du test classique :
mod2 <- aov(response~group, data=mydf2)
summary(mod2)
## Df Sum Sq Mean Sq F value Pr(>F)
## group 2 29713 14857 2.947 0.0696 .
## Residuals 27 136109 5041
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Il est également possible de réaliser une ANCOVA, ou encore un ANOVA à deux facteurs avec interaction, avec la même fonction aovp.
Régression linéaire multiple
set.seed(9876)
fit <-lmp(Sepal.Length ~ Petal.Length +Petal.Width , data=iris)
## [1] "Settings: unique SS : numeric variables centered"
summary(fit)
##
## Call:
## lmp(formula = Sepal.Length ~ Petal.Length + Petal.Width, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.18534 -0.29838 -0.02763 0.28925 1.02320
##
## Coefficients:
## Estimate Iter Pr(Prob)
## Petal.Length 0.5418 5000 <2e-16 ***
## Petal.Width -0.3196 1817 0.0523 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4031 on 147 degrees of freedom
## Multiple R-Squared: 0.7663, Adjusted R-squared: 0.7631
## F-statistic: 241 on 2 and 147 DF, p-value: < 2.2e-16 D’autres exemples de l’utilisation de la fonction lmp sont présents dans le chapitre 12 du livre
“R In action” de Robert I. Kabacoff.
Le package pgirmess
Le package pgirmess propose la fonction PermTestqui permet de réaliser des tests de permutation après l’ajustement de modèles linéaires, des modèles linéaires généralisés, et de modèles linéaires mixtes (à partir d’un objet lme construit avec le package nlme).
Régression linéaire multiple
library(pgirmess)
set.seed(9876)
my_lm <-lm(Sepal.Length ~ Petal.Length +Petal.Width , data=iris)
PermTest(my_lm, B=10000)
##
## Monte-Carlo test
##
## Call:
## PermTest.lm(obj = my_lm, B = 10000)
##
## Based on 10000 replicates
## Simulated p-value:
## p.value
## Petal.Length 0.0000
## Petal.Width 0.0482 L’argument B permet de définir le nombre de permutations réalisées.
Modèle linéaire généralisé
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
glm.D93 <- glm(counts ~ outcome + treatment, family=poisson)
PermTest(glm.D93,B=100)
##
## Monte-Carlo test
##
## Call:
## PermTest.glm(obj = glm.D93, B = 100)
##
## Based on 100 replicates
## Simulated p-value:
## p.value
## outcome 0.12
## treatment 1.00 Modèle linéaire à effet mixte
library(nlme)
my_lme <- lme(score ~ Machine , random=~1| Worker, data = Machines)
PermTest(my_lme,B=5000)
##
## Monte-Carlo test
##
## Call:
## PermTest.lme(obj = my_lme, B = 5000)
##
## Based on 5000 replicates
## Simulated p-value:
## p.value
## (Intercept) 0.9994
## Machine 0.0000 Pour aller plus loin
Si vous voulez en savoir davantage sur les tests de permutation et les procédures de rééchantillonage, je vous conseille le livre de référence “An Introduction to the Bootstrap”, d’Efron et Tibshirani.
Pour une vision plus synthétique, je vous conseille le chapitre 12 du livre “R In action” de Robert I. Kabacoff.
Conclusion
De façon générale, les tests de permutation peuvent être employés chaque fois que :
- les échantillons sont de petites tailles,
- les données ne suivent pas bien les distributions théoriques supposées,
- des outliers sont présents.
De mon point de vue, les tests de permutation sont donc surtout utiles lorsqu’il n’y a pas d’alternative non paramétrique. C’est le cas par exemple de l’ANOVA à deux facteurs, de la régression linéaire multiple, des modèles linéaires généralisés ou encore des modèles linéaires à effets mixtes.
Il me semble que les tests de permutation peuvent également être très utiles dans le cadre de procédures répétées. Par exemple lorsqu’on doit réaliser exactement les mêmes analyses mais pour différentes variables réponses. Vérifier les hypothèses de normalité et / ou d’homogénéité peut vite devenir fastidieux. Du coup, automatiser les analyses en utilisant ce type de tests peut être une bonne solution. Je vais d’ailleurs l’essayer très prochainement.
Et vous, utilisez-vous régulièrement des tests de permutation ? Si oui dans quel contexte ?
Si cet article vous a plu, ou vous a été utile, et si vous le souhaitez, vous pouvez soutenir ce blog en faisant un don sur sa page Tipeee 
Note : Je touche une petite commission (entre 3 et 6%) si vous passez par les liens Amazon de cet article pour acheter le livres mentionné. cela m’aide à entretenir ce blog, merci si vous le faites ! 
Crédits photos : aitoff
Bonjour,
j’ai une question concernant le test par permutation. Est-ce un problème si les groupes ne sont pas de même taille? Par exemple, pour la fonction choose() , je ne vois pas bien comment faire pour déterminer le nombre de permutations possibles (pour mes data : control = 150 et test =76).
De plus, et dans le même cadre, auriez vous la possibilité de faire un article sur le nouveau package “dabestr” qui est vraiment bien pour la visualisation des données avec l’intervalle de confience calculé avec bootstrap?
Merci pour tout ce que vous faites, c’est vraiment une mine d’or pour un étudiant en thèse comme moi !
Bonne continuation,
Yoann
Merci pour les détails, je me retrouve mieux avec vos différentes explications
Bonjour,
Merci pour vos articles. Je voulais vous demander est-ce que le même cas peut s’appliquer dans une manova
Bonjour,
Je viens de lire votre passage sur les Modèle linéaire à effet mixte et j’ai du mal à comprendre comment votre p-value peut être égale à 0. Pourriez vous donner des détails sur ce chiffre?
Cordialement
Bonsoir
la pvalue est arrondie à 0 car elle est très faible (<2.2e-16) , comme vous pouvez le voir ci dessous :
Anova(my_lme)
Analysis of Deviance Table (Type II tests)
Response: score
Chisq Df Pr(>Chisq)
Machine 175.6 2 < 2.2e-16 ***Bonne continuation
Bonjour,
Merci pour cette article. J’ai appliqué une anova non paramétrique en utilisant le package lmperm. Mais maintenant je suis bloqué car je n’arrive pas à comparer les différents groupes entre eux, comme la fonction aovp() ne renseigne que sur la significativité en général. Existe t il des test de comparaisons multiples style TukeyHSD ou autre pour l’ANOVA en permutation ? Merci d’avance
Bonjour,
une piste consiste peut être à utiliser le package pgrimess, vous trouverez un exemple de tests post-hoc à la fin de cette page : https://dmwiig.net/2015/03/24/using-r-in-nonparametric-statistical-analysis-the-kruskal-wallis-test-part-three-post-hoc-pairwise-multiple-comparison-analysis-of-ranked-means/
Bonne continuation
Bonjour, Je travaille actuellement sur des modèles linéaire pour comparer l’effet d’un environnement sur les comportement d’individus , j’ai donc deux groupes ayant eu des environnements différents . J’ai voulu comparer ces deux groupes à l’aide de mes modèles linaires seulement chez des effectifs très différents entre les groupes : N = 5 et N=15. Je voudrais donc faire un test de permutation pour palier à ce problème mais en essayant avec PermTest je me retrouve avec le message suivant :
Error in UseMethod(“PermTest”) :
no applicable method for ‘PermTest’ applied to an object of class “c(‘lmerModLmerTest’, ‘lmerMod’, ‘merMod’)”
Il semble donc que PermTest n’st pas applicable pour mon modèle alors que j’ai utilisé lme comme dans votre exemple… Avez vous une explication à mon problème ?
Merci d’avance
Bonjour,
Est ce que vous parvenez à reproduire l’exemple de l’article ?
Je n’ai pas de réponse immédiate, si votre problème persiste, essayez de contacter le développeur du package, vous trouverez ses coordonnées sur la première page de la documentation du package.
Bonne continuation
Bonjour.
Merci pour le post. Je m’intéresse au test de permutation à la suite des modèles linéaires à effets mixtes avec PermTest du packages pgirmess.
Je voudrais s’il y’a un test post hoc applicable après avoir exécute le test de permutation. Quoi faire à la suite du test de permutation?
Merci.
Je ne connais pas précisément la méthode à employer avec un modèle mixte , mais vous pouvez peut être vous inspirer de cette page : https://rcompanion.org/rcompanion/d_06a.html
Bonne continuation
Bonjour Madame
Je souhaite apporter une petite précision à cet article. J’ai souhaité réaliser des PERMANOVA dans le cadre d’un travail sur les biomasses végétales le long d’un gradient terre-mer. Les stations étaient disposées le long de ce continuum de manière empirique (i.e. non-aléatoirement) afin d’évaluer l’effet des types de végétation sur les biomasses obtenues. La variable d’intérêt est donc la biomasse totale et la variable explicative est donc le type de végétation avec 7 niveaux (sous-types de types de végétation). Si je me refère à votre article, je peux donc employer la fonction lmperm::aovp() pour réaliser un test de permutation qui s’écrirait ainsi : mod <- aovp(BMTOT ~ Typedevegetation, perm="Exact ou Prob", data="X") et faire les comparaisons 2 à 2 avec la fonction glht du package multCompen procédure Tukey de cette manière : Comp <- glht(mod,linfct=mcp(Typedevegetation="Tukey")) et obtenir les lettres de significativité de cette manière multComp::cld(Comp).
J'observe que les résultats sont très proches de ce que j'obtiens avec :
(1) une anova classique avec un léger défaut de la normalité des résidus (mais shapiro.test est très sensible aux écarts quand on a beaucoup de données et le QQ-plot est globalement très bon, hormis 3 "outliers")
(2) un Kruskall-Wallis suivi d'un test de Dunn.
La procédure de la PERMANOVA révèle des différences entre des sous-groupes qui semblaient diverger mais qui pouvaient ne pas être considérés comme différents par les ANOVA classiques et la procédure non paramétrique (peut-être pour une raison de taille d'échantillons). Ceci a tendance à me faire croire que je peux avoir confiance envers les résultats de ma PERMANOVA qui a pu, en faisant des permutations, révéler qu'il y a bien une différence entre les sous-groupes qui étaient "marginalement" significatifs.
Toutefois, à la lecture de l'article de M. Anderson (2001), il apparait que dans une telle situation où les traitements existent par défaut dans la nature la procédure assume que les variances sont identiquement et indépendamment distribuées au sein des groupes. Cette hypothèse ne semble pouvoir être écartées que dans le cas d'un positionnement aléatoire des stations d'échantillonnage ou dans le cadre d'une expérimentation en laboratoire, en condition contrôlée, avec une assignation aléatoire des individus dans les différentes classes. (source : DOI: 10.1139/cjfas-58-3-626 partie "Background and rational for permutation tests").
De plus, M. Anderson (2014), encore lui, précise que la PERMANOVA est très robuste face à l'hétérogénéité SEULEMENT dans le cadre de plans équilibrés (i.e. même nombre de réplicats par groupe) et que le choix de la métrique de distance (e.g. Bray-Curtis, Euclidienne, etc) est de fait très importante. (Source : DOI: 10.1002/9781118445112.stat07841 partie 3.2)
Il en ressort que la fonction lmperm::aovp() ne semble pas être adéquate dans bien des situations (à minima en écologie) et qu'elle peut nous faire conclure des différences à tort. La lecture des deux articles d'Anderson est intéressante pour bien comprendre les subtilités sous-jacente aux PERMANOVA et aux autres tests de permutations (dans le cas de l'article de 2001). Il doit exister plus de biblio sur le sujet mais je n'ai pas le niveau suffisant en mathématiques et en statistiques pour tout bien comprendre. Je recommande également plutôt l'utilisation du package "vegan" et des fonctions "adonis()" et "permutations()" pour réaliser ce genre de test puisque nous pouvons choisir la métrique de distance et ensuite travailler sur l'objet créer pour réaliser des nMDS. Je recommande également la lecture des Mat&Met d'articles d'écologie qui paraissent dans les grosses revues et qui présentent souvent des SI détaillant toutes les procédures utilisées avec les packages respectifs (voir par exemple Rovai et al. 2018).
Les tests de permutations ne sont pas des procédures si simples en fin de compte. Elles font entrer plus de subtilités que je ne le pensais alors si vous travaillez sur des données écologiques avec des plans déséquilibrés et des sites qui ne sont pas disposé aléatoirement, il convient de considérer les quelques mises en garde formulées plus haut.
Bien à vous.
Bonjour,
Un grand Merci pour ce retour qui est très intéressant ! Je vais essayer d’aller lire les deux papier d’Anderson.
Bien à vous.
Bonjour,
Merci pour votre blog qui est souvent très intéressant. Cependant, je vais également aller dans le sens de Corto : mes recherches suggèrent que les tests de permutation sont plus compliqués que cela. Par exemple, si l’on s’en sert pour tester des “égalités de moyennes”, alors ils supposent globalement les mêmes prérequis que les tests plus “classiques”, qu’ils soient paramétriques ou non (les tests non-paramétriques faisant eux-aussi des suppositions concernant la variance des groupes). C’est notamment indiqué ici : https://en.wikipedia.org/wiki/Permutation_test#Limitations, ou encore dans cet article de Huang et al. (2006) : https://doi.org/10.1093/bioinformatics/btl383, qui suggère de ne PAS utiliser ces tests pour tester des égalités de moyenne.
En conséquence, je pense que vous devriez mettre à jour votre article de blog pour prévenir les usagers (quand vos données sont complexes et que vos effectifs sont très faibles, il n’y a aucune solution miracle). Votre idée sur l’automatisation des tests répétés (dernier paragraphe) est du coup impossible, il faut quand même regarder les distributions des données avant d’employer ces tests.
Cordialement,