Paramétrer son ANOVA à deux facteurs
en langage R
L'essentiel de cette page
Faire une anova à deux facteurs sous R, c'est facile !, à première vue, car il y a de nombreux pièges à contrôler pour éviter de se retrouver à faire des interprétations de modèles faux. C'est exactement comment pour les modèles de régressions linéaires, on peut toujours en faire, mais pas nécessairement juste.
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.
1- Démarche générale
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.
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.
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.
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.
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.
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.
2- Un code à copier coller pour simuler un exemple pour faire de l'ANOVA
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 ?
Cliquer ici pour afficher le code simple à copier-coller pour simuler des valeurs...
Lycée <- rep(c("A","B","C","D"), each=60)
établissement <- rep(1:12,each=20)
Sexe <- rep(c("F","M"),120)
Options <- rep(c("Option1","Option2"),each=2,60)
# Le lycée A a une moyenne de 10, contre 12, 8 et 9.5 pour B, C et D
Note <- c(rnorm(60,10,4),rnorm(60,12,4),rnorm(60,8,3),rnorm(60,9.5,3))
# Dans le lycée A, les filles sont bien meilleures (de 2.5 points)
Note[Lycée=="A"&Sexe=="F"] <- Note[Lycée=="A"&Sexe=="F"]+2
Note[Lycée=="B"&Sexe=="F"] <- Note[Lycée=="B"&Sexe=="F"]+3
Note[Lycée=="B"&Sexe=="M"] <- Note[Lycée=="B"&Sexe=="M"]-1
# Dans le lycée C, ce sont les garçons qui sont meilleurs
Note[Lycée=="C"&Sexe=="F"] <- Note[Lycée=="C"&Sexe=="F"]-1.5
Note[Lycée=="C"&Sexe=="M"] <- Note[Lycée=="C"&Sexe=="M"]+1.5
# Créer un biais lié à l'option
# Note[((Lycée=="C")|(Lycée=="A"))&(Options=="Option2")]<- Note[((Lycée=="C")|(Lycée=="A"))&(Classe=="Option2")]-4
# Note[((Lycée=="C")|(Lycée=="A"))&(Options=="Option1")]<- Note[((Lycée=="C")|(Lycée=="A"))&(Options=="Option1")]+4
Note[Options=="Option2"]<- Note[Options=="Option2"]-3
Note[Options=="Option1"]<- Note[Options=="Option1"]+3
# Créer une variable déséquilibrée Sport
# 30% de sportifs aléatoires
Sport <- sample(c("Sportif","Non-Sportif"),60*4,replace=T,p=c(0.3,0.8))
# En réalité 50% : on ne garde comme sportifs que les biens notés.
# Les données sont donc déséquilibrées (plus de sportifs que de non sportifs et un effet sport sur la note)
Sport[Note<median(Note[Sport=="Sportif"])]<-"Non-Sportif"
# Les notes vont de 0 à 20
Note <- ifelse(Note>20,20,Note)
Note <- ifelse(Note<0,0,Note)
Note <- round(Note,0)
# Compilation
data <- data.frame(Note,Lycée,établissement,Sexe,Options,Sport)
library(lattice)
bwplot(Note~Lycée|Sexe+Options,data=data)
La librairie {lattice} permet une prise en main simlpe de R pour croiser les facteurs.
# Charger la bibliothèque ggplot2
library(ggplot2)
# Supposons que 'data' est votre dataframe et qu'il contient les colonnes 'Note', 'Sexe' et 'Lycée'
# data <- votre_dataframe
# Créer un violon plot
p <- ggplot(data, aes(x = Sexe, y = Note, fill = Sexe))
p <- p + geom_violin(trim = FALSE)# ne pas couper sur la plage de données (ici donne des étendus au-delà de 0 ou 20...
p <- p + facet_grid(Options~ Lycée)
#p <- p + facet_grid(~ Lycée)
# Afficher le graphique
print(p)
Cette représentation permet de voir aussitôt un problème de répartition bimodale qui laisse suppose qu'une variable cachée va empêcher les résidus d'être normaux
3- Vérifier que vos données sont équilibrées
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.
# Vérifier l'équilibrage des données
# Ce tableau permet en croisant les facteurs de voir qu'il y a autant d'individus de chaque catégories (10 filles dans le lycée A, 10 garçons dans le lycée B)....
table(data$Lycée,data$Sexe)
F M
A 10 10
B 10 10
C 10 10
D 10 10
# Identifier les réplications, la structure des données
# On voit ici que chaque lycée est représenté 20 fois et chaque sexe 40
replications(Note~Lycée+Sexe,data=data)
Lycée Sexe
20 40
3.1- Si mes données sont équilibrées :
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.
3.2- Si mes données ne sont pas équilibrées :
Mes données ne sont pas équilibrées lorsque les catégories croisées ne donnent pas la même chose.
On le voit aussi lorsque aov() fait des anova (de moindres carrés de type 1) différents selon l'ordre de saisi des facteurs.
table(data$Sport,data$Sexe)
# Selon l'ordre de saisi des facteurs, on n'obtient pas le même modèle anova
aov(Note~Sexe*Sport,data=data)->model; summary(model)
aov(Note~Sport*Sexe,data=data)->model; summary(model)
F M
Non-Sportif 102 105
Sportif 18 15
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/
Remarque : un modèle ANOVA produit avec aov() est impossible à valider (normalité des résidus et homogénéité des variances) car je n'ai pas trouvé de moyen d'accéder à ses résidus. Certains valident en ligne à partir du modèle de type 1, ce qui peut sembler étonnant.
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 1 : Une autre piste à explorer est de rééquilibrer son ANOVA en tirant au sort parmi ses données. Rien n'empêche de renouveler un grand nombre de fois le tirage de sort pour tester plusieurs alternatives. Voici un petit bout de code expérimental (cliquer sur ce texte).
c(rep("F",50),rep("M",30)) -> sx
diplome <- sample(c("A","B","C"),length(sx),p=c(0.50,0.2,0.3),replace=T)
table(sx,diplome)
Y <- rnorm(length(sx),8,4)
uniq <- c("F.A","F.B","F.C","M.A","M.B","M.C")
moy <- 0
for (i in uniq ) {
moy <- moy+1
Y[interaction(sx,diplome)==i] <- Y[interaction(sx,diplome)==i]+1*rnorm(length(Y[interaction(sx,diplome)==i]),moy,1)
}
Y <- ifelse(Y<0,0,Y)
Y <- ifelse(Y>20,20,Y)
Y <- round(Y,0)
boxplot(Y~paste0(sx,diplome))
summary(aov(Y~sx*diplome))
ind <- c()
for (i in uniq) {
ind <- c(ind,sample(which(interaction(sx,diplome)==i),5))
}
compil <- data.frame(sx,diplome,Y) ; summary(aov(Y~sx*diplome,data=compil))
compil <- na.omit(compil[ind,]) ; summary(aov(Y~sx*diplome,data=compil))
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.
4- Etablir un modèle tenant compte des effets aléatoires et des répétitions de mesures
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.
4.1- J'ai potentiellement des effets aléatoires
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.
4.2- J'ai des répétitions sur mes mesures
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.
5- Contrôler mes hypothèses pour m'assurer que mon modèle n'est pas faux
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 :
Normalité des résidus: les résidus (différences entre valeurs prédites et réelles) doivent être normales.
Homogénéité des variances : Les variances des données de chaque groupe doivent être égales.
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)
5.1- Contrôler la qualité du modèle ANOVA
Le Graphique “Residuals vs. Fitted” affiche les résidus en Y en fonctions des valeurs ajustées (prédites par le modèle) par rapport aux valeurs ajustées pour vérifier la normalité des résidus et l'homogénéité de la variance.
Si les résidus ne sont pas répartis de manière aléatoire autour de zéro et ne forment pas une bande horizontale alors le modèle a un problème.
plot(model, which = 1)
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)
5.2- Contrôler la normalité des résidus
Contrôler la normalité des résidus de mon modèle peut se faire par plusieurs approches différentes et complémentaires :
Le Graphique Q-Q (quantile-quantile) exprime la relation entre les quantiles des résidus standardisés (résidus divisés par l'écart-type des résidus) et les quantiles théoriques que devraient avoir ces résidus standardisés s'ils étaient normaux.
Si les points ne sont pas alignés sur la droite en pointillée, cela indique une distribution non normale des résidus.
plot(model, which = 2)
Voici un exemple de modèles dont les résidus ne sont pas normaux ! Il faut aller creuser ici s'il n'y a pas une sous-variable qui fausse les données.
Exemple : si je cherche à prédire des tailles et que j'ai oublié qu'il y avait la variable sexe (avec des filles et garçons), je vais avoir trop systématiquement des résidus positifs pour les garçons (taille sous-estimée) ou négatifs pour les fille (taille surestimée).
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.
Un simple histogramme peut aussi permettre de vérifier que l'on a une répartition normale des résidus.
hist(model$residuals)
Si le graphique dessine une cloche, c'est que les résidus sont répartis de façon normale.
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.
5.3- Contrôler l'homogénéité de la variance
Le Graphique de disperson mettant en relation les carrés des résidus standardisés en fonction des valeurs théoriques permet de contrôler l'homogénéité de la variance.
Si la courbe rouge dessine une droite horizontale tout va bien.
plot(model, which = 3)
library(dae)
residus <- residuals.aovlist(model)
residus_std <- residus/sd(residus)
residus_std_sqrt <- sqrt(residus_std^2)
plot(fitted(model),residus_std_sqrt)
lines(smooth.spline(fitted(model),residus_std_sqrt,df=10,cv=F),col="red",lwd=2)
Les résidus standardisés en fonction du niveau des facteurs (combinés) permet aussi de contrôler l'homogénéité de la variance.
Sur le graphique ci-dessous, la courbe rouge doit dessiné une droite autour de 0 et les résidus standardisés doivent se répartir de la même façon pour chaque groupe.
plot(model, which = 4)
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.
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
6- Quand les hypothèses ne sont pas vérifiées, penser ANOVA non paramétrique ou ANOVA robuste
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.
Cette comparison en simulant des données en faisant varier l'écart-type confirme que le test de Scheirer-Ray-Hare est l'ANOVA la plus conservatrice et que les ANOVA classiques et par mermutation sont plus puissantes. Si on analyse les erreurs cumulées (Type 1 et Type2), l'ANOVA par mermutation et le test de Scheirer-Ray-Hare se substituent bien à l'ANOVA classique.
6.1- Anova sur données trimmées
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)
6.2- Test de Scheirer-Ray-Hare
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}.
6.3- L'anova par permutation.
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)
6.4- Anova sur médianes
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
6.5- ANOVA non paramétrique basée sur le bootstrap (NANOVA)
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.
7- Valider le meilleur modèle
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]
8- Identifier les moyennes marginales, caractériser les effets
Pour caractériser les effets, une fonction est assez sympathique : confint()
Si j'applique confint(model) cela me renvoie ce résultat.
model<- aov(Note~Lycée+Sexe,data=data)
confint(model)
Cela donne l'intervalle de confiance qui traduit la différence entre chaque valeur et une valeur de référence.
Ainsi :
intercept donne la valeur de base indépendamment du lycée et du sexe (donc la note des F du Lycée A (référence par défaut)) : ici cette note serait de 10.74 à 13.22.
LycéeB, traduit l'écart de note entre les individus du LycéeA et du Lycée B. On voit que globalement, le Lycée B semble meilleur et le LycéeC certainement plus mauvais.
SexeM, traduit l'écart des garçons par rapport aux Filles.
2.5 % 97.5 %
(Intercept) 10.7434675 13.223199190
LycéeB -0.1016534 3.034986707
LycéeC -5.0516534 -1.915013293
LycéeD -3.7516534 -0.615013293
SexeM -2.2089697 0.008969735
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")
Lycée Sexe emmean SE df lower.CL upper.CL
A F 11.83 0.8 232 9.62 14.04
B F 15.17 0.8 232 12.96 17.38
C F 6.30 0.8 232 4.09 8.51
D F 10.17 0.8 232 7.96 12.38
A M 9.87 0.8 232 7.66 12.08
B M 10.43 0.8 232 8.22 12.64
C M 9.07 0.8 232 6.86 11.28
D M 9.17 0.8 232 6.96 11.38
On peut retrouver les moyennes marginales en faisant aussi une prédiction ou en récupérant les valeurs fitted (cf. plus haut)
by(predict(model,data),data[,c(2,4)],mean) -> moyennes # prédites sur les colonnes 2 (Lycée) et 4 (Sexe)
# Mise en tableau
ftable(moyennes)
Sexe F M
Lycée
A 11.833333 9.866667
B 15.166667 10.433333
C 6.300000 9.066667
D 10.166667 9.166667
model <- aov(Note~Lycée+Sexe,data=data)
resultat <- emmeans(model,~Lycée+Sexe, adjust = "bonferroni")
matrice <- xtabs(emmean ~ Lycée + Sexe, data = resultat)
barplot(matrice,beside=T,col=rainbow(4))
cf. le code partie suivante permet d'avoir ce graphique complet avec les lettres des groupes déduit d'un test post-hoc
9- Test post-hoc
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.
Ces graphiques sont à titre d'exemple.
Cette modélisation a été faite à partir de la simulation de 4 échantillons normaux de 18 individus chacun dont seul le 4ème avait une moyenne différente.
paramétrique
Très conservateur. Comparaison de tous les groupes deux à deux- faible nombre de groupes
scheffe.test()
paramétrique
conservateur 'mais moins que Scheffé (risque d’erreur de type I plus élevé)), comparaison des moyennes paires par paires en maintenant le risque global alpha d'erreurs de type 1 - grand nombre de groupes
HSD.test() de {agricolae}
paramétrique
Puissant, comparaison des moyennes paires par paires en les triant pour réduire au maximum le nombre de comparaison.
SNK.test() de {agricolae}
Test t par paires
(ex : test de Bonfererroni)
paramétrique
plus ou moins conservateur selon l'ajustement de la valeur p (Bonferonni est conservateur, erreur de type 2 élevé)
pairwise.t.test(p.adjust.method="bonferroni")
Wilcoxon par paires
non-paramétrique
robuste, comparaison sur les rangs, plus ou moins conservateur selon l'ajustement de la valeur p
wilcox.t.test()
Médianes par paires
non-paramétrique
robuste, comparaison sur les médiane avec bootstrap percentile (estimation de 2,5% percentile et 97,5% percentile)
medpb() de {WRS2}
Moyennes rognées
non-paramétrique
robuste, comparaison sur les moyennes avec données extrêmes rognées (% global tr)
lincon() de {WRS2}
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"))
resultat <- emmeans(model,~Lycée+Sexe, adjust = "bonferroni")
pwpp(resultat, method = "pairwise")
Cette figure permet de comparer rapidement les moyennes (en y) et de voir si les différences (représentées par un trait) entre deux moyennes sont significatives.
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}.