Comparaisons multiples et ajustement des p-values avec le logiciel R
Dans un article précédent, je parlais des comparaisons multiples qui sont réalisées après à une ANOVA significative, afin d’identifier les moyennes statistiquement différentes entre elles.
Je disais que lorsque des comparaisons multiples sont réalisées, le risque de se tromper (en déclarant que deux moyennes sont différentes si en réalité, elles ne le sont pas) n’est plus contrôlé, comme dans le cas d’un test unique, mais qu’il augmente avec le nombre de tests.
Dans cette situation, il est souvent admis, qu’il est nécessaire d’ajuster (c’est à dire de relever ) les p-values obtenues pour chacun des multiples tests, afin que le risque de se tromper, pour l’ensemble de ces tests, soit à nouveau contrôlé, à un niveau souhaité.
Les comparaisons multiples sont très répandues en biostatistiques, on les rencontre dans les ANOVA, mais également dans les régressions linéaires multiple (et ses dérivées : modèles linéaires généralisés, modèles linéaires à effets mixtes, etc…) puisque se sont de multiples tests qui sont réalisées pour évaluer la significativité de multiples prédicteurs inclus dans ces modèles.
La semaine dernière je suis revenue sur cette problématique des comparaisons multiples et de l’ajustement des pvalues avec deux étudiants, alors j’ai pensé que ça pourrait intéresser certains d’entre vous de faire le point sur ce sujet.
Dans cet article je vais donc essayer de vous expliquer pourquoi le risque de se tromper augmente lorsqu’on augmente le nombre de tests, quelles sont les approches pour y remédier mais sans entrer dans leur détail, et surtout comment les mettre en œuvre avec le logiciel R. Avant de rentrer dans le vif du sujet, je vais faire un petit rappel sur les tests d’hypothèses.
Rappels sur les tests d’hypothèse
Statistique du test et hypothèses
Tous les tests d’hypothèses, c’est-à-dire les tests statistiques classique ont :
- une statistique : c’est une variable calculée avec les données observée
- et deux hypothèses : H0 l’hypothèse nulle et H1 l’hypothèse alternative.
Par exemple, la statistique du test t de Student employé pour comparer deux moyennes est :
\[T = \frac{m_1 – m_2}{s_p\star\sqrt{{\frac{1}{n_1}}+\frac{1}{n2} }}\]
Avec :
\[s_p=\sqrt\frac{(n_1-1)s_1²+(n_2-1)s_2²}{n_1+n_2-2}\]
Où :
- n1 est la taille de l’échantillon 1
- n2, la taille de l’échantillon 2
- s1² est la variance estimée de l’échantillon 1
- s2² la variance estimée de l’échantillon 2.
Les hypothèses du test t de Student sont :
\[H0 : m_1 – m_2 =0\]
\[H1 : m_1-m_2 \ne 0\]
Distribution de la statistique du test sous l’hypothèse nulle
Sous l’hypothèse nulle, la statistique du test suit une distribution théorique. Par exemple, sous l’hypothèse nulle que les moyennes ne sont pas différentes, la statistique T du test t suit une distribution de Student.
Sur l’axe des abscisses, se trouvent toutes les valeurs possibles que peut prendre la statistique T, sous l’hypothèse H0 que les moyennes sont identiques. L’axe des ordonnées correspond aux probabilités de ces différentes valeurs. Ainsi, plus la densité est forte, plus la probabilité de la valeur de T correspondante est élevée.
L’hypothèse nulle du test T spécifie que les deux moyennes ne sont pas différentes, autrement dit que leur différence est nulle. Cette différence de moyennes se retrouve au numérateur de la statistique T, c’est pour cela que la courbe de densité est centrée sur 0 et que la probabilité est la plus forte pour cette valeur.
NB : La distribution de Student dépend du nombre de données. On dit que la statistique T suit une distribution de Student à n1+n2-2 degrès de liberté.
Règle de décision et risque alpha
Qui dit test, dit règle de décision ! Il faut bien une méthode pour décider si les moyennes des échantillons observées sont différentes ou non, compte tenu de la dispersion des données autour de ces moyennes.
C’est à cela que sert la statistique du test. En connaissant sa distribution théorique sous l’hypothèse nulle (que les moyennes ne sont pas différentes), cela signifie que l’on connait toutes les valeurs qu’elle peut prendre, avec les probabilités correspondantes. On va alors fixer une limite dans ces valeurs. Pour cela, on va choisir une valeur qui a une très faible probabilité de se produire ; en pratique 5%.
Et comme l’hypothèse alternative du test est très souvent bilatérale (cela veut dire que l’on ne présume pas que m1 > m2 ou que m1 < m2) alors on répartie généralement les 5% de chaque côté de la courbe de densité de probabilité. Au final, on rejette l’hypothèse d’égalité des moyennes chaque fois que la valeur absolue de la statistique T calculée avec les données observées est supérieure à la valeur théorique de cette statistique correspondant à une probabilité de 2.5% !
Cette probabilité de 5%, qui correspond donc à 2.5% à droite et à gauche de la courbe de densité de probabilité, est appelé risque alpha. Il correspond à la probabilité de se tromper, en concluant à la différence significative des moyennes, si en réalité elles ne le sont pas. Ce type d’erreur est appelé “faux positif”. Pour exprimer les choses autrement : même si les moyennes ne sont pas différentes, la statistique T calculée peut prendre des valeurs extrêmes. Mais si elle dépasse une limite au-delà de laquelle seulement 2.5% des statistiques théoriques sont attendues, alors l’hypothèse d’égalité est rejetée, et les deux moyennes sont considérées comme significativement différentes.
La pvalue
En pratique, lorsqu’on réalise un test statistique avec un logiciel, comme R par exemple, on ne se soucie pas directement de la valeur seuil. En revanche, on s’intéresse de près à la p-value fournie par le logiciel. Cette p-value correspond à la probabilité de la statistique calculée (avec les données observées) sur la courbe théorique définie sous H0. Dit plus simplement, la p-value correspond à la probabilité, si les moyennes ne sont pas différentes, d’observer une valeur de T égale à la statistique calculée avec les données mesurées. D’un point de vue pratique-t-il s’agit de l’aire sous la courbe de densité, au-delà de la valeur de la statistique T calculée avec les données observée.
Cette p-value est ensuite comparée au risque alpha que l’on s’est fixé avant le test (il est quasiment toujours égal à 5%) :
- si la probabilité est inférieure au seuil alpha que l’on s’est fixé, cela signifie que la statistique T calculée se situe au-delà du seuil de décision. On conclura alors au rejet de l’hypothèse nulle et donc au fait que les moyennes sont significativement différentes.
- au contraire, si la probabilité est supérieure au seuil alpha que l’on s’est fixé, cela signifie que la statistique T se situe avant le seuil de décision. On conclura alors à l’acceptation de l’hypothèse nulle, et donc au fait que les moyennes ne sont pas significativement différentes.
Les comparaisons multiples et l’augmentation de la proabbilité de se tromper
Imaginons que nous avons comparé globalement 3 moyennes (dont un contrôle) à l’aide d’une ANOVA, que celle-ci est significative, et qu’à présent, nous souhaitons comparer les deux moyennes à la moyennes contrôle pour savoir si elles sont différentes. Nous allons donc réaliser deux tests post-hocs.
Comme nous l’avons vu précédemment, lorsqu’on ne considère qu’un seul test, la probabilité de faire une erreur en concluant que les deux moyennes sont différentes alors qu’elles ne le sont pas , autrement dit de faire un faux positif, correspond au risque alpha:
\[\text {P(faire 1 erreur)} = {\alpha }\]
La probabilité complémentaire, correspond donc à la probabilité de ne pas commettre d’erreur :
\[\text {P(NE PAS faire d’erreur)} = 1- {\alpha } \]
Lorsque deux tests sont réalisés, la probabilité de ne pas faire d’erreur (de type faux positif) lors du premier test ET lors du second test, est la probabilité de ne commettre aucune erreur. En pratique cette probabilité de ne pas faire d’erreur lors des deux tests est le produit de la probabilité de ne pas faire d’erreur lors d’un test unique :
\[ \text {P(NE PAS faire d’erreur dans les 2 tests)} = (1- {\alpha }) {\star}(1- {\alpha }) \]
\[\text {P(NE PAS faire d’erreur dans les 2 tests)} = (1- {\alpha })^2\]
Ne pas faire d’erreur lors des deux tests signifie “ne faire aucune erreur lors des 2 tests”, On peut donc aussi écrire :
\[\text {P(faire aucune erreur lors des 2 tests)} = (1- {\alpha })^2\]
La probabilité complémentaire de cette probabilité correspond, elle à la probabilité de commettre au moins une erreur, de type faux positif, lors des deux tests. On peut l’écrire :
\[\text {P(faire au moins une erreur lors des 2 tests)} = 1-(1- {\alpha })^2\]
Lorsqu’on réalise k comparaisons multiples, l’équation se généralise en :
\[\text {P(faire au moins une erreur lors de k tests)} = 1-(1- {\alpha })^k \]
Cette probabilité est appelée risqua alpha global:
\[ {\alpha}_{global} = 1-(1- {\alpha })^k \]
Ce qu’il faut retenir ici c’est que le risque de commettre au moins une erreur (de type faux positif) augmente lorsque le nombre de comparaisons multiples augmente. Avec une simulation, nous pouvons obtenir le graph ci-dessous :
alpha = 0.05
k <- 1:100
Pr <- 1-(1-alpha)^k
mydf <- data.frame(k=k, Pr=Pr)
library(ggplot2)
ggplot(mydf, aes(x=k, y=Pr))+
geom_point(colour=c("#66CDAA"), size=2)+
theme_classic()+
ylab("Risque alpha global")+
xlab("Nombre de tests réalisés") 
C’est assez marquant, non ?
Imaginons qu’on fasse une ANOVA pour comparer 4 moyennes, et qu’ensuite on veuille faire toutes les comparaisons 2 à 2, ça fait 3+2+1, soit 6 tests. Eh bien, si on fixe le risque de se tromper pour un test unique à 5%, alors le risque de se tromper au moins une fois sur les 6 tests n’est plus de 5% mais de 26 % !
On est donc tenté de faire quelque chose !
Les approches pour prendre en compte l’augmentation du risque global de se tromper
On peut diminuer le risque alpha pour qu’ensuite le risque alpha global soit égal à 5%. Mais cela n’est pas très pratique puisque le niveau du risque alpha dépend du nombre de tests à réaliser dans chaque situation.
L’autre possibilité est d’agir au niveau des p-values. Il s’agit, en quelque sorte, de les recalibrer, en pratique de les rehausser, pour pouvoir continuer à les comparer dans tous les cas, à un risque alpha de 5%. C’est cela qu’on appelle l’ajustement des p-values !
Sans trop entrer dans les détails, il existe, grosso modo, deux types d’approches :
- celles qui contrôlent le risque alpha global de se tromper au moins une fois ; on l’appelle le “Family-wise error rate” (FEWR) en anglais
- celles qui contrôle le risque alpha global de se tromper au moins une fois, mais en ne considérant que les tests qui ont rejetté H0 ; On l’apelle le “False discovery rate” (FDR) en anglais.
Les méthodes qui contrôlent le Family-wise error rate (FEWR)
On peut diminuer le risque alpha pour qu’ensuite le risque alpha global soit égal à 5%. Mais cela n’est pas très pratique puisque le niveau du risque alpha dépend du nombre de tests à réaliser dans chaque situation.
L’autre possibilité est d’agir au niveau des p-values. Il s’agit, en quelque sorte, de les recalibrer, en pratique de les rehausser, pour pouvoir continuer à les comparer dans tous les cas, à un risque alpha de 5%. C’est cela qu’on appelle l’ajustement des p-values !
Sans trop entrer dans les détails, il existe, grosso modo, deux types d’approches :
- celles qui contrôlent le risque alpha global de se tromper au moins une fois ; on l’appelle le “Family-wise error rate” (FEWR) en anglais
- celles qui contrôlent le risque alpha global de se tromper au moins une fois, mais en ne considérant que les tests qui ont rejetté H0 ; On l’apelle le “False discovery rate” (FDR) en anglais.
Les méthodes qui contrôlent le Family-wise error rate (FEWR)
Là encore, on peut les diviser en 2 catégories :
- les approches “single step” qui corrigent toutes les p-values de la même façon. La plus connue est la méthode de Bonferroni. Pour plus de détails sur la correction vous pouvez consulter cette page)
- les approches séquentielles avec lesquelles la correction s’adapte à la valeur de la p-value, la plus connue est la méthode d’Holm. Là encore, pour plus de détails sur la méthode vous pouvez consulter ce lien (le calcul est décrit dans la partie Extension).
Les méthodes qui contrôlent le False discovery rate
Les plus connues sont les approches de Benjamini et Hochberg. Elles sont parfois préférées car moins conservatives (elles augmentent moins les p-values) que les approches qui contrôlent le FEWR.
Quelle méthode choisir ?
Il est plutôt admis que la méthode de Bonferronni est trop restrictive et que les autres méthodes sont préférables. En écologie la méthode d’Holm est souvent utilisée. De mon point de vu, le choix de la méthode d’ajustement n’a pas beaucoup d’importance (en dehors de la méthode de Bonferronni), le plus important étant d’en appliquer une !
Néanmoins, si l’ajustement des pvalues est généralement réalisé dans le cas des tests post-hoc après ANOVA, son utilisation dans les modèles linéaires complets ne fait pas consensus, du moins pas en écologie. C’est pourquoi on trouve souvent des publications qui fournissent à la fois les values brutes et les pvalues ajustées.
Ajustement des pvalues avec le logiciel R
Il y a deux possibilités, soit vous ajustez vous- même les pvalues, en dehors de l’ANOVA ou de la régression linéaire multiple avec la fonction p.adjust, soit vous utilisez un package qui le fait pour vous, comme le package multcomp.
Ajuster les p-values avec la fonction p.adjust
La fonction “p.adjust” du package “stats” (qui est chargé par défaut) propose 7 méthodes d’ajustement : holm, hochberg, hommel, bonferroni, Benjamini & Hochberg , Benjamini & Yekutiel et fdr (qui semble indentique à l’approche Benjamini & Hochberg)
Imaginons que vous disposez d’un vecteur de pvalues, et que vous voulez les ajuster. Ici, j’utilise le jeu de données “Robey” du package “car” et je compare toutes les moyennes de la variable “tfr” deux à deux, en utilisant des tests t de Student :
options(digits = 3, scipen = -5)
library(car)
mc <- pairwise.t.test(Robey$tfr,Robey$region)
pval_brutes <- mc$p.value
pval_brutes
## Africa Asia Latin.Amer
## Asia 2.93e-05 NA NA
## Latin.Amer 1.49e-04 3.48e-01 NA
## Near.East 1.63e-01 1.47e-01 3.48e-01 Il suffit alors d’utiliser la fonction p.adjust en utilisant comme argument les pvalues brutes, ainsi que la méthode d’ajustement désirée. Pour plus de détails sur les approches disponibles, utilisez l’aide : ?p.adjust
pval_ajust <- p.adjust(pval_brutes, "holm")
pval_ajust <- matrix(pval_ajust , nc=3) # met les p-values sous le forme d'une matrice.
pval_ajust
## [,1] [,2] [,3]
## [1,] 1.76e-04 NA NA
## [2,] 7.47e-04 6.96e-01 NA
## [3,] 5.89e-01 5.89e-01 6.96e-01 Dans le cas d’une régression linéaire multiple :
options(digits = 3, scipen = -5)
library(car)
rlm <- lm(prestige~., data=Prestige)
summary(rlm)
##
## Call:
## lm(formula = prestige ~ ., data = Prestige)
## ## Residuals:
## Min 1Q Median 3Q Max
## -1.30e+01 -4.98e+00 6.98e-01 4.87e+00 1.92e+01
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.21e+01 8.02e+00 -1.51e+00 1.3e-01
## education 3.93e+00 6.53e-01 6.02e+00 3.6e-08 ***
## income 9.95e-04 2.60e-04 3.82e+00 2.4e-04 ***
## women 1.31e-02 3.02e-02 4.30e-01 6.7e-01
## census 1.16e-03 6.18e-04 1.87e+00 6.5e-02 .
## typeprof 1.08e+01 4.68e+00 2.30e+00 2.4e-02 *
## typewc 2.88e-01 3.14e+00 9.00e-02 9.3e-01
## ---
## Signif. codes: 0e+00 '***' 1e-03 '**' 1e-02 '*' 5e-02 '.' 1e-01 ' ' 1e+00
##
## Residual standard error: 7.04e+00 on 91 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.841, Adjusted R-squared: 0.831
## F-statistic: 80.2 on 6e+00 and 91 DF, p-value: <2e-16 On récupère alors les p-values non ajustés, comme cela :
rlm_pval_brut <- summary(rlm)$coeff[,4]
rlm_pval_brut
## (Intercept) education income women census typeprof
## 1.34e-01 3.64e-08 2.40e-04 6.65e-01 6.47e-02 2.35e-02
## typewc
## 9.27e-01 Puis on les injecte dans la fonction p.adjust() :
rlm_pval_ajust <- p.adjust(rlm_pval_brut, "fdr")
rlm_pval_ajust
## (Intercept) education income women census typeprof
## 1.87e-01 2.55e-07 8.40e-04 7.76e-01 1.13e-01 5.49e-02
## typewc
## 9.27e-01 Utiliser un package pour ajuster les pvalues
J’utilise généralement le package multcomp. Les méthodes d’ajustement utilisées par ce package sont un peu complexes, je ne vais pas entrer dans le détail. On peut juste retenir que tout en étant basées sur les procédures classiques de Bonferroni, Holm, Hochberg etc.., elles sont plus puissantes car elles prennent en compte la corrélation entre les tests. Les pvalues ajustées par les procédures du pcakages sont alors moins corrigées, c’est à dire plus faibles que celles obtenues par la fonction p.adjust
La fonction du package multcomp qui permet d’obtenir des p-values ajustées est la fonction glht.
Voici un exemple de l’ajustement des p-values pour les comparaisons multiples après une ANOVA. POur plus de détails, vous pouvez consulter cet article.
library(multcomp)
my_anova <-lm(response~trt, data=cholesterol)
mc_tukey <- glht(my_anova, linfct=mcp(trt="Tukey"))
summary(mc_tukey)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = response ~ trt, data = cholesterol)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 2times - 1time == 0e+00 3.44e+00 1.44e+00 2.39e+00 1.4e-01
## 4times - 1time == 0e+00 6.59e+00 1.44e+00 4.57e+00 3.7e-04 ***
## drugD - 1time == 0e+00 9.58e+00 1.44e+00 6.64e+00 < 1e-04 ***
## drugE - 1time == 0e+00 1.52e+01 1.44e+00 1.05e+01 < 1e-04 ***
## 4times - 2times == 0e+00 3.15e+00 1.44e+00 2.18e+00 2.0e-01
## drugD - 2times == 0e+00 6.14e+00 1.44e+00 4.25e+00 9.7e-04 ***
## drugE - 2times == 0e+00 1.17e+01 1.44e+00 8.12e+00 < 1e-04 ***
## drugD - 4times == 0e+00 2.99e+00 1.44e+00 2.07e+00 2.5e-01
## drugE - 4times == 0e+00 8.57e+00 1.44e+00 5.94e+00 < 1e-04 ***
## drugE - drugD == 0e+00 5.59e+00 1.44e+00 3.87e+00 3.1e-03 **
## ---
## Signif. codes: 0e+00 '***' 1e-03 '**' 1e-02 '*' 5e-02 '.' 1e-01 ' ' 1e+00
## (Adjusted p values reported -- single-step method) La méthode utilisée par défaut est une méthode single step, mais vous pourrez constater que les p-values obtenues sont légèrement plus faibles que celles obtenues avec la fonction p.adjust() employée avec la procédure de Holm.
On peut également utiliser la fonction glht() dans le cadre d’une régression linéaire multiple. Elle permet alors d’obtenir la pvalue ajustée du test de significativité de tous les coefficients des prédicteurs inclus dans le modèle.
rlm <- lm(prestige~., data=Prestige)
summary(rlm) # pvalues non ajustées
##
## Call:
## lm(formula = prestige ~ ., data = Prestige)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.30e+01 -4.98e+00 6.98e-01 4.87e+00 1.92e+01
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.21e+01 8.02e+00 -1.51e+00 1.3e-01
## education 3.93e+00 6.53e-01 6.02e+00 3.6e-08 ***
## income 9.95e-04 2.60e-04 3.82e+00 2.4e-04 ***
## women 1.31e-02 3.02e-02 4.30e-01 6.7e-01
## census 1.16e-03 6.18e-04 1.87e+00 6.5e-02 .
## typeprof 1.08e+01 4.68e+00 2.30e+00 2.4e-02 *
## typewc 2.88e-01 3.14e+00 9.00e-02 9.3e-01
## ---
## Signif. codes: 0e+00 '***' 1e-03 '**' 1e-02 '*' 5e-02 '.' 1e-01 ' ' 1e+00
##
## Residual standard error: 7.04e+00 on 91 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.841, Adjusted R-squared: 0.831
## F-statistic: 80.2 on 6e+00 and 91 DF, p-value: <2e-16
summary(glht(rlm)) # pvalues ajustées
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = prestige ~ ., data = Prestige) ## ## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) == 0e+00 -1.21e+01 8.02e+00 -1.51e+00 5.0e-01
## education == 0e+00 3.93e+00 6.53e-01 6.02e+00 < 1e-02 ***
## income == 0e+00 9.95e-04 2.60e-04 3.82e+00 < 1e-02 **
## women == 0e+00 1.31e-02 3.02e-02 4.30e-01 9.9e-01
## census == 0e+00 1.16e-03 6.18e-04 1.87e+00 2.9e-01
## typeprof == 0e+00 1.08e+01 4.68e+00 2.30e+00 1.2e-01
## typewc == 0e+00 2.88e-01 3.14e+00 9.00e-02 1.0e+00
## ---
## Signif. codes: 0e+00 '***' 1e-03 '**' 1e-02 '*' 5e-02 '.' 1e-01 ' ' 1e+00
## (Adjusted p values reported -- single-step method)
Pour aller plus loin
Vous pouvez aussi aller voir du côté des vignettes , en copiant collant ce code dans la console R
vignette("generalsiminf", "multcomp")
vignette("chfls1", "multcomp") J’espère que cet article permettra au plus grand nombre de comprendre l’enjeu derrière les comparaisons multiples, et de savoir ajuster les pvalues avec le logiciel R.
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 livre mentionnés cela m’aide à entretenir ce blog, merci si vous le faites ! 
Retrouvez ici 4 de mes articles les plus consultés:
