ANOVA à 2 facteurs avec R : Tutoriel
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.
Liste des packages et jeu de données
Les packages utilisés dans ce tutoriel sont les suivants :
library(tidyverse)
library(car)
library(Rmisc)
library(multcomp) 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).
N1 <- 9
Tension <-rep(c("L", "M", "H"), each=2*N1)
Wool <- rep(rep(c("A", "B"),each= N1),3)
set.seed(42)
L_A<- rnorm(N1,12,1)
L_B <- rnorm(N1,9,1)
M_A <- rnorm(N1,7,1)
M_B <- rnorm(N1,10,1)
H_A <- rnorm(N1,8.5,1)
H_B <- rnorm(N1,7,1)
Usure <- c(L_A, L_B, M_A, M_B, H_A, H_B)
mydata <- data.frame(Usure, Tension, Wool)
mydata$Tension <- as.factor(mydata$Tension)
mydata$Tension <- factor(mydata$Tension, levels=c("L", "M", "H")) Voici un aperçu de ces données :
head(mydata)
## Usure Tension Wool
## 1 13.37096 L A
## 2 11.43530 L A
## 3 12.36313 L A
## 4 12.63286 L A
## 5 12.40427 L A
## 6 11.89388 L A 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.
ggplot(mydata, aes(x=Tension, y=Usure, colour=Wool, fill=Wool))+
geom_point(position=position_jitterdodge(dodge.width=0.7), size=2) +
geom_boxplot(alpha=0.5, position = position_dodge(width=0.8), fatten=NULL)+
scale_colour_manual(values=c(c("#1E90FF", "#D02090", "#FFFFFF")))+ scale_fill_manual(values=c("#1E90FF", "#D02090", "#FFFFFF"))+
stat_summary(fun.y = mean, geom = "errorbar", aes(ymax = ..y.., ymin = ..y..), width=0.65,size = 1.5, linetype = "solid",position = position_dodge(width=0.7))+
ylab("Usure")+
theme_classic() 
Les traits horizontaux représentent les moyennes.
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.
mydata_avg <- summarySE(mydata,measurevar="Usure",groupvars=c("Tension","Wool"))
mydata_avg
## Tension Wool N Usure sd se ci
## 1 L A 9 12.615076 0.8564690 0.2854897 0.6583404
## 2 L B 9 8.935897 1.4346080 0.4782027 1.1027374
## 3 M A 9 6.893546 1.4235480 0.4745160 1.0942359
## 4 M B 9 9.825706 1.0491770 0.3497257 0.8064689
## 5 H A 9 7.888296 0.9320269 0.3106756 0.7164193
## 6 H B 9 7.338493 0.8753724 0.2917908 0.6728708
# pour éviter le chevauchement sur le graph
p <- position_dodge(0.1)
ggplot(mydata_avg, aes(x=Tension, y=Usure, colour=Wool, group=Wool)) +
geom_errorbar(aes(ymin=Usure-ci, ymax=Usure+ci), width=.1, position=p) +
geom_line(position=p, size=2) +
geom_point(position=p, size=1)+
scale_colour_manual(values=c(c("#1E90FF", "#D02090", "#FFFFFF")))+
scale_fill_manual(values=c("#1E90FF", "#D02090", "#FFFFFF"))+
theme_classic() 
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 .
Remarque:
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”.
mydata$grp <- interaction(mydata$Tension, mydata$Wool, sep="_")
library(boot)
# création de la fonction moyenne pour la fonction boot
moyenne<-function(data,indice)
{
data.star <- data[indice]
moy <- mean(data.star,na.rm=TRUE)
}
#permet de fixer la graine des tirages aléatoires et donc d'obtenir des résultats toujours identiques
set.seed(1234)
# realisation des échantillons bootstrap et estimation des moyennes
bL_A <- boot(mydata$Usure[mydata$grp=="L_A"], statistic= moyenne, R=1000)
# Enfin, l'objet `bL_A` que nous venons de créer est passé en arugment de la fonction `boot.ci`, qui calcule l'intervalle de confiance.
boot.ci(bL_A)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = bL_A)
##
## Intervals :
## Level Normal Basic
## 95% (12.09, 13.15 ) (12.08, 13.13 )
##
## Level Percentile BCa
## 95% (12.10, 13.15 ) (12.13, 13.21 )
## Calculations and Intervals on Original Scale La fonction renvoie quatre types d’intervalle de confiance. On utilise généralement les méthodes “percentile” ou “BCa”.
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 :
mod1 <- lm(Usure~Tension * Wool, contrasts=list(Tension=contr.sum, Wool=contr.sum), data=mydata) 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 !
library(car)
Anova(mod1, type=3)
## Anova Table (Type III tests)
##
## Response: Usure
## Sum Sq Df F value Pr(>F)
## (Intercept) 4292.9 1 3409.8331 < 2.2e-16 ***
## Tension 98.4 2 39.0605 8.529e-11 ***
## Wool 2.5 1 2.0037 0.1634
## Tension:Wool 98.4 2 39.0954 8.416e-11 ***
## Residuals 60.4 48
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 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.
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.
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.
L’absence de corrélation entre les résidus et le facteur étudié peut également se vérifier 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.
plot(mod1,1) 
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.
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.
plot(mod1,2) 
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.
shapiro.test(residuals(mod1))
##
## Shapiro-Wilk normality test
##
## data: residuals(mod1)
## W = 0.98451, p-value = 0.7081 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).
plot(mod1,3) 
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.
bartlett.test(residuals(mod1)~mydata$grp)
##
## Bartlett test of homogeneity of variances
##
## data: residuals(mod1) by mydata$grp
## Bartlett's K-squared = 4.3834, df = 5, p-value = 0.4956
leveneTest(residuals(mod1)~mydata$grp)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 0.4043 0.8434
## 48
fligner.test(residuals(mod1)~mydata$grp)
##
## Fligner-Killeen test of homogeneity of variances
##
## data: residuals(mod1) by mydata$grp
## Fligner-Killeen:med chi-squared = 2.3529, df = 5, p-value = 0.7985 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 :
# permet de séparer la fenêtre graphique en 4 parties (2 lignes, et 2 colonnes)
par(mfrow=c(2,2))
plot(mod1) 
Démarche en cas d'interaction significative
Les hypothèses étant validées, les résultats peuvent être interprétés.
Anova(mod1, type=3)
## Anova Table (Type III tests)
##
## Response: Usure
## Sum Sq Df F value Pr(>F)
## (Intercept) 4292.9 1 3409.8331 < 2.2e-16 ***
## Tension 98.4 2 39.0605 8.529e-11 ***
## Wool 2.5 1 2.0037 0.1634
## Tension:Wool 98.4 2 39.0954 8.416e-11 ***
## Residuals 60.4 48
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 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.
ANOVA à un facteur et comparaisons deux à deux
Pour plus de détail sur cette démarche, vous pouvez consulter mon article
mod_grp <- lm(Usure~grp, data=mydata)
Anova(mod_grp)
## Anova Table (Type II tests)
##
## Response: Usure
## Sum Sq Df F value Pr(>F)
## grp 199.315 5 31.663 4.131e-14 ***
## Residuals 60.431 48
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 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”.
library(multcomp)
mc_tukey <- glht(mod_grp, linfct=mcp(grp="Tukey"))
summary(mc_tukey)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = Usure ~ grp, data = mydata) ## ## Linear Hypotheses: ## Estimate Std. Error t value Pr(>|t|)
## M_A - L_A == 0 -5.7215 0.5289 -10.817 < 0.001 ***
## H_A - L_A == 0 -4.7268 0.5289 -8.936 < 0.001 ***
## L_B - L_A == 0 -3.6792 0.5289 -6.956 < 0.001 ***
## M_B - L_A == 0 -2.7894 0.5289 -5.274 < 0.001 ***
## H_B - L_A == 0 -5.2766 0.5289 -9.976 < 0.001 ***
## H_A - M_A == 0 0.9948 0.5289 1.881 0.42613
## L_B - M_A == 0 2.0424 0.5289 3.861 0.00437 **
## M_B - M_A == 0 2.9322 0.5289 5.544 < 0.001 ***
## H_B - M_A == 0 0.4449 0.5289 0.841 0.95813
## L_B - H_A == 0 1.0476 0.5289 1.981 0.36820
## M_B - H_A == 0 1.9374 0.5289 3.663 0.00774 **
## H_B - H_A == 0 -0.5498 0.5289 -1.039 0.90211
## M_B - L_B == 0 0.8898 0.5289 1.682 0.54982
## H_B - L_B == 0 -1.5974 0.5289 -3.020 0.04397 *
## H_B - M_B == 0 -2.4872 0.5289 -4.702 < 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method) Il est également possible de visualiser ces comparaions multiples sur un graph.
par(mar=c(3,7,3,3))
plot(mc_tukey) 
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.
tuk.cld <- cld(mc_tukey) tuk.cld
## L_A M_A H_A L_B M_B H_B
## "d" "a" "ab" "bc" "c" "a" 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.
letters <- tuk.cld$mcletters$Letters
myletters_df <- data.frame(grp=levels(mydata$grp),letters=letters)
myletters_df
## grp letters
## L_A L_A d
## M_A M_A a
## H_A H_A ab
## L_B L_B bc
## M_B M_B c
## H_B H_B a
ggplot(mydata, aes(x=grp, y=Usure, colour=grp, fill=grp))+
geom_boxplot(outlier.alpha = 0, alpha=0.5)+
geom_jitter(width=0.25)+
stat_summary(fun.y=mean, colour="black", geom="point",
shape=18, size=3) +
theme_classic()+
theme(legend.position="none")+
theme(axis.text.x = element_text(angle=30, hjust=1, vjust=1))+
geom_text(data = myletters_df, aes(label = letters, y = 30 ), colour="black", size=5) 
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éepour garder un risque alpha global de 5%, cela peut empêcher la mise en évidence de différences significatives.
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.
Ajustement du modèle ANOVA à un facteur, sans intercept
mod_grp2 <- lm(Usure ~ grp - 1, data = mydata) On peut visualiser les moyennes avec la fonction summary :
summary(mod_grp2)
##
## Call:
## lm(formula = Usure ~ grp - 1, data = mydata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.59235 -0.64873 -0.06734 0.74193 2.35075
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## grpL_A 12.615 0.374 33.73 <2e-16 ***
## grpM_A 6.894 0.374 18.43 <2e-16 ***
## grpH_A 7.888 0.374 21.09 <2e-16 ***
## grpL_B 8.936 0.374 23.89 <2e-16 ***
## grpM_B 9.826 0.374 26.27 <2e-16 ***
## grpH_B 7.338 0.374 19.62 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.122 on 48 degrees of freedom
## Multiple R-squared: 0.9867, Adjusted R-squared: 0.9851
## F-statistic: 594.7 on 6 and 48 DF, p-value: < 2.2e-16 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 :
Tukey <- contrMat(table(mydata$Tension), "Tukey")
K1 <- cbind(Tukey, matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)))
rownames(K1) <- paste(levels(mydata$Wool)[1], rownames(K1), sep = ":")
K1
## L M H
## A:M - L -1 1 0 0 0 0
## A:H - L -1 0 1 0 0 0
## A:H - M 0 -1 1 0 0 0 Puis, de la même façon, on construit la matrice des contrastes pour la laine de type B :
K2 <- cbind(matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)), Tukey)
rownames(K2) <- paste(levels(mydata$Wool)[2], rownames(K2), sep = ":")
K2
L M H
B:M - L 0 0 0 -1 1 0
B:H - L 0 0 0 -1 0 1
B:H - M 0 0 0 0 -1 1 Enfin, on assemble les deux matrices :
K <- rbind(K1, K2)
colnames(K) <- c(colnames(Tukey), colnames(Tukey))
K
## L M H L M H
## A:M - L -1 1 0 0 0 0
## A:H - L -1 0 1 0 0 0
## A:H - M 0 -1 1 0 0 0
## B:M - L 0 0 0 -1 1 0
## B:H - L 0 0 0 -1 0 1
## B:H - M 0 0 0 0 -1 1 Réalisation des comparaisons souhaitées
Enfin, on obtient les comparaisons souhaitées :
summary(glht(mod_grp2, linfct = K))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = Usure ~ grp - 1, data = mydata)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## A:M - L == 0 -5.7215 0.5289 -10.817 < 1e-04 ***
## A:H - L == 0 -4.7268 0.5289 -8.936 < 1e-04 ***
## A:H - M == 0 0.9948 0.5289 1.881 0.283048
## B:M - L == 0 0.8898 0.5289 1.682 0.390740
## B:H - L == 0 -1.5974 0.5289 -3.020 0.021907 *
## B:H - M == 0 -2.4872 0.5289 -4.702 0.000127 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method) Ci-dessous les commandes pour comparer les moyennes d’usure de laines Aet B pour chaque niveau de tension :
Tukey <- contrMat(table(mydata$Wool), "Tukey")
K3 <- cbind(Tukey,
matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)),
matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)))
rownames(K3) <- paste(levels(mydata$Tension)[1], rownames(K3), sep = ":")
K3
## A B
## L:B - A -1 1 0 0 0 0
K4 <- cbind(matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)),
Tukey,
matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)))
rownames(K4) <- paste(levels(mydata$Tension)[2], rownames(K4), sep = ":")
K4
## A B
## M:B - A 0 0 -1 1 0 0
K5 <- cbind(matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)),
matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)),
Tukey)
rownames(K5) <- paste(levels(mydata$Tension)[3], rownames(K5), sep = ":")
K6 <- rbind(K3, K4, K5)
colnames(K6) <- c(colnames(Tukey), colnames(Tukey), colnames(Tukey))
K6
## A B A B A B
## L:B - A -1 1 0 0 0 0
## M:B - A 0 0 -1 1 0 0
## H:B - A 0 0 0 0 -1 1
summary(glht(mod_grp2, linfct = K6))
##
## Simultaneous Tests for General Linear Hypotheses ## ## Fit: lm(formula = Usure ~ grp - 1, data = mydata) ## ## Linear Hypotheses: ## Estimate Std. Error t value Pr(>|t|)
## L:B - A == 0 -5.7215 0.5289 -10.817 <1e-04 ***
## M:B - A == 0 1.0476 0.5289 1.981 0.15
## H:B - A == 0 -2.4872 0.5289 -4.702 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method) Démarche en cas d'interaction quantitative 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.
Démarche en cas d'interaction non significative
Pour illustrer cette démarche, je vais utiliser une autre simulation de données:
N1 <- 9
Tension <-rep(c("L", "M", "H"), each=2*N1)
Wool <- rep(rep(c("A", "B"),each= N1),3)
set.seed(42)
L_A<- rnorm(N1,12,1)
L_B <- rnorm(N1,10,1)
M_A <- rnorm(N1,8,1)
M_B <- rnorm(N1,6,1)
H_A <- rnorm(N1,10,1)
H_B <- rnorm(N1,7,1)
Usure <- c(L_A, L_B, M_A, M_B, H_A, H_B)
mydata2 <- data.frame(Usure, Tension, Wool)
mydata2$Tension <- as.factor(mydata2$Tension)
mydata2$Tension <- factor(mydata2$Tension, levels=c("L", "M", "H"))
mydata2_avg <- summarySE(mydata2,measurevar="Usure",groupvars=c("Tension","Wool"))
mydata2_avg
## Tension Wool N Usure sd se ci
## 1 L A 9 12.615076 0.8564690 0.2854897 0.6583404
## 2 L B 9 9.935897 1.4346080 0.4782027 1.1027374
## 3 M A 9 7.893546 1.4235480 0.4745160 1.0942359
## 4 M B 9 5.825706 1.0491770 0.3497257 0.8064689
## 5 H A 9 9.388296 0.9320269 0.3106756 0.7164193
## 6 H B 9 7.338493 0.8753724 0.2917908 0.6728708
# pour éviter le chevauchement sur le graph
p <- position_dodge(0.1)
ggplot(mydata2_avg, aes(x=Tension, y=Usure, colour=Wool, group=Wool)) +
geom_errorbar(aes(ymin=Usure-ci, ymax=Usure+ci), width=.1, position=p) +
geom_line(position=p, size=2) +
geom_point(position=p, size=1)+
scale_colour_manual(values=c(c("#1E90FF", "#D02090", "#FFFFFF")))+
scale_fill_manual(values=c("#1E90FF", "#D02090", "#FFFFFF"))+
theme_classic() 
Ici, les profils ne se croisent pas, on s’attend donc à une interaction non significative.
mod2 <- lm(Usure~Tension * Wool, contrasts=list(Tension=contr.sum, Wool=contr.sum), data=mydata2) Anova(mod2, type=3) ## Anova Table (Type III tests) ## ## Response: Usure ## Sum Sq Df F value Pr(>F)
## (Intercept) 4213.0 1 3346.3922 < 2.2e-16 ***
## Tension 181.4 2 72.0620 3.498e-15 ***
## Wool 69.3 1 55.0409 1.681e-09 ***
## Tension:Wool 1.2 2 0.4588 0.6348
## Residuals 60.4 48
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 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.
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“.
mod3 <- lm(Usure~Tension + Wool, data=mydata2) Anova(mod3) ## Anova Table (Type II tests) ## ## Response: Usure ## Sum Sq Df F value Pr(>F)
## Tension 181.449 2 73.656 1.246e-15 ***
## Wool 69.295 1 56.259 9.916e-10 ***
## Residuals 61.586 50
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 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.
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éunit dans une seule matrice , qui est donnée en argument de la fonction glht().
K1 <- glht(mod3, mcp(Wool = "Tukey"))$linfct
K1
## (Intercept) TensionM TensionH WoolB
## B - A 0 0 0 1
## attr(,"type")
## [1] "Tukey"
K2 <- glht(mod3, mcp(Tension = "Tukey"))$linfct
K2
## (Intercept) TensionM TensionH WoolB
## M - L 0 1 0 0
## H - L 0 0 1 0
## H - M 0 -1 1 0
## attr(,"type")
## [1] "Tukey"
K <- rbind(K1,K2) summary(glht(mod3, linfct = K)) ## ## Simultaneous Tests for General Linear Hypotheses ## ## Fit: lm(formula = Usure ~ Tension + Wool, data = mydata2) ## ## Linear Hypotheses: ## Estimate Std. Error t value Pr(>|t|)
## B - A == 0 -2.2656 0.3021 -7.501 < 1e-04 ***
## M - L == 0 -4.4159 0.3699 -11.937 < 1e-04 ***
## H - L == 0 -2.9121 0.3699 -7.872 < 1e-04 ***
## H - M == 0 1.5038 0.3699 4.065 0.000646 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
plot(summary(glht(mod3, linfct = K))) 
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.
Si cet article vous a plu, ou vous a été utile, et si vous le souhaitez, vous pouvez soutenir ce blog en faisant un don sur sa page Tipeee 
Poursuivez votre lecture
- Introduction à l’ANOVA à 2 facteurs
- ANOVA à 2 facteurs : principe
- ANOVA à un facteur : partie 1
- ANOVA à un facteur : partie 2 – la pratique
- ANOVA à un facteur : quand les hypothèses ne sont pas satisfaites
- Tutoriel : comparaison de 2 moyennes avec le logiciel R
- Comparaison de moyennes : indiquez la significativité des différences sur le graph
Très très très instructif et enrichissant cet article merci.
C’est très intéressant, Merci Bien.
Merci Salima.
Très instructif. Merci infiniment.
Je viens de lire votre article, clair et précis. Je réalise justement des Anova à deux facteurs (des Ancova plus précisément) et il y a un cas de figure que je ne retrouve pas dans vos exemples. Peut-on expliquer la présence d’un effet d’interaction significatif lors du test du modèle, et à la fois aucune paire significative lors du POST HOC ? Cela me semble étrange mais je me retrouve face à ce cas de figure plus d’une fois dans mes analyses. Je vous remercie pour les informations.
Bonjour Juliette,
oui je pense que c’est possible. Les tests posthoc ajustent les p-values pour garder un risque alpha global de 5%. Plus il ya de comparaisons plus les p values sont “réhaussées”, et au finall elles peuvent toutes être >0.05.
Bonne continuation.
L’excellente qualité de tous vos didacticiels permettent de réactualiser avec efficacité tous mes anciens cours et facilitent la mise en œuvre sous R des tests statistiques toujours plus complets avec les nombreuses et nouvelles librairies. Bravo et merci pour ce travail.
Une question en lien avec l’anova 2 voies (facteur A à deux modalités a1, a2 et facteur B à deux modalités b1, b2). En absence d’interaction A:B, si le facteur B n’a pas d’effet peut-on regrouper les deux modalités b1 et b2 pour ne considérer que le facteur A ?
Merci de votre aide
Bonjour,
merci pour vos tutoriels très clairs,
j’avais une question sur l’anova à deux facteurs avec interaction : dans le cas où les résidus du modèle ne sont pas normaux, est-il possible de réaliser un test de Kruskal-Walis avec interaction? et de réaliser un test post-hoc par la suite?
merci pour votre réponse,
Bonjour,
le test de Kruskal Wallis est une anova non paramétrique à 1 facteur, on ne peut donc pas tester l’interaction. Si vos résidus ne sont pas normaux vous pouvez essayer une transfo log de la réponse.
Bonne continuation.
Bonjour Claire
Merci pour le sacrifice.
Moi j’ai un soucis, mes données ne sont pas normaux et j’ai fait une transformation logarithmique mais la normalité ne se vérifie pas toujours. Or je veux tester l’interaction, que dois-je faire?
Merci de m’aider
Bonjour,
vous trouverez des pistes ici : https://www.researchgate.net/post/Is_there_a_non-parametric_equivalent_of_a_2-way_ANOVA
Bonne continuation
Bonjour,
Que faire dans le cas de figures où les hypothèses de validité de l’ANOVA à 2 facteurs ne sont pas satisfaites? Quelles sont les solutions envisageables ? Je suis dans l’impasse dans mes analyses … Merci par avance
Bonjour, vous pouvez consulter cet article : https://statistique-et-logiciel-r.com/anova-a-2-facteurs-quand-les-hypotheses-ne-sont-pas-satisfaites/
Bonjour,
Merci pour ce tutoriel.
Je dois réaliser une ANOVA à 3 facteur (type carré latin) avec des données déséquilibrées (manque de données). J’ai essayé de faire l’ANOVA avec les “*” pour prendre en compte les interactions même si je sais qu’avec un carré latin les interactions sont supposées nulles mais on m’a demandé de les tester quand même (enfin de tester une des interactions en particulier) … Cependant, j’ai un message d’erreur en me disant “residual df = 0″. J’ai alors essayé avec les ” ” et donc en supposant nulles les interactions et dans ce cas l’ANOVA fonctionne.
Pourquoi ne puis je pas tester les interactions comme dans votre exemple ?
Merci d’avance
Cordialement
Super intéressant merci !!
Bonjour,
Merci beaucoup pour ce tutoriel ! Il est extrêmement bien fait. J’avais deux questions à propos de l’ANOVA à deux facteurs :
-Est-ce normal à l’étape 6.1, lors des comparaisons multiples avec la fonction glht, d’avoir toutes les erreurs standards (SE) égales entre elles ?
-Peut-on utiliser la fonction lsmeans à la place de glht, afin de mieux contrôler les comparaisons multiples que l’on souhaite réaliser ?
Merci d’avance
Bonjour Claire-Anne,
oui, c’est normal que toutes les SE soit identique, car le calcul est fait de manière globale.
Je ne connais pas bien lsmeans, mais par contre, actuellement j’utilise le package emmeans pour faire les comparaisons multiples, et je le trouve très bien.
Vous trouverez sa vignette ici : https://cran.r-project.org/web/packages/emmeans/vignettes/comparisons.html.
Bonne continuation
Bonjour Claire,
Je vous remercie pour votre blog et vos articles très pédagogiques et facile à comprendre. Je m’inspire souvent de votre blog et je le fais connaître sans aucune retenue!!
Une question par rapport à l’ANOVA en général. Sauf erreur de ma part, vous n’abordez pas la partie “Modélisation” et donc des coefficients estimés pour chaque niveau des facteurs. Je doute sur l’écriture de mon modèle lorsque l’effet de l’interaction est significative . Dois-je mettre les coefficients associés à chacun des facteurs A et B pris indépendamment l’un de l’autre?
Connaissez vous une référence où l’on en parle?
Je vous remercie .
Très bonne continuation.
Céline
Bonjour,
oui, il faut mettre les coefficients associés à chacun des facteurs A et B pris indépendamment l’un de l’autre?
Pour éviter toute erreur, je vous recommande d’employer le package equatiomatic et la fonction extract_eq, comme cela :
extract_eq(mod1, use_coefs=TRUE)Vous trouverez plus d’infos sur l’utilisation de ce package dans cet article : https://delladata.fr/equatiomatic-package/.
Bonne continuation
Merci infiniment pour la mise à jour.
Bonjour, je ne comprends pas la ligne (lorsque les résultats de l’anova à deux facteurs apparaissent) tension:wool.
Merci d’avance si pouvez répondre.
Bonjour
il s’agit de l’interaction tension * wool, R l’écrit avec la syntaxe tension:wool.
Bonne continuation.
Bonjour, du fond de mon coeur j’aimerais vius dire merci pour vrai ce tutoriel qui m’a grandement aidé à comprendre les analyses. Mais svp, j’aimerais savoir comment faire dans le cas d’une analyse Multivariée où j’ai plusieurs réponses comment faire pour avoir les résultats pour toute les réponses. Ou bien faut-il seulement procéder réponse par réponse
[…] Les analyses de variance ANOVA : tutoriel approfondi sur le blog de Claire Della Vedova. […]
Svp, aussi j’aimerais savoir comment utiliser le test SNK (student newman keul) pour effectuer le test post hoc
Bonjour,
peut être SNK.test() du package Agricolae, mais je ne l’ai jamais employé.
Bonne continuation
Bonjour, du fond du coeur j’aimerais vous dire merci pour ce tutoriel qui m’a grandement aidé à comprendre les analyses. Mais svp, j’aimerais savoir comment faire dans le cas d’une analyse Multivariée où j’ai plusieurs réponses..? comment faire pour avoir les résultats pour toute les réponses. Ou bien faut-il seulement procéder réponse par réponse..?
Cordialement
Bonjour,
Regardez du côté des MANOVA, je ne connais pas bien cette approche, mais cela pourrait correspondre à votre besoin.
Sinon, oui, vous pouvez faire réponse par réponse.
Bonne continuation.
Bonjour et Merci pour toute cette mine d’informations très precieuses. Si je dois faire une Anova à 1 ou plusieurs facteurs pour analyser 10 variables par exemple, en “enchainant” les Anovas, j’augmente mon risque d’erreur alpha, je dois donc corriger cette inflation en appliquant une correction multiple. Oui mais dans la pratique comment faire pour des anova à 1 ou 2 facteurs ? Bien à vous.
Bonjour Emmanuelle,
à ma connaissance, il n’y a pas vraiment de consensus sur le sujet. A priori, ce que je ferai, c’est un ajustement par ANOVA. Sinon, vous pouvez aussi faire 3 types de comparaisons : 1) sans correction, 2) avec correction par anova 3) pour toutes les anova. Puis faire un tableau recapitulatif avec les 3 pvalues obtenues.
J’espère que ma réponse vous aide un peu.
Bonne continuation
Bonjour Claire,
J’ai essayé votre tutoriel mais en faisant un test de tukey unilatéral. Je ne comprends pas les résultats que j’obtiens avec :
mc_tukey <- glht(mod_grp, linfct=mcp(grp="Tukey"), alternative ="greater")
summary(mc_tukey)
ou bien
mc_tukey = 0 (avec pval <0.001 ***, et t-value =-6.956) ?
Par avance, je vous remercie pour l'aide que vous pourrez m'apporter.
Bonne journée
Bonjour,
je ne comprends pas votre question, pourriez vous copier-coller l’ensemble de la sortie, s’il vous plait ?
Merci.
Bonjour Claire,
Merci, voici ce que j’obtiens :
> library(multcomp)
> mc_tukey summary(mc_tukey)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lm(formula = Usure ~ grp, data = mydata)
Linear Hypotheses:
Estimate Std. Error t value Pr(= 0 -5.7215 0.5289 -10.817 = 0 -4.7268 0.5289 -8.936 = 0 -3.6792 0.5289 -6.956 = 0 -2.7894 0.5289 -5.274 = 0 -5.2766 0.5289 -9.976 = 0 0.9948 0.5289 1.881 1.0000
L_B – M_A >= 0 2.0424 0.5289 3.861 1.0000
M_B – M_A >= 0 2.9322 0.5289 5.544 1.0000
H_B – M_A >= 0 0.4449 0.5289 0.841 1.0000
L_B – H_A >= 0 1.0476 0.5289 1.981 1.0000
M_B – H_A >= 0 1.9374 0.5289 3.663 1.0000
H_B – H_A >= 0 -0.5498 0.5289 -1.039 0.7513
M_B – L_B >= 0 0.8898 0.5289 1.682 1.0000
H_B – L_B >= 0 -1.5974 0.5289 -3.020 0.0243 *
H_B – M_B >= 0 -2.4872 0.5289 -4.702 <0.001 ***
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported — single-step method)
Ma question est : est-ce que ce sont bien les hypothèses alternatives que l'on voit ici ?
Je suis vraiment vraiment novice en stats, en R…
merci par avance pour votre aide.
Bonjour,
Je pense qu’il s’agit plutôt des hypothèses nulles.
Ma compréhension de la ligne suivante
H_B – L_B >= 0 -1.5974 0.5289 -3.020 0.0243 *
H0 : HB -LB >=0, H1 : HB- LB <0
Ici pvalue <0.05 donc H0 est rejetée donc HB - LB <0 autrement dit HB < LB Vous pouvez conclure que HB est significativement inférieur à LB. J'espère que cela vous aide. Bonne continuation
Bonjour,
Ah d’accord ! Je vous remercie beaucoup.
Cordialement,
Pauline
Bonsoir ,
je n’ai pas compris la difference entre interception qualitative et quantitative . De plus , le lien qui devait l’expliquer echoue sur une page a vendre .
merci par avance
Bonjour,
quand les profils se croisent certains auteurs parlent d’interaction qualitative (les 2 variables intragissent et conduisent à des profils de réponse “opposé”)
Quand les profils ne se croisent pas, mais ne sont pas parallèles certains auteurs parlent d’interaction quantitative (les 2 variables interagissent ensemble, mais dans le même sens).
J’espère que cela vous aide.
Bien cordialement.
Bonjour.
J’essaie d’appliquer cet exercice en utilsant votre cod sauf que je bloque à la 1ere ligne.
N1 n’est pas reconnue par R.
Merci pour l’eclaircissement
Bonjour,
c’est parce que les deux premières lignes étaient collées, il faut lire ;
N1 <- 9 Tension <-rep(c("L", "M", "H"), each=2*N1) Wool <- rep(rep(c("A", "B"),each= N1),3bonjour,
oui je ne l’ai remarqué qu’après et désolé pour mon etourderie. En tous cas merci pour votre réponse.
Pour m’excuser et à toute fin utile, je vous livre un résultat de ma recherche internet afin de connaitre le calcul de p-value du test de tukey (et que je ne trouve pas dans votre code):
D’abord il faut avoir sous la main la table “Studentized Range q Table” (site: https://real-statistics.com/statistics-tables/studentized-range-q-table/)
et pour le calcul de la p-value d’une paire de modalité donnée dans R c’est :
p41 = 1-ptukey(4.334627,4,16) #pour le couple (4;1) pv= 0.0336561<0.05 rejet Ho
p41 : paire des 2 modalités 4 et 1
ptukey = probabilité selon les 3 arguments (la stat du test tukey pour le couple 4 et 1, le DDL de l'inter-groupe; le DDL le l'intra(résidus))
Dans R le ptukey est donné à 95%
Bonjour,
Dans le cadre d’une analyse de la variance à deux facteurs (sans interaction), si les conditions (homoscédasticité et normalité) ne sont pas vérifiées et si les deux facteurs sont significatives (test de permutation pour l’anova à deux facteurs).
Quel test permet de réaliser les comparaisons multiples pour ces deux facteurs ?
Sachant qu’on est pas dans le cas de test de Friedman (avec mesure répétée ou par bloc) pour pouvoir appliquer frdAllPairsConoverTest.
Merci par avance,
Bonjour Eva,
Dans cette situation, j’envisagerai 3 pistes :
1) essayer de comprendre pourquoi les conditions ne sont pas vérifiées : la variable à expliquer est-elle bien numérique continue ? Est-ce qu’un autre type de modélisation ne serait pas plus adaptée (par exemple GLM de type poisson si les réponses à expliquer sont des comptages)
2) employer une transformation BoxCox pour améliorer la normalité et l’homogénéité des résidus, et si cela fonctionne, utiliser la fonction emmeans() pour obtenir les comparaisons
3) réaliser des anova non paramétriques à un facteur (Kruskall Wallis) en fixant la modalité d’un des facteurs, puis les comparaison multiples avec des testx de wilcoxon et en ajustant la -value pour prendre en compte l’inflation du risque alpha , ou des tests de Dunn.
J’espère que cela vous aide.
Merci beaucoup pour votre retour! Effectivement transformation de type BoxCox a résolu le problème.
Bonne journée,
Bonjour,
merci pour cet excellent travail
un eclaircissement svp : le test de leven appliqué à 2 facteurs est conditionné par la necessité d’inclure dans le modele
de regression le terme de l’interaction sinon cela marche pas. Moi j’ai calculé le test de leven de chaque facteur à part et j’ai obtenu les résultats suivants:
1/leveneTest(mod1) #pour le modele complet avec interaction
Levene’s Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 5 0.4043 0.8434
48
2/ leveneTest(residuals(mod1)~mydata$Tension) # pour facteur Tension uniquement
Levene’s Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 2 0.7971 0.4562
51
3/ leveneTest(residuals(mod1)~mydata$Wool) #pour le facteur wool
Levene’s Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 0.1153 0.7356
52
On arrive à la meme conclusion : Acceptation de la constance de la variance.
Mais est ce qu’on est obligé de se limiter à tester seulement l’égalité des variances
au niveau uniquement des croisements des 2 facteurs comme vous l’aviez fait.
Et donc est ce que l’égalité des variances au niveau de l’interaction implique celle presente au niveau de
chaque facteur à part?
Merci pour réponse