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

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")