Correction des p-values en fonction du nombre de tests statistiques

en langage R

L'essentiel de cette page

Lorsque l'on multiplie les tests et surtout le même test (ex : contrôle de la normalité des échantillons), il est bon de contrôler le taux d'erreur sur l'ensemble des tests (= par famille).

Plusieurs corrections sont utilisables avec p.adjust() de R, la plus commune étant "holm" pour les tests post-hocs ou "bonferroni" qui présente l'inconvénient d'être très conservateur (il passe à côté des différences significatives).

Lorsque l'on fait des tests indépendants, contrôle des variances ou des normalités, la correction de Sidak est très pertinente (plus puissante, moins conservatrice), mais nécessite la fonction Sidak.p.adjust() du package {MHTdiscrete}.

Holm

Ordre les p-values par ordre croissant, puis pour chaque p-value, multiplie par le nombre de tests restants et compare à alpha. Si la p-value est supérieure à alpha, arrête le processus. Sinon, continue avec la prochaine p-value.

Très conservateur (mais moins que Bonferroni)

Bonferroni

Multiplie chaque p-value par le nombre total de tests et compare à alpha.

Très conservateur

Sidak

Multiplie chaque p-value par (1 - (1 - alpha)^(1/m)) and compare to alpha 2.

Puissant et bien adaptés aux échantillons indépendants

Hommel

Ordre les p-values par ordre croissant, puis pour chaque p-value, multiplie par le nombre total de tests et compare à alpha/(m+1-i). Si la p-value est supérieure à ce seuil et toutes les autres p-values précédentes sont également supérieures à leur seuil respectif, arrête le processus. Sinon, continue avec la prochaine p-value.

Puissant

Ordre les p-values par ordre croissant, puis pour chaque p-value, multiplie par le nombre total de tests et compare à alpha*i/m. Si la p-value est supérieure à ce seuil, arrête le processus. Sinon, continue avec la prochaine p-value.

Puissant, moins conservateur que Holm et Bonferroni, mais plus lent.

Comme BH mais utilise une valeur plus conservatrice pour l’estimation du nombre de vraies hypothèses nulles.

Puissant, moins conservateur que Holm et Bonferroni, mais plus lent.

FDR

Contrôle le taux de fausses découvertes plutôt que le taux d’erreur par famille. Multiplie chaque p-value par m/i et compare à alpha.

Puissant, moins conservateur que Holm et Bonferroni, mais ne contrôle pas le taux d'erreur par famille.

Quelques exemples d'applications :

# Simulons des variables suivant la loi de poisson...

variables <- c(rpois(40, 2)+rnorm(40,0,1),rpois(40, 6)+rnorm(40,0,1),rpois(40, 7)+rnorm(40,0,1),rpois(40, 6)+rnorm(40,0,1)) ; catégories <- rep(c("A","B","C","D"),each=40)

data <- data.frame(variables,catégories)

data$catégories <- factor(data$catégories)

# Corrections de p-values de variables indépendantes (Sidak)

pval <- by(data$variables,data$catégories,function(x){shapiro.test(x)$p.value})

# Je regarde les p-values, si une est inférieur à 0.05 alors l'une variable ne suit pas la loi normale

print(pval)

# Je vais tout de même faire une correction

# Appliquer la correction de Sidak

library("MHTdiscrete")

# tout de suite mes p-values semblent moins significatives

Sidak.p.adjust(pval)

# On peut aussi comparer nos p-values brutes à un seuil corrigée avec la fonction sidak ci-dessous :

sidak <- function(p, n = length(p)) {

n_reel <- length(p)

for (i in 1:n_reel) {

p[i] <- 1-(1-p[i])^(1/n)

}

return(p)

}

sidak(0.05 , n=4) #correction pour 4 p-values (4 tests indépendants)

Et si je souhaite comparer par un wilcoxon.test() les variables une à une, il y a cette fois-ci une dépendance entre mes tests. En effetsi var2 est supérieur à var1 et var3 supérieur à var2, on s'attend bien à avoir var3 supérieur à var1.

# Une boucle qui fait manuellement des Wilcoxon 2 à 2

pval <- c()

for (i in 1:3) {

for (j in (i+1):4) {

pval <- c(pval , wilcox.test(data$variable[data$catégories==levels(data$catégories)[i]],

data$variable[data$catégories==levels(data$catégories)[j]])$p.value)

}

# Les p-values

round(pval,2)

#Les p-values ajustées

round(p.adjust(pval,"holm"),2) # Les plus values corrigées, plus conservatrices

round(p.adjust(pval,"bonferroni"),2) # Les plus values corrigées, plus conservatrices

# On peut en réalité dans cas là faire un pairwise.wilcox.test()

pairwise.wilcox.test(data$variables,data$catégories,"holm")