Comment effectuer une ANOVA bidirectionnelle dans R

Une ANOVA bidirectionnelle (« analyse de variance ») est utilisée pour déterminer s’il existe ou non une différence statistiquement significative entre les moyennes de trois groupes indépendants ou plus qui ont été répartis sur deux facteurs.

Ce tutoriel explique comment effectuer une ANOVA bidirectionnelle dans R.

Exemple : ANOVA bidirectionnelle dans R

Supposons que nous souhaitions déterminer si l’intensité de l’exercice et le sexe ont un impact sur la perte de poids. Dans ce cas, les deux facteurs que nous étudions sont l’exercice et le sexe , et la variable de réponse est la perte de poids, mesurée en livres.

Nous pouvons effectuer une ANOVA bidirectionnelle pour déterminer si l’exercice et le sexe ont un impact sur la perte de poids et pour déterminer s’il existe une interaction entre l’exercice et le sexe sur la perte de poids.

Nous recrutons 30 hommes et 30 femmes pour participer à une expérience dans laquelle nous en assignons au hasard 10 de chacun pour suivre un programme sans exercice, d’exercice léger ou d’exercice intense pendant un mois.

Le code suivant crée le bloc de données avec lequel nous allons travailler :

#make this example reproducible
set.seed(10)

#create data frame
data <- data.frame(gender = rep(c("Male", "Female"), each = 30),
                   exercise = rep(c("None", "Light", "Intense"), each = 10, times = 2),
                   weight_loss = c(runif(10, -3, 3), runif(10, 0, 5), runif(10, 5, 9),
                                   runif(10, -4, 2), runif(10, 0, 3), runif(10, 3, 8)))

#view first six rows of data frame
head(data)

#  gender exercise weight_loss
#1   Male     None  0.04486922
#2   Male     None -1.15938896
#3   Male     None -0.43855400
#4   Male     None  1.15861249
#5   Male     None -2.48918419
#6   Male     None -1.64738030

#see how many participants are in each group
table(data$gender, data$exercise)

#         Intense Light None
#  Female      10    10   10
#  Male        10    10   10

Explorer les données

Avant même d’ajuster le modèle ANOVA bidirectionnel, nous pouvons mieux comprendre les données en trouvant la moyenne et l’écart type de la perte de poids pour chacun des six groupes de traitement à l’aide du package dplyr :

#load dplyr package
library(dplyr)

#find mean and standard deviation of weight loss for each treatment group
data %>%
  group_by(gender, exercise) %>%
  summarise(mean = mean(weight_loss),
            sd = sd(weight_loss))

# A tibble: 6 x 4
# Groups:   gender [2]
#  gender exercise   mean    sd
#          
#1 Female Intense   5.31  1.02 
#2 Female Light     0.920 0.835
#3 Female None     -0.501 1.77 
#4 Male   Intense   7.37  0.928
#5 Male   Light     2.13  1.22 
#6 Male   None     -0.698 1.12 

Nous pouvons également créer un boxplot pour chacun des six groupes de traitement afin de visualiser la répartition de la perte de poids pour chaque groupe :

#set margins so that axis labels on boxplot don't get cut off
par(mar=c(8, 4.1, 4.1, 2.1))

#create boxplots
boxplot(weight_loss ~ gender:exercise,
data = data,
main = "Weight Loss Distribution by Group",
xlab = "Group",
ylab = "Weight Loss",
col = "steelblue",
border = "black", 
las = 2 #make x-axis labels perpendicular
)

Nous pouvons immédiatement constater que les deux groupes qui ont participé à des exercices intenses semblent avoir des valeurs de perte de poids plus élevées. Nous pouvons également constater que les hommes ont tendance à avoir des valeurs de perte de poids plus élevées que les femmes dans les groupes d’exercices intenses et légers .

Ensuite, nous adapterons le modèle ANOVA bidirectionnel à nos données pour voir si ces différences visuelles sont réellement statistiquement significatives.

Ajustement du modèle ANOVA bidirectionnel

La syntaxe générale pour ajuster un modèle ANOVA bidirectionnel dans R est la suivante :

aov (variable de réponse ~ prédicteur_variable1 * prédicteur_variable2, données = ensemble de données)

Notez que le * entre les deux variables prédictives indique que nous souhaitons également tester un effet d’interaction entre les deux variables prédictives.

Dans notre exemple, nous pouvons utiliser le code suivant pour ajuster le modèle ANOVA bidirectionnel, en utilisant poids_perte comme variable de réponse et sexe et exercice comme deux variables prédictives.

Nous pouvons ensuite utiliser la fonction summary() pour afficher le résultat de notre modèle :

#fit the two-way ANOVA model
model <- aov(weight_loss ~ gender * exercise, data = data)

#view the model output
summary(model)

#                Df Sum Sq Mean Sq F value Pr(>F)    
#gender           1   15.8   15.80  11.197 0.0015 ** 
#exercise         2  505.6  252.78 179.087 <2e-16 ***
#gender:exercise  2   13.0    6.51   4.615 0.0141 *  
#Residuals       54   76.2    1.41                   
#---
#Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

À partir des résultats du modèle, nous pouvons voir que le sexe , l’exercice et l’interaction entre les deux variables sont tous statistiquement significatifs au niveau de signification de 0,05.

Vérification des hypothèses du modèle

Avant d’aller plus loin, nous devons vérifier que les hypothèses de notre modèle sont remplies afin que nos résultats du modèle soient fiables. En particulier, une ANOVA bidirectionnelle suppose :

1. Indépendance – les observations de chaque groupe doivent être indépendantes les unes des autres. Puisque nous avons utilisé un plan randomisé , cette hypothèse devrait être satisfaite, nous n’avons donc pas à trop nous en soucier.

2. Normalité – la variable dépendante doit avoir une distribution approximativement normale pour chaque combinaison des groupes des deux facteurs.

Une façon de vérifier cette hypothèse consiste à créer un histogramme des résidus du modèle. Si les résidus sont à peu près normalement distribués, cette hypothèse devrait être satisfaite.

#define model residuals
resid <- model$residuals

#create histogram of residuals
hist(resid, main = "Histogram of Residuals", xlab = "Residuals", col = "steelblue")

Les résidus sont à peu près normalement distribués, nous pouvons donc supposer que l’hypothèse de normalité est remplie.

3. Variance égale – les variances pour chaque groupe sont égales ou approximativement égales.

Une façon de vérifier cette hypothèse consiste à effectuer un test de Levene pour l’égalité des variances en utilisant le package car :

#load car package
library(car)

#conduct Levene's Test for equality of variances
leveneTest(weight_loss ~ gender * exercise, data = data)

#Levene's Test for Homogeneity of Variance (center = median)
#      Df F value Pr(>F)
#group  5  1.8547 0.1177
#      54  

Étant donné que la valeur p du test est supérieure à notre niveau de signification de 0,05, nous pouvons supposer que notre hypothèse d’égalité des variances entre les groupes est remplie.

Analyser les différences de traitement

Une fois que nous avons vérifié que les hypothèses du modèle sont satisfaites, nous pouvons alors effectuer un test post hoc pour déterminer exactement quels groupes de traitement diffèrent les uns des autres.

Pour notre test post hoc, nous utiliserons la fonction TukeyHSD() pour effectuer le test de Tukey pour des comparaisons multiples :

#perform Tukey's Test for multiple comparisons
TukeyHSD(model, conf.level=.95) 

#  Tukey multiple comparisons of means
#    95% family-wise confidence level
#
#Fit: aov(formula = weight_loss ~ gender * exercise, data = data)
#
#$gender
#                diff       lwr      upr     p adj
#Male-Female 1.026456 0.4114451 1.641467 0.0014967
#
#$exercise
#                   diff       lwr       upr   p adj
#Light-Intense -4.813064 -5.718493 -3.907635 0.0e+00
#None-Intense  -6.938966 -7.844395 -6.033537 0.0e+00
#None-Light    -2.125902 -3.031331 -1.220473 1.8e-06
#
#$`gender:exercise`
#                                  diff        lwr         upr     p adj
#Male:Intense-Female:Intense  2.0628297  0.4930588  3.63260067 0.0036746
#Female:Light-Female:Intense -4.3883563 -5.9581272 -2.81858535 0.0000000
#Male:Light-Female:Intense   -3.1749419 -4.7447128 -1.60517092 0.0000027
#Female:None-Female:Intense  -5.8091131 -7.3788841 -4.23934219 0.0000000
#Male:None-Female:Intense    -6.0059891 -7.5757600 -4.43621813 0.0000000
#Female:Light-Male:Intense   -6.4511860 -8.0209570 -4.88141508 0.0000000
#Male:Light-Male:Intense     -5.2377716 -6.8075425 -3.66800066 0.0000000
#Female:None-Male:Intense    -7.8719429 -9.4417138 -6.30217192 0.0000000
#Male:None-Male:Intense      -8.0688188 -9.6385897 -6.49904786 0.0000000
#Male:Light-Female:Light      1.2134144 -0.3563565  2.78318536 0.2185439
#Female:None-Female:Light    -1.4207568 -2.9905278  0.14901410 0.0974193
#Male:None-Female:Light      -1.6176328 -3.1874037 -0.04786184 0.0398106
#Female:None-Male:Light      -2.6341713 -4.2039422 -1.06440032 0.0001050
#Male:None-Male:Light        -2.8310472 -4.4008181 -1.26127627 0.0000284
#Male:None-Female:None       -0.1968759 -1.7666469  1.37289500 0.9990364

La valeur p indique s’il existe ou non une différence statistiquement significative entre chaque groupe.

Par exemple, dans la dernière ligne ci-dessus, nous voyons que le groupe d’hommes sans exercice n’a pas connu de différence statistiquement significative en termes de perte de poids par rapport au groupe de femmes sans exercice (valeur p : 0,990364).

Nous pouvons également visualiser les intervalles de confiance à 95 % résultant du test de Tukey en utilisant la fonction plot() dans R :

#set axis margins so labels don't get cut off
par(mar=c(4.1, 13, 4.1, 2.1))

#create confidence interval for each comparison
plot(TukeyHSD(model, conf.level=.95), las = 2)

Rapport des résultats de l’ANOVA bidirectionnelle

Enfin, nous pouvons rapporter les résultats de l’ANOVA bidirectionnelle de manière à résumer les résultats :

Une ANOVA bidirectionnelle a été réalisée pour examiner les effets du sexe ( homme, femme) et du programme d’exercice (aucun, léger, intense) sur la perte de poids (mesure en livres). Il y avait une interaction statistiquement significative entre les effets du sexe et de l’exercice sur la perte de poids (F(2, 54) = 4,615, p = 0,0141). Des tests post-hoc HSD de Tukey ont été effectués.

Pour les hommes, un programme d’exercices intenses entraîne une perte de poids significativement plus élevée qu’un programme léger (p < 0,0001) ou qu’aucun programme d’exercice (p < 0,0001). De plus, chez les hommes, un régime léger entraîne une perte de poids significativement plus élevée que l’absence de régime d’exercice (p < 0,0001).

Pour les femmes, un programme d’exercices intenses entraîne une perte de poids significativement plus élevée qu’un programme léger (p < 0,0001) ou qu’aucun programme d’exercice (p < 0,0001).

Des contrôles de normalité et le test de Levene ont été effectués pour vérifier que les hypothèses de l’ANOVA étaient respectées.