Les tests post hoc sous R

Les tests post-hoc se divisent entre deux objectifs : être puissant ou être conservateur.

• Les tests conservateur évitent les erreurs de type 1 en évitant d'écarter l'hypothèse nulle à tort en voyant de fausses différences significatives.

• Les tests puissants évitent les erreurs de type 2 en évitant de passer à côte de différences significatives.

Une comparaison entre les tests les plus classiques confirme bien que Newan-Keuls est puissant (plus d'erreur de type1) et que Tukey ou Scheffé sont conservateurs (plus d'erreur de type 2). En additionnant les erreurs (type 1 et type 2) en fonctions de différents échantillons testés (plus ou moins variables), on voit que le test de Newman-Keuls fait globalement moins d'erreur.

Une remarque intéressante, un bootstrap sur moyenne avec la fonction pairwise.boot(mu="mean") de {KefiR} donne de meilleurs résultats par simple boostrap que l'ensemble des tests (type 1 et 2 cumulé). Il est puissant,  et s'il se trompe plus sur les données qui varient peu (faible écart-type), il s'en sort mieux dans les conditions où les données varient beaucoup.

• Comment faire un Newman-Keuls  ?

data$interaction <- interaction(Lycée,Sexe)

model <- aov(Note~interaction,data=data)

library(agricolae)

print(SNK.test(model,"interaction"))

• Comment faire un Tukey ?

HSD.test(model, "interaction", group = TRUE)$groups

• Un bout de code qui fait de la comparaison mutiple et renvoie des groupes

library(multcomp)

m <- lm(Note~interaction, data = data)

r <- glht(m, linfct = mcp(interaction = "Tukey"))

summary(r)

mult<- summary(glht(m, linfct = mcp(interaction = "Tukey")), test = adjusted("holm"))

mult

cld(mult)

• Faire un test de Dunnett (comparaison à un contrôle)

library(DescTools)

DunnettTest(Note ~ Lycée, data = data,control="B")

On ajoutera que d'autres tests sont à creuser : LSD (faible nombre de groupes), Dunnett (comparaison à un contrôle), Duncan (comme Scheffé : sur un petit nombre de groupes aux variances pouvant être différentes - travaille en revanche sur les étendus). Source : IWV93009.pdf (univ-irem.fr) 

D'autres, enfin, après avoir travaillé sur les moyennes marginales avec emmeans() de {emmeans}, préféreront utiliser la fonction emmeans_test() de {rstatix} pour établir des comparaison paire à paire.

emmeans_test(data,Note~interaction, p.adjust.method = "bonferroni")

# ou encore

pairs(emmeans(model,"Lycée"))

pwpm(emmeans(model,"Lycée"))

model <- aov(Note~Lycée+Sexe,data=data)

resultat <- emmeans(model,~Lycée+Sexe, adjust = "bonferroni")

# Création d'une matrice de moyennes estimées pour chaque combinaison de 'Lycée' et 'Sexe'

matrice <- xtabs(emmean ~ Lycée + Sexe, data = resultat)

# Création de matrices pour les limites inférieures et supérieures des intervalles de confiance

matriceIC_l <- xtabs(lower.CL~ Lycée + Sexe, data = resultat)

matriceIC_u <- xtabs(upper.CL~ Lycée + Sexe, data = resultat)

# Création d'un diagramme à barres avec les moyennes estimées

barplot(matrice,beside=T,col=rainbow(4))

# Chargement de la bibliothèque 'gplots'

library(gplots)

# Création d'un diagramme à barres avec les intervalles de confiance

mybar <- barplot2(matrice,beside=T,col=rainbow(4),

plot.ci=T,ci.l=matriceIC_l,ci.u=matriceIC_u,ylim=c(0,20))

# Ajout d'une grille au graphique

grid()

# Réalisation d'une analyse de variance (ANOVA) sur les notes en fonction de l'interaction entre 'Lycée' et 'Sexe'

data$interaction <- interaction(Lycée,Sexe)

model <- aov(Note~interaction,data=data)

# Réalisation d'un test post-hoc HSD (Honest Significant Difference) sur le modèle ANOVA

HSD.test(model, "interaction", group = TRUE) -> myhsd

# Extraction des interactions à partir des résultats du test HSD

mes_interactions <- rownames(myhsd$groups)

# Chargement de la bibliothèque 'stringr'

library(stringr)

# Séparation des interactions en deux facteurs

str_split_fixed(mes_interactions, "\\.",2)->mes_fac

# Tri des interactions en fonction du deuxième facteur puis du premier facteur

order(mes_fac[,2])->ordre # Trier les valeurs de sexe FFFFMMMM

ordre[order(mes_fac[ordre,1])]->ordre # Trier les valeurs de Lycée AABBCCDD

ordre[order(mes_fac[ordre,2])]->ordre # Retrier en fonction du sexe pour avoir FA,FB,FC,FD,MA,MB,MC,MD

rownames(myhsd$groups)[ordre] # contrôler le tri

text(mybar,matriceIC_u+1,myhsd$groups$groups[ordre]) # Afficher les lettres de significativité

On peut aussi indiquer les p-values sur des graphiques réalisés avec {ggplot2} grâce à la suite {ggpubr}. 

Sources externes  : 

• ANOVA dans R: Excellente Référence - Datanovia 

• ANOVA à 2 facteurs : principe - DellaData