ANOVA à 2 facteurs avec R : Tutoriel

tuto anova avec r

Ce tutoriel dédié à la réalisation d’ANOVA à deux facteurs avec le logiciel R, fait suite à deux premiers articles consacrés à cette approche statistique, l’un d’introduction, l’autre détaillant son fonctionnement.

En lisant ce tutoriel vous apprendrez :

  • comment explorer vos données,
  • comment réaliser une ANOVA à deux facteurs,
  • comment vérifier les hypothèses de validité,
  •  comment interpréter les résultats,
  • que faire lorsque l’interaction est significative, et lorsqu’elle ne l’est pas.

Table des matières

 

0. Liste des packages et jeu de données

Les packages utilisés dans ce tutoriel sont les suivants :

 

Pour illustrer ce tutoriel consacré à l’ANOVA à deux facteurs, je vais utiliser des données simulées inspirées du jeu de données warprbreaks (présent dans le package “dataset”), chargé par défaut. Nous allons imaginer que ces données sont issues d’un plan d’expérience visant à observer l’usure de la laine (dans une unité de votre choix) au sein de métiers à tisser, en fonction de la tension du fil (Low / Medium / High) et du type de laine (Type A / Type B).

 

Voici un aperçu de ces données :

 

1. Exploration visuelle

En premier lieu, il est toujours utile de représenter les données pour se faire une première idée avant de réaliser l’ANOVA. Pour cette étape, j’ai l’habitude d’utiliser des boxplots, en ajoutant par-dessus les données observées. Cela permet de se rendre compte du nombre de données par modalité, ainsi que de leur distribution (plutôt normal, plutôt asymétrique ?), et de la présence éventuelle d’outliers (valeurs extrêmes). Ces deux derniers points pouvant biaiser les résultats de l’analyse.

 

tuto anova avec r

Les traits horizontaux représentent les moyennes.

 

2.Calculer les moyennes et leurs intervalles de confiance

Il peut également être intéressant, avant de réaliser l’ANOVA,  de calculer les moyennes et leur intervalle de confiance. Une façon rapide de le faire est d’employer la fonction “summarySE” du package “Rmisc”. Les estimations des intervalles de confiances sont basées sur une distribution t. Elles sont donc biaisées si les données ne suivent pas une loi Normale. Je vous recommande de les considérer uniquement de façon approximative, elles ne servent qu’à se faire une première idée. L’intérêt d’utiliser la fonction summarySE c’est qu’on peut ensuite utiliser la sortie pour représenter les moyennes et leur intervalle de confiance sur un graph.

 

 

 

tutoriel anova à deux facteurs avec R

 

 

Ici, les profils se croisent, on s’attend donc à ce que l’interaction soit de type qualitative et significative. Pour plus d’informations sur les interactions, consultez cet article .

Pour calculer des intervalles de confiance robustes, par bootstrap, on peut aussi utiliser les commandes décrites dans l’article “Analyses statistiques descriptives de données numériques – partie 2″ Il faut le faire pour chaque groupe. Ici un exemple avec les données du croisement des modalités Tension=Low et Wool=A.

Dans un premier temps, on créé une variable “grp” correspondant au croisement des modalités des deux facteurs, grâce à la fonction “interaction”.

 

La fonction renvoie quatre types d’intervalle de confiance. On utilise généralement les méthodes “percentile” ou “BCa”.

 

3. Réalisation de l’ANOVA à deux facteurs

En première approche, on ajuste toujours un modèle ANOVA à deux facteurs avec un terme d’interaction. Il existe à ce niveau là, une petite difficulté liée au fait que lorsque les effectifs ne sont pas égaux dans chaque groupe (croisement des modalités), la part de variance de l’interaction peut se calculer de plusieurs façons. On parle de carrés de type II ou de type III.

Lorsqu’on ajuste le modèle complet, c’est à dire avec l’interaction, on utilise généralement des carrés de type III. Pour cela, il est nécessaire de changer les contrastes des deux facteurs du format par défaut de type “contr.treatment”, vers le format “contr.sum”. C’est ce qui permet au logiciel de calculer correctement ces carrés de type III.

Pour plus d’infos sur les contrastes, vous pouvez consulter de R book pages 368 à 386. Cette modifications des contrastes peut se faire pendant l’ajustement du modèle, ou en amont de celui-ci. Je préfère le faire au moment de l’ajustement.

Remarque : lorsque les effectifs sont équilibrés, les résultats des carrés de type III et de type II sont strictement identiques. ici, il ne serait donc pas nécessaire de changer les contrastes. Je vous conseille, néaNmoins, de le faire systématiquement pour en prendre l’habitude.

Comme expliqué dans l’article sur l’ANOVA à un facteur, l’ajustement d’une ANOVA peut se faire avec les fonctions lm ou aov, mais j’utilise sytématiquement “lm”, par habitude.

L’interaction est incluse dans le modèle en employant le signe * entre les deux facteurs :

 

4. Visualisation des résultats

Le package “car” dispose d’une fonction Anova (avec un A majuscule) qui permet d’obtenir les résultats avec les carrés de type III.

Attention, il faut utiliser la fonction Anova et pas anova, c’est très important, car les deux fonctions ne fournissent pas les mêmes résultats !

 

Avant de passer à l’interprétation des résultats, il est nécessaire de vérifier que les hypothèses de validité de l’ANOVA à deux facteurs sont satisfaites, car si cela n’est oas le cas, les résultats ne sont pas valides.

 

5 Vérification des hypothèses de validité

Comme évoqué dans le premier article décortiquant le principe de l’ANOVA à deux facteurs , les résultats de cette méthode sont valides (on peut avoir confiance dans les résultats), si ces trois hypothèses sont vérifiées :

  • Les résidus sont indépendants,
  • Les résidus suivent une loi normale de moyenne 0,
  • Les résidus relatifs aux différentes modalités sont homogènes (ils ont globalement la même dispersion), autrement dit leur variance est constante.

5.1 Indépendance des résidus

L’indépendance des résidus signifie que les résidus ne doivent pas être corrélés entre eux. Par exemple, il ne doit pas avoir de lien entre un résidu et celui de la donnée suivante, ou précédente. On voit cela, lorsque des données sont répétées sur des sujets identiques. On parle alors d’autocorrélation des résidus. De la même façon, les résidus ne doivent pas être corrélés au facteur étudié.

L’absence d’autocorrélation se valide par l’étude du plan expérimental : pour réaliser une ANOVA à deux facteurs, il ne doit pas avoir de données répétées. Si c’est le cas, il faut utiliser un autre type de modèle, comme un modèle linéaire à effet mixte.

Il est aussi possible de faire un test statistique, le test de Dubin-Watson. L’hypothèse nulle de ce test spécifie un coefficient d’autocorrélation = 0, alors que l’hypothèse alternative spécifie un coefficient d’autocorrélation différent de 0 . On conclut donc à l’absence d’autocorrélation lorsque la pvalue du test est supérieure à 0.05.

 

Ici la p-value est < 0.05, l’hypothèse H0 n’est donc pas rejetée, et on conclut à l’absence d’auto-corrélation.

L’absence de corrélation entre les résidus et le facteur étudié se vérifie généralement de façon visuelle lors du diagnostic de régression, par un plot des résidus vs fitted values, ou encore des résidus vs les modalités du facteur.

 

tutoriel anova 2 avec R

Ici, la valeur des résidus ne semble pas dépendre du groupe (croisement des modalités) puisqu’ils sont tous globalement centrés sur 0.

 

5.2 Normalité des résidus

Pour vérifier cette hypothèse, on utilise généralement un QQplot et/ou un test de normalité comme le test de Shapiro-Wilk.

 

tutoriel anova 2 avec R

Ici, les points sont bien répartis le long de la ligne, cela signifie que les résidus sont distribués selon une loi normale. Le fait que les points soient centrés sur 0 (sur l’axe des y), montre que leur moyenne est égale à 0.

L’hypothèse nulle du test de normalité de Shapiro Wilk spécifie que les résidus suivent une loi normale, alors que son hypothèse alternative spécifie qu’ils suivent une autre distribution quelconque. Pour accepter la normalité des résidus, il est donc nécessaire d’obtenir une p-value > 0.05.

 

 5.3 Homogénéité des variances

L’hypothèse d’homogénéité des variances, c’est-à-dire l’hypothèse que les résidus ont une variance constante, peut s’évaluer graphiquement et/ou à l’aide d’un test statistique.

La méthode graphique consiste à représenter les résidus standardisés en fonction des valeurs prédites (les moyennes des différents groupes).

 

anova 2 ways with R

On voit ici que les dispersions des résidus (leurs écartements verticaux) relatives à chaque groupe (croisement des modalités des 2 facteurs) sont globalement identiques, l’hypothèse d’homogénéité des résidus est acceptée.

On peut également utiliser le test de Bartlett, le test de Levene, ou encore le test de Fligner-Killeen. Leurs hypothèses nulles spécifient que les variances des différents groupes sont globalement identiques. A l’inverse, leurs hypothèses alternatives spécifient qu’au moins 2 variances (les variances de 2 modalités) sont différentes.

Pour accepter l’hypothèse d’homogénéité des résidus, il est donc nécessaire d’obtenir une p-value > 0.05.

 

Ici, la p-value est largement supérieure à 0.05, l’hypothèse d’homogénéité des résidus est donc acceptée.

NB : Pour visualiser tous les plots du diagnostic de régression en une fois, il est possible d’utiliser les commande suivantes :

 

anova 2 ways with r

6 Démarche en cas d’interaction significative :

Les hypothèses étant validées, les résultats peuvent être interprétés.

 

Ici, sans surprise, compte tenu du croisement des profils de l’usure en fonction de la tension, pour les laines de type A et B, l’interaction est qualitative et significative.

Dans ce cas, comme expliqué dans cet article , il n’est pas possible d’interpréter les effets propres des facteurs “Tension” et “Wool”.

Dans ce cas de figure, deux solutions sont envisageables. La première consiste à faire une ANOVA à un facteur sur la variable “grp” créée précédemment (croisement des modalités des 2 facteurs). Puis, si l’effet est significatif, des comparaisons multiples peuvent être réalisées pour mettre en évidence les moyennes significativement différentes deux à deux.

La seconde solution consiste à réaliser les comparaisons des moyennes relatives aux modalités d’un facteur, séparément pour chacune des modalités de l’autre facteur. Par exemple, comparer les moyennes des usures des tensions L, M et H pour les laines de type A d’une part, et de type B d’autre part.

Dans l’esprit, c’est un peu comme si on faisait une ANOVA à un facteur (qui serait la Tension) et ses comparaisons multiples subséquentes (pour chaque modalité de laine (A ou B)).

En pratique, cette approche nécessite de construire une matrice de contrastes afin de définir les comparaisons souhaitées.

6.1 ANOVA à un facteur et comparaions deux à deux

Pour plus de détail sur cette démarche, vous pouvez consulter mon article ANOVA à un facteur : partie 2 – la pratique .

 

L’ANOVA à un facteur montre un effet significatif de la variable “grp” (croisement des modalités des facteurs Tension et Wool). Les moyennes sont ensuite comparées deux à deux selon l’approche de Tukey, en employant la fonction “glht” du package “multcomp”.

 

 

Il est également possible de visualiser ces comparaions multiples sur un graph.

 

comparaisons multiples avec r

 

Le package multcomp, contient également une fonction cld()  qui permet, dans le cadre du test de Tukey, d’indiquer par des lettres la significativité des comparaisons. Lorsque deux modalités partagent une même lettre, cela signifie que leurs différences ne sont pas significativement différentes. A l’inverse, lorsque deux modalités ne partagent pas de lettres en commun, alors cela signifie que leurs moyennes sont significativement différentes.

 

 

Comme expliqué dans un de mes précédents articles, on peut alors utiliser ces lettres pour les ajouter sur un graph réalisé avec ggplot2.

 

significativité

Cette méthode peut aboutir à réaliser un grand nombre de comparaisons, dont certaines ne sont pas intéressantes. Comme les p-values sont ajustées (<https://wp.me/p93iR1-qf>) pour garder un risque alpha global de 5%, cela peut empêcher la mise en évidence de différences significatives.

6.2 Comparaisons à l’intérieur d’une modalité

Cette approche est plus technique. Elle consiste dans un premier temps à ajuster un modèle ANOVA à un facteur sur la variable “grp” (croisement des modalités), en omettant l’intercept, afin que les paramètres du modèle correspondent aux moyennes des différentes conditions. Dans unsecond temps, il s’agit de construite une matrice de contrastes, correspondant aux comparaisons de moyennes souhaitées. Enfin, le modèle et la matrice sont fournis en argument de la fonction “glht” du package “multcomp”, pour obtenir ces comparaisons.

Dans cette approche, on limite donc les comparaisons à celles qui nous intéressent. Dans l’exemple ci-dessous, on va comparer deux à deux les moyennes de l’usure pour les tensions L, M et H, pour la laine de type Ad’un coté, puis la laine B de l’autre.

 

 6.2.1 Ajustement du modèle ANOVA à un facteur, sans intercept

 

On peut visualiser les moyennes avec la fonction summary :

 

6.2.2 Construction de la matrice de contratses

On commence par construire la matrice des contrastes permettant les comparaisons deux à deux pour la laine de type A :

 

Puis, de la même façon, on construit la matrice des contrastes pour la line de type B :

 

Enfin, on assemble les deux matrices :

 

6.2.3 réalisation des comparaisons souhaitées

Enfin, on obtient les comparaisons souhaitées :

 

Ci-dessous les commandes pour comparer les moyennes d’usure de laines A
et B pour chaque niveau de tension :

 

7. Démarche en cas d’interaction quantitive significative

Dans cette situation, les effets propres des facteurs sont généralement interprétés ( à partir du modèle contenant l’interaction). En fonction de la p-value (inférieure ou supérieure au seuil de significativité choisi) on conclura à la présence d’un effet significatif, ou à l’absence de mise en évidence d’un effet significatif.

Si au moins un des effets est significatif, on réalisera les comparaisons multiples correspondantes, comme décrit dans le paragraphe 6.

 

 

8. Démarche en cas d’interaction non significative

Pour illustrer cette démarche, je vais utiliser une autre simulation de données:

 

interaction

 

Ici, les profils ne se croisent pas, on s’attend donc à une interaction non significative.

 

Ici, l’interaction n’est pas significative. Avant d’interpréter les résultats, on va ajuster à nouveau le modèle de l’ANOVA à deux facteurs, mais sans le terme d’interaction, puisque celle-ci n’est pas significative.

8.1 Ajustement du modèle sans le terme d’interaction

Pour réaliser une ANOVA à deux facteurs, sans terme d’interaction, il suffit de remplacer, dans la formule du modèle,  le signe ” * ” par le signe ” + “. Par ailleurs, lorsque le modèle ne contient pas de terme d’interaction, on utilise les carrés de type II. Pour cela, il suffit simplement d’utiliser les contrastes par défaut qui sont de type “contr.treatment“.

 

Ici, les effets des deux facteurs sont significatifs. Concernant le facteur Tension, cela signifie qu‘au moins deux moyennes d’usure sont différentes entre les tensions L, M et H. Concernant le facteur laine, la même conclusion serait tirée s’il existait plus de deux modalités. Dans notre cas de figure, cela signifie alors que la moyenne d’usure de la laine A est significativement supérieure à celle de la laine B.

8.2 Comparaisons multiples

Là encore, le package multcomp, permet de réaliser toutes les comparaisons en une seule fois, et donc d’ajuster les p-values de façon parfaitement adéquate. Pour cela, on réalise deux matrices de contrastes (K1 et K2), une pour chaque facteur, afin de définir les comparaisons souhaitées. Puis on les réunie dans une seule matrice , qui est donnée en argument de la fonction glht.

 

comparaisons anova 2 r

Comme attendu, les résultats nous montrent que la moyenne d’usure de la laine A et supérieure à celle de la laine B. Et que les moyennes d’usure des tensions sont toutes significativement différentes deux.

Conclusion

J’espère que ce tutoriel permettra au plus grand nombre de réaliser sereinement vos ANOVA à deux facteurs avec R. Si quelque chose vous freine ou vous pose problème lorsque vous faites ce type d’analyse, indiquez le moi en commentaire, et j’essaierais de vous apporter une réponse.

Dans un prochain article, j’aborderai les solutions envisageables (transformation, ANOVA basée sur les tests de permutation, etc..) lorsque les hypothèses de validité de l’ANOVA à 2 facteurs ne sont pas satisfaites.

Et si cet tutoriel vous a plu : partagez-le !

 

 

Retrouvez ici d’autres articles en lien avec celui que vous venez de lire :

 

Partager l'article
  •  
  •  
  •  
  •  
  •  
    7
    Partages
  • 7
  •  
  •  
  •  
  •  

3 commentaires

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *