Paramétrer son ANOVA à deux facteurs

Un jargon à maîtriser :

• Hypothèse de sphéricité, cf. anova sur mesures répétées.

• Indépendance

• Interaction

• Contrastes, les contrastes sont utilisés pour répondre à des questions très précises telles que savoir si la différence entre deux ensembles de moyennes distincts est significative (source externe).

• Effet aléatoire

• Moyennes marginales

• Projections, les projections sont utilisées pour déterminer la contribution de chaque groupe à la variance totale des données. 

L’ANOVA (Analyse de la variance) est une méthode statistique utilisée pour comparer les moyennes de trois groupes ou plus. Lorsque vous effectuez une ANOVA, il y a plusieurs paramètres et concepts importants à prendre en compte.

1. Complexité du modèle : La complexité du modèle fait référence au nombre de facteurs et d’interactions que vous incluez dans votre modèle ANOVA. Un modèle plus complexe peut fournir des informations plus détaillées sur les effets des différents facteurs, mais il peut également être plus difficile à interpréter.

2. Assomptions : Les ANOVA reposent sur certaines assomptions, telles que la normalité des résidus, l’homogénéité des variances et l’indépendance des observations (cf. singularité du modèle). Il est important de vérifier ces assomptions avant d’interpréter les résultats de votre analyse.

3. Singularité du modèle : La singularité du modèle fait référence à la présence de corrélations entre les variables indépendantes dans votre modèle ANOVA. Cela peut se produire lorsque vous incluez des facteurs qui sont dépendants les uns des autres. La singularité du modèle peut entraîner des problèmes lors de l’estimation des paramètres du modèle et de l’interprétation des résultats. Si vous rencontrez des problèmes de singularité, vous devrez peut-être réexaminer votre modèle et exclure certains facteurs ou interactions.

4. Post-hoc : Les tests post-hoc sont utilisés pour comparer les moyennes des différents groupes après avoir effectué une ANOVA. Ils permettent d’identifier quels groupes diffèrent significativement les uns des autres. Il existe plusieurs méthodes de test post-hoc, telles que le test de Tukey, le test de Bonferroni et le test de Scheffé. Un choix (pas toujours simple) doit être fait.

5. Visualisation des résultats : Il est souvent utile de visualiser les résultats d’une ANOVA à l’aide de graphiques, tels que des diagrammes en boîte ou des graphiques en barres. Cela peut vous aider à mieux comprendre les différences entre les groupes et à communiquer vos résultats de manière claire et concise.

J'ai fait passer un examen à des filles et des garçons issus de 4 types de lycées différents. Pour chaque type, 3 établissements ont été testés.  Les élèves n'ont pas la même option selon les situations. Certains font du sport, d'autres non. Ont-ils des notes équivalentes ? Qu'est-ce que joue sur la note ?

La fonction standard aov() de base-R utilise des sommes de carrés de type I. Elle n'est donc appropriée que si vos données sont équilibrées (autant de données dans chaque catégories). Mieux vaut donc contrôler en amont.

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

summary(model)

Attention, les modèles d'aov() dépendent de nombreuses formules telles les régressions linéaires avec lm().

model <- aov(Note~Lycée*Sexe) # Comme précédant en ajoutant l'interaction Sexe-Lycée : car les individus, selon leur sexe pourraient se comporter différemment d'un lycée à l'autre. 

summary(model)

Pour comprendre la notion d'interaction - source externe

On notera que aov() met au format traditionnel "anova" une sortie de la fonction lm(). On peut donc utiliser les écritures propres à lm() avec des I() pour forcer des formules, poly()... C'est d'autant plus utile si on ajoute une variable numérique à son modèle pour faire une ancova.

Dans ce cas, il faut faire une anova de type II ou III (explications à ces sources : source 1, source 2).

Si vos données ne sont pas équilibrées, vous devez effectuer une ANOVA avec des sommes de carrés de type II ou de type III.  En général, on prend le type 3 qui garde les interactions et ne calcule plus à partir d'une moyenne de référence mais d'une moyenne centrée). L’ANOVA de type 1 évalue les effets des variables dans l’ordre où elles sont entrées dans le modèle, tandis que l’ANOVA de type 3 évalue les effets de chaque variable après avoir pris en compte l’effet de toutes les autres variables, indépendamment de l’ordre d’entrée.

Pour ce faire, vous pouvez utiliser la fonction Anova() du package {car}.

# 1 je fais mon anova avec aov() - mais je change la méthode de calcul des contrastes.

model <- aov(Note~Lycée*Sexe,

                  contrasts=list(Sexe=contr.sum,Lycée=contr.sum), # changement de contraste nécessaire

                  data=data)

# 2 je fais une anova de type 3

library(car) # penser à installer {car}

Anova(model,type=3) # Attention ! Anova avec majuscule !

Cela revient aussi à faire, plus simple car changement de contraste non nécessaire

library(emmeans)

joint_tests(model)

Remarque sur les contrastes : il faut faire ce changement pour avoir une vraie anova de type 3. Par défaut, l'on a "contr.treatment" (chaque niveau d'un facteur est comparé avec le niveau de base) alors que "contr.sum" fait en sorte que la somme des contrastes fasse 0 (chaque niveau est comparé à un niveau moyen).

On peut aussi faire une anova de ce type avec cette ligne de commande : drop1(model, .~., test="F" ) 

Sources externes pour aller plus loin : YaRrr! The Pirate’s Guide to R (bookdown.org) et https://delladata.fr/carres-de-type-iii-analyse-de-variance/

Le mieux est de faire le modèle d'ANOVA à partir de lm()

library(car) ; 

# Définir les contrastes pour les facteurs

options(contrasts = c("contr.sum", "contr.poly"))

# Pour un modèle linéaire :

model <- lm(response ~ factor1 * factor2, data = your_data)

Anova(model, type = "III")

# Vérification des résidus pour le modèle linéaire :

qqPlot(model, main="QQ Plot")

leveneTest(residuals(model) ~ factor1*factor2, data = your_data)

On peut aussi le faire sur un modèle mixte (pensons aussi à ajuster les contrastes) :

library(lme4) ; library(nlme)

# Pour un modèle linéaire mixte :

model_mixed <- lmer(response ~ factor1 * factor2 + (1|random_effect), data = your_data)

anova(model_mixed)  # Type III SS n'est pas directement spécifié ici, ajustez selon le logiciel

# Vérification des résidus pour le modèle linéaire mixte :

qqnorm(resid(model_mixed))

qqline(resid(model_mixed))

leveneTest(resid(model_mixed) ~ factor1*factor2, data = your_data)

Solution de secours 2 : On peut aussi basculer sur les ANOVA non-paramétriques (robustes) car certaines sont insensibles à l'ordre des facteurs tel le test de Schreier-Ray-Hare.

Je ne suis pas expert de ce domaine, mais il y a quelques éléments à creuser pour ne pas foncer dans les ANOVAs en pensant que tout ira facilement dans la construction du modèle et son interprétation.

Corriger le modèle en tenant compte de l'effet aléatoire lié au choix des 3 établissements pour chaque type de lycée et en ignorant le facteur Sexe.

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

Corriger le modèle en tenant compte de l'effet aléatoire lié au choix des 3 établissements pour chaque type de lycée.

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

Corriger le modèle en tenant compte de l'effet aléatoire lié au choix des 3 établissements pour chaque type de lycée, mais aussi d'éventuels interactions qui feraient que filles et garçons n'auraient pas le même niveau selon les établissements (un biais de sélection par exemple liées à des matières plus ou moins attractives pour les filles/garçons).

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

Une réflexion dans un forum de source externe : ce sujet est complexe et nécessite une véritable expertise lorsqu'on augmente le nombre de facteurs.

Imaginons que je suis la taille de plusieurs individus dans le temps. J'ai une répétition des mesures sur ces mêmes individus.

Je dois donc construire un modèle qui aurait cette forme :

model <- aov(Taille~Temps+Error(Individus/Temps), data=data) # Cela va me permettre de corriger les effets liés aux individus et les interaction individus:Temps.

summary(model) 

Dans une situation équivalente, on va pouvoir traiter chaque individu avec un modèle mixte pour établir s'il y a un effet du Temps.

Si on applique ce modèle à nos données, on va pouvoir construire le modèle suivant :

# ANOVA

# Approche 1 : voir l'effet du type de lycée en ignorant les variations aléatoires liées aux établissements

model <- aov(Note~Lycée+Error(établissement), data=data) # effet du type lycée sur la note en corrigeant l'effet lié au tirage au sort des établissement

# Approche 2 : l'effet du lycée sur la note en corrigeant en plus les interactions qui font que les établissements peuvent changer de niveau selon le type de Lycée: ici on ne voit plus d'effet lycée.

model <- aov(Note~Lycée+Error(établissement/Lycée), data=data)

# favorisons l'approche 1.

summary(model) 

Et si l'on souhaite construire un modèle mixte : il faut installer {car}, {lme4} et {lmerTest}.

library(lme4)

model <- lmer(Note~Lycée+(1|établissement), data=data) # Construiction d'un modèle mixte difficile à lire

summary(model )

library(lmerTest)

model  <- lmerTest::lmer(Note~Lycée+(1|établissement), data=data) # Mise en forme du modèle pour avoir des p-values

ranova(model) # Pour avoir un anova avec effet aléatoire et confirmer que l'établissement n'a pas d'effet. On obtient la même chose en faisant rand(model)

summary(mod_rep)

library(car)

Anova(mod_rep)  # Conversion du modèle pour avoir une anova classique et voir l'effet du lycée. 

On peut aussi se protéger en présence de données répétées en utilisant les commandes de {rstatix} suivantes :

identify_outliers() #pour identifier d'éventuelles valeurs extrêmes

anova_test() # pour tester différents types d'anova 

get_anova_table()  # pour obtenir la table d'anova/

Une invitation à creuser le sujet. Source pour creuser.

Comprendre ce que signifie des mesures répétées - source externe.

Allons plus loin sur les modèles mixtes grâce au site de DellaData.

Dans une anova à deux facteurs, on doit contrôler quelques hypotyèses pour valider le modèle. On parle d'assomptions.

Ces hypothèses sont les suivantes :

1. Normalité des résidus: les résidus (différences entre valeurs prédites et réelles) doivent être normales.

2. Homogénéité des variances : Les variances des données de chaque groupe doivent être égales.

3. Indépendance : Les observations de chaque groupe doivent être indépendantes les unes des autres et les observations au sein des groupes doivent être obtenues par un échantillon aléatoire. Exemple : la longueur des cheveux "long" "court" n'est pas un facteur indépendant du sexe, puisque que les "F" ont tendance a les avoir plus longs que "M".

Si ces hypothèses ne sont pas vérifiées, on doit alors faire des ANOVA non paramétriques.

Attention, les conditions de l’ANOVA pour mesures répétées diffèrent un peu :

• Normalité des résidus

• Homogénéité des variances et des co-variances (puisque les observations ne sont plus indépendantes puisqu'on mesure plusieurs fois les individus).

Si l'on cite le site de DellaData :

L’hypothèse d’homogénéité des variances et des co-variances se nomme aussi hypothèse de Sphéricité. Elle est évaluée par le test de Mauchly. 

En cas de défaut de sphéricité, les résultats doivent être corrigés par les corrections de Huynh-Feldt ou de Greenhouse-Geisser.

Attention, lorsque  l'on réalise des modèles anova avec effet aléatoire (Error()), on ne peut obtenir les graphiques de validation automatiquement, ni faire de prédiction avec predict(), ni extraire simplement les résidus avec la fonction resid() ou en ajoutant $residuals à la suite du modèle. Il faut donc faire appel à la librairie {dae}.

library(dae)

residuals.aovlist(model) # Pour récupérer les résidus d'un modèle à effet aléatoire.

fitted(model) # Pour récupérer les données prédites (fitted) d'un modèle à effet  aléatoire.

Astuce : si nous n'avons pas de modèle avec effet aléatoire, il est possible de contrôler rapidement son modèle avec la librairie {performance} ou de suivre le tuto ci-après.

library(performance)

check_model(model)  ; check_normality(model) ; check_homogeneity(model)

Pour construire "manuellement" ce graphique (si par exemple, l'on a un modèle avec effet aléatoire).

library(dae)

plot(fitted(model),residuals.aovlist(model))

lines(smooth.spline(fitted(model),residuals.aovlist(model),df=10,cv=F),col="red",lwd=2) 

Contrôler la normalité des résidus de mon modèle peut se faire par plusieurs approches différentes et complémentaires :

Pour construire "manuellement" ce graphique si l'on a un modèle avec effet aléatoire :

residus <- residuals.aovlist(model)

residus_std <- residus/sd(residus)

qqnorm(residus_std) ; qqline(residus_std)

Remarque : les résidus standardisés sont calculable automatiquement à partir d'un modèle avec la fonction rstandard() , mais qui, une fois de plus, ne marche pas sur les modèles avec effet aléatoire.

Si on veut se donner le sentiment d'objectivité, on peut faire un test de contrôle de la normalité.

shapiro.test(model$residuals)$p.value

Attention, si on a plus de 50 résidus, il faut envisager de faire un test de Jarque-Bera ou de Kolmogorov-Smirnov :

# Kolmogorov-Smirnov - peut planter si redondance des résidus

ks.test(model$residuals,"pnorm")$p.value 

# Jarque-Bera - nécessite d'installer {tseries}

library(tseries) 

jarque.bera.test(model$residuals)$p.value

Pour tous ces tests, si la p-value reste au-dessus de 0,05 (seuil à nuancer en regardant les graphiques) alors on peut considérer une distribution normale des résidus.

On peut aussi contrôler l'homogénéité des variances par des tests statistiques tels que :

• Le test de Bartlett : si les données de chaque catégorie ont une distribution normale

• Le test de Levene ou le test de Brown-Forsythe (si les données ne sont pas normales).

• On parle aussi du test robuste de Fligner-Killeen en dernier espoir pour voir vers quelle anova paramétrique se tourner.

1. Vérifier que mes échantillons sont normaux

# Dans cet exemple on va croiser les catégories de la variable 2 (Lycée) et 4 (Sexe) pour voir si les notes suivent la loi normale.

# Remarque : attention, on fait ici un Shapiro, mais il faut penser au test de Jarque-Bera si on a plus de 50 notes par catégories.

by(data$Note, data[,c(2,4)],function(x){shapiro.test(x)$p.value})->pval

# Ajustement de la valeur p pour éviter les erreurs de type 1

# pour en savoir plus sur les ajustements des p-values : cliquer ici

p.adjust(pval,"bonferroni") -> pval# On multiplie les tests alors on ajuste la p-value

# Si une p-value est inférieur à 0.05 (TRUE), c'est qu'un échantillon n'est pas normale, il faut alors penser au test de Levene, sinon faire un test de Bartlett.

any(pval<0.05)

2. Si les groupes sont normaux, faire un de Bartlett

library(stats)

bartlett.test(data$Note,paste(data$Sexe,data$Lycée))

# Si p-value > 0.05, alors les variances sont équivalents

3. Si les groupes ne sont pas normaux, faire un de Levene ou Brown-Forsyth

# Si mes données sont fortement asymétriques (cf. Skewness ou Kurtosis), on va préférer le test de comparaison des médianes de Brown-Forsyth (comparaison par rapport ayx médianes) ou le test de Levene (comparaison par rapport aux moyennes)

# Levene (sur les moyennes)

# Charger le package lawstat

library(lawstat)

# Effectuer le test de Levene qui bascule automatiqueet sur Brown-Forsyth si nécessaire

levene.test(data$Note, paste(data$Sexe,data$Lycée))

# Brown-Forsyth (sur les médianes)

# Charger le package onewaytests

library(onewaytests)

# Effectuer le test de Brown-Forsythe

bf.test(Note ~ paste(Sexe,Lycée), data = data)

Ces deux tests ne sont pas toujours d'accord, bf.test() semblant particulièrement sensible.

4. Si les variances persistent à sembler différentes, faire un contrôle robuste sur les rangs avec le test de Fligner-Killeen

# Fligner-Killeen (sur les rangs)

fligner.test(Note ~ paste(Sexe,Lycée)) # pval>0.05, variances équivalentes

Quand les hypothèses ne sont pas vérifiées, il faut aller creuser vers les anovas non paramétriques telles le test de Scheirer-Ray-Hare (équivalent à deux facteurs du test de Kruskal-Wallis), l'anova par permutation, l'anova sur données trimmées  (t2way()) et l'anova sur médianes (med2way).

Il est très difficile de dire quelle anova sera la plus pertinente (ce n'est pas gravé dans le marbre). Si le test Fligner-Killeen trouve des variances équivalentes, c'est que, au fond, la variance est grosso modo homogène et peut-être faussée par quelques données extrêmes, on peut alors travailler sur des données trimmées (mais il en faut beaucoup puisque  cela revient à éliminer un pourcentage des données hautes et un pourcentages des données basses).

Si on cherche un test conservateur, le test de Scheirer-Ray-Hare peut être une bonne solution.

Si en revanche, cela ne passe pas, il va falloir envisager de travailler sur les médianes ou appliquer une anova par permutation si l'échantillon est petit. 

Trimmer, c'est éliminer les données les plus élevées et les plus basses (10% de chaquecôté  ici)/

Cela permet de passer en  force face aux données  extrêmes avec un test conservateur.

• Attention, t2way() ne marche que pour 2 facteurs

• Il semble données des résultats indépendamment de l'équilibre des données.

# Charger le package

library(WRS2)

data$Sexe <- factor(data$Sexe) # A faire si on veut pas d'erreur

data$Lycée<- factor(data$Lycée) 

t2way(Note~Sexe+Lycée,data=data,tr=0.2) -> model ; print(model)

Le test de Scheirer-Ray-Hare. est basé sur les rangs. Si certaines sources le décrivent comme conservateur, je le vois plutôt à mi-chemin.

• Attention, le test de Scheirer-Ray-Hare ne marche que pour 2 facteurs uniquement (il faut changer de méthodes en présence de 3 facteurs ou plus).

• Ce test est adapté à des données déséquilibrées.

library(rcompanion)

scheirerRayHare(Note~Sexe*Lycée, data = data, type=2)

En sortie :

• Sum Sq illustre la somme des carrés et donc l'influence du facteur sur la variable : si élevé, effet fort.

• H : mesure la différence moyenne entre les groupes. Un facteur qui a un effet a donc un effet H fort.

• p-value : indique la significativité de l'effet du facteur sur la variable.

Ainsi, p-value confirme un effet du facteur, H aide à le quantifier.

Après une telle ANOVA, il est possible de faire un test de Dunn ou un Wilcoxon par paires, puis de retrouver les lettres des groupes avec la fonction cldList() pour Dunn et catego() de {KefiR}.

Source externe : ANOVA non paramétrique à 2 facteurs ou le test de Scheirer-Ray-Hare - DellaData 

L’ANOVA par permutation ne nécessite pas la vérification des hypothèses de normalité et d'homogénéité des variances en ne comparant pas les distributions à des distributions théories mais en échantillonnant pour comparer les distributions entre elles.

Elle s'adapte bien aux petits échantillons non normaux et me donne le sentiment d'être plus puissante que les tests précédents (il faudrait une étude comparative).

Le package {lmPerm} permet de contrôler les effets de permutations

• L'ANOVA par permutation marche avec plus de deux facteurs !

• Ce test est en revanche sensible aux données déséquilibrées.

library(lmPerm)

model <- aovp(Note~Lycée*Sexe, 

              contrasts = list(Lycée=contr.sum,

                             Sexe=contr.sum),

              data=data)

summary(model)

Un test conservateur pour passer en force.

# Charger le package

library(WRS2)

# Penser à convertir les facteurs avec factor()

med2way(Note~Sexe+Lycée,data=data) -> model

D'après Consensus, Baiyu Zhou et W. Wong ont proposé une méthode ANOVA non paramétrique basée sur le bootstrap (NANOVA) pour les données de microarray avec des designs factoriels. Bien que développée pour les données de microarray, cette méthode peut être appliquée à d'autres types de données et permet l'analyse de designs factoriels, y compris ceux avec des designs équilibrés et déséquilibrés (Zhou & Wong, 2011).

Il faudrait pouvoir programmer cela sous R bien qu'il n'existe pas de package utile à ce sujet.

La fonction anova() est parfaite pour comparer différents modèles (source pour en savoir plus).

Un regret, elle ne fonctionne pas sur les ANOVA.s avec effet aléatoire ou les ANOVA.s de type II ou III. Certains reviennent à des ANOVA.s de type I pour comparer les modèles, mais cela manque de pertinence, ou alors il faut tester plusieurs ordres de saisi des facteurs pour vérifier que le déséquilibre n'a pas trop d'importance.

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

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

model3 <- aov(Note~Lycée,data=data)

model4 <- aov(Note~Lycée*Sexe*Options,data=data)

anova(model,model2,model3,model4)

Analysis of Variance Table

Model 1: Note ~ Lycée * Sexe

Model 2: Note ~ Lycée + Sexe

Model 3: Note ~ Lycée

Model 4: Note ~ Lycée * Sexe * Options

  Res.Df    RSS Df Sum of Sq       F    Pr(>F)    

1    232 4126.8                                   

2    235 4708.4 -3   -581.61 17.6849 2.443e-10 ***

3    236 4718.4 -1    -10.00  0.9126    0.3405    

4    224 2455.6 12   2262.78 17.2009 < 2.2e-16 ***

Voici comment interpréter la sortie :

• Res.Df : Il s’agit des degrés de liberté résiduels pour chaque modèle. Les degrés de liberté résiduels sont calculés comme le nombre total d’observations moins le nombre de paramètres dans le modèle.

• RSS : Il s’agit de la somme des carrés résiduels pour chaque modèle. C’est une mesure de la quantité d’erreur restante après l’ajustement du modèle. RSS doit donc être bas.

• Df : Il s’agit de la différence dans les degrés de liberté entre le modèle actuel et le modèle précédent.

• Sum of Sq : Il s’agit de la différence dans la somme des carrés résiduels entre le modèle actuel et le modèle précédent.

• F : Il s’agit de la statistique F pour le test de l'hypothèse que le modèle actuel est significativement meilleur que le modèle précédent.

• Pr(>F) : Il s’agit de la valeur p pour le test d’hypothèse. Une valeur p inférieure à 0,05 indique généralement que le modèle actuel est significativement meilleur que le modèle précédent.

Dans votre cas, le tableau ANOVA indique que le model 4 (Note ~ Lycée * Sexe * Options) est significativement meilleur que le mode l1 (Note ~ Lycée * Sexe) car la valeur p (Pr(>F)) est inférieure à 0,001 (indiqué par ***). Pas de différence en revanche entre les modèle 2 et 3?

On peut aussi extraire la p-value de l'ANOVA de la façon suivante :

unlist(summary(model)[[1]]["Pr(>F)"])[1]

Pa défaut,le niveau de référence est le lycée A et  le sexe F. On peut changer cela en faisant :

data$Sexe <- factor(data$Sexe) # Il faut s'assurer que le sexe soit un facteur

data$Sexe <- relevel(data$Sexe, ref = "M")

Les moyennes marginales (estimées ou prédites) permettent de comparer les effets prédits par le modèle.

La fonction emmeans() de {emmeans} permet de retrouver les moyennes marginales estimées (EMMs) avec leurs intervalles de confiance. Ces derniers pouvant être corrigés par un ajustement de la valeur p. Une source extérieure : Moyennes marginales - DellaData 

library(emmeans)

# Moyenne marginales sans ajustement

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

# Moyennes marginales avec ajustement de la valeur p et donc correction de l'intervalle de confiance

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

# exemple ci-contre

# Les Filles du lycée A ont une moyenne de 11.83 avec une confiance allant de 9.62 à 14.04.

# Regroupement des moyennes par Sexe (on ne va pas étudier toutes les combinaisons de facteurs)

emmeans(model,~Lycée|Sexe, adjust = "bonferroni")

Lorsque l'on a un modèle correctement validé, on peut envisager de faire des tests post-hoc pour identifier quelles catégories ou catégories croisées (ex : Lycée A, Sexe "F", Option 1).

Il existe plusieurs tests post-hoc : Tukey, Student paire à paire avec ajustement de la valeurs p, Wilcoxon sur les rangs paire à paire avec ajustement de la valeurs p, lincon() de {WRS2} pour effectuer des comparaisons multiples après une anova non paramétriques ou encore medpb() de {WRS2} pour effectuer des comparaisons multiples après une anova non paramétriques

Les test 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