ANOVA à deux facteurs en R

Two-way ANOVA in R

Apprenez comment effectuer une ANOVA à deux facteurs dans R. Vous apprendrez également son objectif, ses hypothèses, ses prémisses et comment interpréter les résultats.

Introduction

L’ANOVA à deux facteurs (analyse de variance) est une méthode statistique qui permet d’évaluer l’effet simultané de deux variables catégorielles sur une variable continue quantitative.

L’ANOVA à deux facteurs est une extension de l’ANOVA à une voie car elle permet d’évaluer les effets sur une réponse numérique de deux variables catégorielles au lieu d’une seule. L’avantage d’une ANOVA à deux facteurs par rapport à une ANOVA à une voie est que nous testons la relation entre deux variables, tout en tenant compte de l’effet d’une troisième variable. De plus, elle permet également d’inclure l’interaction possible des deux variables catégorielles sur la réponse.

L’avantage d’une ANOVA à deux facteurs par rapport à une ANOVA à une voie est assez similaire à l’avantage d’une corrélation par rapport à une régression linéaire multiple :

• La corrélation mesure la relation entre deux variables quantitatives. La régression linéaire multiple mesure également la relation entre deux variables, mais cette fois en tenant compte de l’effet potentiel d’autres covariables.
• L’ANOVA à une voie teste si une variable quantitative est différente entre les groupes. L’ANOVA à deux facteurs teste également si une variable quantitative est différente entre les groupes, mais cette fois en tenant compte de l’effet d’une autre variable qualitative.

Auparavant, nous avons discuté de l’ANOVA à une voie dans R. Maintenant, nous montrons quand, pourquoi et comment effectuer une ANOVA à deux facteurs dans R.

• Comprendre les lois sur la propriété intellectuelle dans le trading algorithmique et l’IA en finance.
• Les logiciels d’IA dans les centres d’appels révolutionnent le service client.
• Rapport McKinsey Qu’est-ce que cela signifie pour les ventes B2B ?

Avant d’aller plus loin, j’aimerais mentionner et décrire brièvement quelques méthodes statistiques et tests connexes afin d’éviter toute confusion :

Un test t de Student est utilisé pour évaluer l’effet d’une variable catégorielle sur une variable continue quantitative, lorsque la variable catégorielle a exactement 2 niveaux :

• Test t de Student pour échantillons indépendants si les observations sont indépendantes (par exemple : si nous comparons l’âge entre les femmes et les hommes)
• Test t de Student pour échantillons appariés si les observations sont dépendantes, c’est-à-dire lorsqu’elles sont appariées (c’est le cas lorsque les mêmes sujets sont mesurés deux fois, à deux moments différents dans le temps, avant et après un traitement par exemple)

Pour évaluer l’effet d’une variable catégorielle sur une variable quantitative, lorsque la variable catégorielle a 3 ou plusieurs niveaux :

• ANOVA à une voie (souvent simplement appelée ANOVA) si les groupes sont indépendants (par exemple un groupe de patients ayant reçu le traitement A, un autre groupe de patients ayant reçu le traitement B et le dernier groupe de patients n’ayant reçu aucun traitement ou un placebo)
• ANOVA à mesures répétées si les groupes sont dépendants (lorsque les mêmes sujets sont mesurés trois fois, à trois moments différents dans le temps, avant, pendant et après un traitement par exemple)

Une ANOVA à deux facteurs est utilisée pour évaluer les effets de 2 variables catégorielles (et leur interaction potentielle) sur une variable continue quantitative. C’est le sujet du post.

La régression linéaire est utilisée pour évaluer la relation entre une variable dépendante continue quantitative et une ou plusieurs variables indépendantes :

• Régression linéaire simple s’il n’y a qu’une seule variable indépendante (qui peut être quantitative ou qualitative)
• Régression linéaire multiple s’il y a au moins deux variables indépendantes (qui peuvent être quantitatives, qualitatives ou un mélange des deux)

Une ANCOVA (analyse de covariance) est utilisée pour évaluer l’effet d’une variable catégorielle sur une variable quantitative, tout en contrôlant l’effet d’une autre variable quantitative (connue sous le nom de covariable). ANCOVA est en fait un cas particulier de la régression linéaire multiple avec un mélange d’une variable indépendante qualitative et d’une variable indépendante quantitative.

Dans ce post, nous commençons par expliquer quand et pourquoi une ANOVA à deux facteurs est utile, nous effectuons ensuite des analyses descriptives préliminaires et présentons comment effectuer une ANOVA à deux facteurs dans R. Enfin, nous montrons comment interpréter et visualiser les résultats. Nous mentionnons également brièvement et illustrons comment vérifier les hypothèses sous-jacentes.

Données

Pour illustrer comment effectuer une ANOVA à deux facteurs avec R, nous utilisons le jeu de données penguins, disponible dans le package {palmerpenguins}.

Nous n’avons pas besoin d’importer le jeu de données, mais nous devons d’abord charger le package et ensuite appeler le jeu de données:

# install.packages("palmerpenguins")library(palmerpenguins)

dat <- penguins # renommer le jeustr(dat) # structure du jeu de données

## tibble [344 × 8] (S3: tbl_df/tbl/data.frame)##  $ species          : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...##  $ island           : Factor w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ...##  $ bill_length_mm   : num [1:344] 39.1 39.5 40.3 NA 36.7 39.3 38.9 39.2 34.1 42 ...##  $ bill_depth_mm    : num [1:344] 18.7 17.4 18 NA 19.3 20.6 17.8 19.6 18.1 20.2 ...##  $ flipper_length_mm: int [1:344] 181 186 195 NA 193 190 181 195 193 190 ...##  $ body_mass_g      : int [1:344] 3750 3800 3250 NA 3450 3650 3625 4675 3475 4250 ...##  $ sex              : Factor w/ 2 levels "female","male": 2 1 1 NA 1 2 1 2 NA NA ...##  $ year             : int [1:344] 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...

Le jeu de données contient 8 variables pour 344 manchots, résumées ci-dessous:

summary(dat)

##       species          island    bill_length_mm  bill_depth_mm  ##  Adelie   :152   Biscoe   :168   Min.   :32.10   Min.   :13.10  ##  Chinstrap: 68   Dream    :124   1st Qu.:39.23   1st Qu.:15.60  ##  Gentoo   :124   Torgersen: 52   Median :44.45   Median :17.30  ##                                  Mean   :43.92   Mean   :17.15  ##                                  3rd Qu.:48.50   3rd Qu.:18.70  ##                                  Max.   :59.60   Max.   :21.50  ##                                  NA's   :2       NA's   :2      ##  flipper_length_mm  body_mass_g       sex           year     ##  Min.   :172.0     Min.   :2700   female:165   Min.   :2007  ##  1st Qu.:190.0     1st Qu.:3550   male  :168   1st Qu.:2007  ##  Median :197.0     Median :4050   NA's  : 11   Median :2008  ##  Mean   :200.9     Mean   :4202                Mean   :2008  ##  3rd Qu.:213.0     3rd Qu.:4750                3rd Qu.:2009  ##  Max.   :231.0     Max.   :6300                Max.   :2009  ##  NA's   :2         NA's   :2

Dans ce post, nous nous concentrons sur les trois variables suivantes:

• species : l’espèce du manchot (Adelie, Chinstrap ou Gentoo)
• sex : le sexe du manchot (femelle ou mâle)
• body_mass_g : la masse corporelle du manchot (en grammes)

Si nécessaire, plus d’informations sur ce jeu de données peuvent être trouvées en exécutant ?penguins dans R.

body_mass_g est la variable quantitative continue et sera la variable dépendante, tandis que species et sex sont toutes deux des variables qualitatives.

Ces deux dernières variables seront nos variables indépendantes, également appelées facteurs. Assurez-vous qu’elles sont lues comme des facteurs par R. Si ce n’est pas le cas, elles devront être transformées en facteurs.

Objectif et hypothèses d’une ANOVA à deux facteurs

Comme mentionné précédemment, une ANOVA à deux facteurs est utilisée pour évaluer simultanément l’effet de deux variables catégorielles sur une variable quantitative continue.

Elle est appelée ANOVA à deux facteurs car nous comparons des groupes formés par deux variables catégorielles indépendantes.

Ici, nous voulons savoir si la masse corporelle dépend de l’espèce et/ou du sexe. En particulier, nous nous intéressons à :

1. mesurer et tester la relation entre l’espèce et la masse corporelle,
2. mesurer et tester la relation entre le sexe et la masse corporelle, et
3. vérifier éventuellement si la relation entre l’espèce et la masse corporelle est différente pour les femelles et les mâles (ce qui équivaut à vérifier si la relation entre le sexe et la masse corporelle dépend de l’espèce)

Les deux premières relations sont appelées effets principaux, tandis que le troisième point est connu sous le nom d’effet d’interaction.

Les effets principaux testent si au moins un groupe est différent d’un autre (tout en contrôlant l’autre variable indépendante). D’autre part, l’effet d’interaction vise à tester si la relation entre deux variables diffère en fonction du niveau d’une troisième variable.

Lors de la réalisation d’une ANOVA à deux facteurs, il n’est pas obligatoire de tester l’effet d’interaction. Cependant, omettre un effet d’interaction peut conduire à des conclusions erronées si l’effet d’interaction est présent.

Si nous revenons à notre exemple, nous avons les tests d’hypothèses suivants :

Effet principal du sexe sur la masse corporelle :

• H0 : la masse corporelle moyenne est égale entre les femelles et les mâles
• H1 : la masse corporelle moyenne est différente entre les femelles et les mâles

Effet principal de l’espèce sur la masse corporelle :

• H0 : la masse corporelle moyenne est égale pour les 3 espèces
• H1 : la masse corporelle moyenne est différente pour au moins une espèce

Interaction entre le sexe et l’espèce :

• H0 : il n’y a pas d’interaction entre le sexe et l’espèce, ce qui signifie que la relation entre l’espèce et la masse corporelle est la même pour les femelles et les mâles (de même, la relation entre le sexe et la masse corporelle est la même pour les 3 espèces)
• H1 : il y a une interaction entre le sexe et l’espèce, ce qui signifie que la relation entre l’espèce et la masse corporelle est différente pour les femelles que pour les mâles (de même, la relation entre le sexe et la masse corporelle dépend de l’espèce)

Prémisses d’une ANOVA à deux facteurs

La plupart des tests statistiques nécessitent certaines hypothèses pour que les résultats soient valides, et une ANOVA à deux facteurs ne fait pas exception.

Les hypothèses d’une ANOVA à deux facteurs sont similaires à celles d’une ANOVA à un facteur. Pour résumer :

• Type de variable : la variable dépendante doit être quantitative continue, tandis que les deux variables indépendantes doivent être catégorielles (avec au moins deux niveaux).
• Indépendance : les observations doivent être indépendantes entre les groupes et à l’intérieur de chaque groupe.
• Normalité :
• Pour les petits échantillons, les données doivent suivre approximativement une distribution normale
• Pour les grands échantillons (généralement n ≥ 30 dans chaque groupe/échantillon), la normalité n’est pas requise (grâce au théorème central limite)
• Égalité des variances : les variances doivent être égales entre les groupes.
• Outliers : Il ne devrait pas y avoir d’outliers significatifs dans aucun groupe.

Plus de détails sur ces hypothèses peuvent être trouvés dans les hypothèses d’une ANOVA à un facteur.

Maintenant que nous avons vu les hypothèses sous-jacentes de l’ANOVA à deux facteurs, nous les examinons spécifiquement pour notre ensemble de données avant d’appliquer le test et d’interpréter les résultats.

Type de variable

La variable dépendante masse corporelle est quantitative continue, tandis que les deux variables indépendantes sexe et espèce sont des variables qualitatives (avec au moins 2 niveaux).

Par conséquent, cette hypothèse est respectée.

Indépendance

L’indépendance est généralement vérifiée en fonction de la conception de l’expérience et de la manière dont les données ont été collectées.

Pour simplifier, les observations sont généralement:

• indépendantes si chaque unité expérimentale (ici un manchot) n’a été mesurée qu’une seule fois et que les observations sont collectées à partir d’une partie représentative et sélectionnée de manière aléatoire de la population, ou
• dépendantes si chaque unité expérimentale a été mesurée au moins deux fois (comme c’est souvent le cas dans le domaine médical, par exemple, avec deux mesures sur les mêmes sujets; une avant et une après le traitement).

Dans notre cas, la masse corporelle n’a été mesurée qu’une seule fois sur chaque manchot, et sur un échantillon représentatif et aléatoire de la population, donc l’hypothèse d’indépendance est respectée.

Normalité

Nous avons un grand échantillon dans tous les sous-groupes (chaque combinaison des niveaux des deux facteurs, appelée cellule):

table(dat$species, dat$sex)

##            ##             female male##   Adelie        73   73##   Chinstrap     34   34##   Gentoo        58   61

donc la normalité n’a pas besoin d’être vérifiée.

Pour plus de précision, nous montrons encore comment vérifier la normalité, comme si nous avions de petits échantillons.

Il existe plusieurs méthodes pour tester l’hypothèse de normalité. Les méthodes les plus courantes étant:

• un graphique QQ par groupe ou sur les résidus, et/ou
• un histogramme par groupe ou sur les résidus, et/ou
• un test de normalité (par exemple le test de Shapiro-Wilk) par groupe ou sur les résidus.

La manière la plus simple/la plus courte de vérifier la normalité est de vérifier la normalité avec un graphique QQ sur les résidus. Pour dessiner ce graphique, nous devons d’abord sauvegarder le modèle:

# save modelmod <- aov(body_mass_g ~ sex * species,  data = dat)

Ce morceau de code sera expliqué plus en détail.

Maintenant, nous pouvons dessiner le graphique QQ sur les résidus. Nous montrons deux façons de le faire, d’abord avec la fonction plot() et deuxièmement avec la fonction qqPlot() du package {car}:

# méthode 1plot(mod, which = 2)

# méthode 2library(car)

qqPlot(mod$residuals,  id = FALSE # remove point identification)

Le code pour la méthode 1 est légèrement plus court, mais il manque l’intervalle de confiance autour de la ligne de référence.

Si les points suivent la ligne droite (appelée ligne de Henry) et tombent dans la bande de confiance, nous pouvons supposer la normalité. C’est le cas ici.

Si vous préférez vérifier la normalité sur la base d’un histogramme des résidus, voici le code:

# histogrammehist(mod$residuals)

L’histogramme des résidus montre une distribution gaussienne, ce qui est en ligne avec la conclusion du graphique QQ.

Bien que le graphique QQ et l’histogramme soient largement suffisants pour vérifier la normalité, si vous souhaitez le tester plus formellement avec un test statistique, le test de Shapiro-Wilk peut également être appliqué sur les résidus:

# test de normalitéshapiro.test(mod$residuals)

## ##  Shapiro-Wilk normality test## ## data:  mod$residuals## W = 0.99776, p-value = 0.9367

⇒ Nous ne rejetons pas l’hypothèse nulle selon laquelle les résidus suivent une distribution normale (p-value = 0,937).

À partir du graphique QQ-plot, de l’histogramme et du test de Shapiro-Wilk, nous concluons que nous ne rejetons pas l’hypothèse nulle de normalité des résidus.

L’hypothèse de normalité est donc vérifiée, nous pouvons maintenant vérifier l’égalité des variances.

Homogénéité des variances

L’égalité des variances, également appelée homogénéité des variances ou homoscédasticité, peut être vérifiée visuellement avec la fonction plot():

plot(mod, which = 3)

Comme l’étalement des résidus est constant, la ligne lisse rouge est horizontale et plate, il semble donc que l’hypothèse de variance constante soit satisfaite ici.

Le graphique de diagnostic ci-dessus est suffisant, mais si vous préférez, il peut également être testé plus formellement avec le test de Levene (également du package {car}):

leveneTest(mod)

## Test de Levene pour l'homogénéité des variances (center = médiane)##        Df valeur F Pr(>F)## group   5  1.3908 0.2272##       327

⇒ Nous ne rejetons pas l’hypothèse nulle selon laquelle les variances sont égales (p-valeur = 0,227).

Les approches visuelles et formelles donnent la même conclusion ; nous ne rejetons pas l’hypothèse d’homogénéité des variances.

Outliers

La manière la plus simple et la plus courante de détecter les outliers est visuelle grâce aux boxplots par groupes.

Pour les femelles et les mâles :

library(ggplot2)

# boxplots par sexe
ggplot(dat) +  aes(x = sex, y = body_mass_g) +  geom_boxplot()

Pour les trois espèces :

# boxplots par espèces
ggplot(dat) +  aes(x = species, y = body_mass_g) +  geom_boxplot()

Il y a, selon le critère de l’écart interquartile, deux outliers pour l’espèce Chinstrap. Ces points ne sont cependant pas assez extrêmes pour biaiser les résultats.

Par conséquent, nous considérons que l’hypothèse d’absence d’outliers significatifs est respectée.

ANOVA à deux facteurs

Nous avons montré que toutes les hypothèses sont vérifiées, nous pouvons donc maintenant procéder à la mise en œuvre de l’ANOVA à deux facteurs dans R.

Cela nous permettra de répondre aux questions de recherche suivantes :

• En contrôlant l’espèce, la masse corporelle est-elle significativement différente entre les deux sexes ?
• En contrôlant le sexe, la masse corporelle est-elle significativement différente pour au moins une espèce ?
• La relation entre l’espèce et la masse corporelle est-elle différente entre les pingouins femelles et mâles ?

Analyses préliminaires

Avant de réaliser tout test statistique, il est bon de faire des statistiques descriptives afin d’avoir un premier aperçu des données, et peut-être, avoir un aperçu des résultats à attendre.

Cela peut être fait via des statistiques descriptives ou des graphiques.

Statistiques descriptives

Si nous voulons simplifier, nous pouvons calculer seulement la moyenne pour chaque sous-groupe :

# moyenne par groupe
aggregate(body_mass_g ~ species + sex,  data = dat,  FUN = mean)

##     species    sex body_mass_g## 1    Adelie female    3368.836## 2 Chinstrap female    3527.206## 3    Gentoo female    4679.741## 4    Adelie   male    4043.493## 5 Chinstrap   male    3938.971## 6    Gentoo   male    5484.836

Ou éventuellement, la moyenne et l’écart-type pour chaque sous-groupe en utilisant le package {dplyr} :

# moyenne et écart-type par groupebibliothèque(dplyr)

group_by(dat, sex, species) %>%  summarise(    mean = round(mean(body_mass_g, na.rm = TRUE)),    sd = round(sd(body_mass_g, na.rm = TRUE))  )

## # A tibble: 8 × 4## # Groups:   sex [3]##   sex    species    mean    sd##   <fct>  <fct>     <dbl> <dbl>## 1 female Adelie     3369   269## 2 female Chinstrap  3527   285## 3 female Gentoo     4680   282## 4 male   Adelie     4043   347## 5 male   Chinstrap  3939   362## 6 male   Gentoo     5485   313## 7 <NA>   Adelie     3540   477## 8 <NA>   Gentoo     4588   338

Graphiques

Si vous êtes un lecteur régulier du blog, vous savez que j’aime dessiner des graphiques pour visualiser les données avant d’interpréter les résultats d’un test.

Le graphique le plus approprié lorsque nous avons une variable quantitative et deux variables qualitatives est un diagramme en boîte par groupe. Cela peut facilement être réalisé avec le package {ggplot2} :

# diagramme en boîte par groupebibliothèque(ggplot2)

ggplot(dat) +  aes(x = species, y = body_mass_g, fill = sex) +  geom_boxplot()

Certaines observations sont manquantes pour le sexe, nous pouvons les supprimer pour avoir un graphique plus concis :

dat %>%  filter(!is.na(sex)) %>%  ggplot() +  aes(x = species, y = body_mass_g, fill = sex) +  geom_boxplot()

Notez que nous aurions également pu réaliser le graphique suivant :

dat %>%  filter(!is.na(sex)) %>%  ggplot() +  aes(x = sex, y = body_mass_g, fill = species) +  geom_boxplot()

Mais pour un graphique plus lisible, je préfère mettre la variable avec le plus petit nombre de niveaux en couleur (qui est en fait l’argument fill dans la couche aes()) et la variable avec le plus grand nombre de catégories sur l’axe des x (c’est-à-dire l’argument x dans la couche aes()).

À partir des moyennes et des diagrammes en boîte par sous-groupe, nous pouvons déjà voir que, dans notre échantillon :

• les manchots femelles ont tendance à avoir une masse corporelle inférieure à celle des mâles, et c’est le cas pour toutes les espèces considérées, et
• la masse corporelle est plus élevée pour les manchots de Gentoo que pour les deux autres espèces.

Gardez à l’esprit que ces conclusions ne sont valables que dans notre échantillon ! Pour généraliser ces conclusions à la population, nous devons effectuer l’ANOVA à deux facteurs et vérifier la signification des variables explicatives. C’est l’objectif de la prochaine section.

ANOVA à deux facteurs en R

Comme mentionné précédemment, inclure un effet d’interaction dans une ANOVA à deux facteurs n’est pas obligatoire. Cependant, afin d’éviter des conclusions erronées, il est recommandé de vérifier d’abord si l’interaction est significative ou non, et en fonction des résultats, de l’inclure ou non.

Si l’interaction n’est pas significative, il est sûr de la supprimer du modèle final. Au contraire, si l’interaction est significative, elle doit être incluse dans le modèle final qui sera utilisé pour interpréter les résultats.

Nous commençons donc avec un modèle qui inclut les deux effets principaux (c’est-à-dire le sexe et l’espèce) et l’interaction:

# ANOVA à deux facteurs avec interaction# enregistrer le modèlemod <- aov(body_mass_g ~ sex * species,  data = dat)

# afficher les résultatssummary(mod)

##              Df    Sum Sq  Mean Sq F value   Pr(>F)    ## sex           1  38878897 38878897 406.145  < 2e-16 ***## species       2 143401584 71700792 749.016  < 2e-16 ***## sex:species   2   1676557   838278   8.757 0.000197 ***## Residuals   327  31302628    95727                     ## ---## Signif. codes:  0 '***' 0,001 '**' 0,01 '*' 0,05 '.' 0,1 ' ' 1## 11 observations deleted due to missingness

La somme des carrés (colonne Sum Sq) montre que l’espèce explique une grande partie de la variabilité de la masse corporelle. C’est le facteur le plus important pour expliquer cette variabilité.

Les valeurs p sont affichées dans la dernière colonne de la sortie ci-dessus (Pr(>F)). À partir de ces valeurs p, nous concluons que, au niveau de signification de 5%:

• en contrôlant pour l’espèce, la masse corporelle est significativement différente entre les deux sexes,
• en contrôlant pour le sexe, la masse corporelle est significativement différente pour au moins une espèce, et
• l’interaction entre le sexe et l’espèce (affichée à la ligne sex:species dans la sortie ci-dessus) est significative.

Ainsi, grâce à l’effet d’interaction significatif, nous venons de voir que la relation entre la masse corporelle et l’espèce est différente entre les mâles et les femelles. Comme c’est significatif, nous devons le conserver dans le modèle et nous devrions interpréter les résultats à partir de ce modèle.

Si, au contraire, l’interaction n’était pas significative (c’est-à-dire si la valeur p ≥ 0,05), nous aurions supprimé cet effet d’interaction du modèle. À des fins d’illustration, ci-dessous le code pour une ANOVA à deux facteurs sans interaction, appelée modèle additif :

# ANOVA à deux facteurs sans interactionaov(body_mass_g ~ sex + species,  data = dat)

Pour les lecteurs habitués à effectuer des régressions linéaires en R, vous remarquerez que la structure du code pour une ANOVA à deux facteurs est en fait similaire :

• la formule est variable dépendante ~ variables indépendantes
• le signe + est utilisé pour inclure des variables indépendantes sans interaction
• le signe * est utilisé pour inclure des variables indépendantes avec une interaction

La ressemblance avec une régression linéaire n’est pas une surprise car une ANOVA à deux facteurs, comme toutes les ANOVA, est en fait un modèle linéaire.

Remarquez que le code suivant fonctionne également et donne les mêmes résultats :

# méthode 2mod2 <- lm(body_mass_g ~ sex * species,  data = dat)

Anova(mod2)

## Table ANOVA (tests de type II)## ## Réponse : body_mass_g##                Sum Sq  Df F value    Pr(>F)    ## sex          37090262   1 387.460 < 2.2e-16 ***## species     143401584   2 749.016 < 2.2e-16 ***## sex:species   1676557   2   8.757 0.0001973 ***## Residuals    31302628 327                      ## ---## Signif. codes:  0 '***' 0,001 '**' 0,01 '*' 0,05 '.' 0,1 ' ' 1

Remarquez que la fonction aov() suppose une conception équilibrée, c’est-à-dire que nous avons des tailles d’échantillon égales dans les niveaux de nos variables de regroupement indépendantes.

Pour les conceptions non équilibrées, c’est-à-dire des nombres inégaux de sujets dans chaque sous-groupe, les méthodes recommandées sont :

• la ANOVA de type II lorsqu’il n’y a pas d’interaction significative, qui peut être effectuée en R avec Anova(mod, type = "II"), 5 et
• la ANOVA de type III lorsqu’il y a une interaction significative, qui peut être effectuée en R avec Anova(mod, type = "III").

Cela dépasse le cadre du post et nous supposons ici une conception équilibrée. Pour le lecteur intéressé, voir cette discussion détaillée sur les ANOVA de type I, II et III.

Comparaisons par paires

À travers les deux effets principaux étant significatifs, nous avons conclu que :

• en contrôlant pour l’espèce, la masse corporelle est différente entre les femelles et les mâles, et
• en contrôlant pour le sexe, la masse corporelle est différente pour au moins une espèce.

Si la masse corporelle est différente entre les deux sexes, étant donné qu’il y a exactement deux sexes, cela doit être parce que la masse corporelle est significativement différente entre les femelles et les mâles.

Si l’on veut savoir quel sexe a la plus grande masse corporelle, cela peut être déduit des moyennes et/ou des boxplots par sous-groupe. Ici, il est clair que les mâles ont une masse corporelle significativement plus élevée que les femelles.

Cependant, ce n’est pas aussi simple pour l’espèce. Laissez-moi expliquer pourquoi ce n’est pas aussi facile que pour les sexes.

Il y a trois espèces (Adélie, Chinstrap et Gentoo), donc il y a 3 paires d’espèces :

1. Adélie et Chinstrap
2. Adélie et Gentoo
3. Chinstrap et Gentoo

Si la masse corporelle est significativement différente pour au moins une espèce, cela pourrait être que :

• la masse corporelle est significativement différente entre Adélie et Chinstrap mais pas significativement différente entre Adélie et Gentoo, et non significativement différente entre Chinstrap et Gentoo, ou
• la masse corporelle est significativement différente entre Adélie et Gentoo mais pas significativement différente entre Adélie et Chinstrap, et non significativement différente entre Chinstrap et Gentoo, ou
• la masse corporelle est significativement différente entre Chinstrap et Gentoo mais pas significativement différente entre Adélie et Chinstrap, et non significativement différente entre Adélie et Gentoo.

Ou, il pourrait également être que :

• la masse corporelle est significativement différente entre Adélie et Chinstrap, et entre Adélie et Gentoo, mais pas significativement différente entre Chinstrap et Gentoo, ou
• la masse corporelle est significativement différente entre Adélie et Chinstrap, et entre Chinstrap et Gentoo, mais pas significativement différente entre Adélie et Gentoo, ou
• la masse corporelle est significativement différente entre Chinstrap et Gentoo, et entre Adélie et Gentoo, mais pas significativement différente entre Adélie et Chinstrap.

Enfin, il pourrait également être que la masse corporelle est significativement différente entre toutes les espèces.

Comme pour une ANOVA à un facteur, nous ne pouvons pas, à ce stade, savoir précisément quelle espèce est différente de laquelle en termes de masse corporelle. Pour le savoir, nous devons comparer chaque espèce deux par deux grâce à des tests post-hoc (également appelés comparaisons par paires).

Il existe plusieurs tests post-hoc, le plus courant étant le test de Tukey HSD, qui teste toutes les paires possibles de groupes. Comme mentionné précédemment, ce test ne doit être effectué que sur la variable d’espèce car il n’y a que deux niveaux pour le sexe.

Comme pour l’ANOVA à un facteur, le test de Tukey HSD peut être effectué en R comme suit :

# méthode 1TukeyHSD(mod, which = "species")

##   Tukey multiple comparisons of means##     95% family-wise confidence level## ## Fit: aov(formula = body_mass_g ~ sex * species, data = dat)## ## $species##                        diff       lwr       upr     p adj## Chinstrap-Adelie   26.92385  -80.0258  133.8735 0.8241288## Gentoo-Adelie    1377.65816 1287.6926 1467.6237 0.0000000## Gentoo-Chinstrap 1350.73431 1239.9964 1461.4722 0.0000000

ou en utilisant le package {multcomp}:

# méthode 2library(multcomp)

summary(glht(  aov(masse_corporelle_g ~ sexe + espèce,    data = dat  ),  linfct = mcp(espèce = "Tukey")))

## ##   Tests simultanés pour des hypothèses linéaires générales## ## Comparaisons multiples des moyennes : Contrastes de Tukey## ## ## Ajustement de la donnée : aov(formule = masse_corporelle_g ~ sexe + espèce, data = dat)## ## Hypothèses linéaires :##                         Estimate Std. Error t value Pr(>|t|)    ## Manchot de Chinstrap - Manchot d'Adélie == 0    26,92      46,48   0,579     0,83    ## Manchot papou - Manchot d'Adélie == 0     1377,86      39,10  35,236   <1e-05 ***## Manchot papou - Manchot de Chinstrap == 0  1350,93      48,13  28,067   <1e-05 ***## ---## Codes signif. :  0 '***' 0,001 '**' 0,01 '*' 0,05 '.' 0,1 ' ' 1## (Valeurs p ajustées rapportées - méthode à étape unique)

ou en utilisant la fonction pairwise.t.test() avec la méthode d’ajustement du p-value de votre choix :

# méthode 3pairwise.t.test(dat$masse_corporelle_g, dat$espèce,  p.adjust.method = "BH")

## ## Comparaisons deux à deux utilisant des tests t avec SD groupés## ## données :  dat$masse_corporelle_g et dat$espèce ## ##           Manchot d'Adélie Manchot de Chinstrap## Manchot de Chinstrap 0,63   -        ## Manchot papou    <2e-16 <2e-16   ## ## Méthode d'ajustement de la valeur p : BH

Notez que lorsque vous utilisez la seconde méthode, c’est le modèle sans interaction qui doit être spécifié dans la fonction glht(), même si l’interaction est significative. De plus, n’oubliez pas de remplacer mod et espèce dans mon code par le nom de votre modèle et le nom de votre variable indépendante.

Les deux méthodes donnent les mêmes résultats, à savoir :

• la masse corporelle n’est pas significativement différente entre le Manchot de Chinstrap et le Manchot d’Adélie (p-valeur ajustée = 0,83),
• la masse corporelle est significativement différente entre le Manchot papou et le Manchot d’Adélie (p-valeur ajustée < 0,001), et
• la masse corporelle est significativement différente entre le Manchot papou et le Manchot de Chinstrap (p-valeur ajustée < 0,001).

N’oubliez pas que ce sont les valeurs de p ajustées qui sont rapportées, pour éviter le problème des tests multiples qui se produit lors de la comparaison de plusieurs paires de groupes.

Si vous souhaitez comparer toutes les combinaisons de groupes, cela peut être fait avec la fonction TukeyHSD() en spécifiant l’interaction dans l’argument which :

# toutes les combinaisons de sexe et d'espèceTukeyHSD(mod,  which = "sexe:espèce")

##   Comparaisons multiples de Tukey des moyennes##     Niveau de confiance global de 95%## ## Ajustement de la donnée : aov(formule = masse_corporelle_g ~ sexe * espèce, data = dat)## ## $`sexe:espèce`##                                      diff       lwr       upr     p adj## male:Adelie-female:Adelie        674,6575  527,8486  821,4664 0,0000000## female:Chinstrap-female:Adelie   158,3703  -25,7874  342,5279 0,1376213## male:Chinstrap-female:Adelie     570,1350  385,9773  754,2926 0,0000000## female:Gentoo-female:Adelie     1310,9058 1154,8934 1466,9181 0,0000000## male:Gentoo-female:Adelie       2116,0004 1962,1408 2269,8601 0,0000000## female:Chinstrap-male:Adelie    -516,2873 -700,4449 -332,1296 0,0000000## male:Chinstrap-male:Adelie      -104,5226 -288,6802   79,6351 0,5812048## female:Gentoo-male:Adelie        636,2482  480,2359  792,2606 0,0000000## male:Gentoo-male:Adelie         1441,3429 1287,4832 1595,2026 0,0000000## male:Chinstrap-female:Chinstrap  411,7647  196,6479  626,8815 0,0000012## female:Gentoo-female:Chinstrap  1152,5355  960,9603 1344,1107 0,0000000## male:Gentoo-female:Chinstrap    1957,6302 1767,8040 2147,4564 0,0000000## female:Gentoo-male:Chinstrap     740,7708  549,1956  932,3460 0,0000000## male:Gentoo-male:Chinstrap      1545,8655 1356,0392 1735,6917 0,0000000## male:Gentoo-female:Gentoo        805,0947  642,4300  967,7594 0,0000000

Ou avec la fonction HSD.test() du package {agricolae}, qui indique les sous-groupes qui ne sont pas significativement différents les uns des autres avec la même lettre :

library(agricolae)

HSD.test(mod,  trt = c("sex", "species"),  console = TRUE # imprimer les résultats)

## ## Étude : mod ~ c("sex", "species")## ## Test HSD pour body_mass_g ## ## Erreur quadratique moyenne : 95726.69 ## ## sex:species,  moyennes## ##                  body_mass_g      std  r  Min  Max## female:Adelie       3368.836 269.3801 73 2850 3900## female:Chinstrap    3527.206 285.3339 34 2700 4150## female:Gentoo       4679.741 281.5783 58 3950 5200## male:Adelie         4043.493 346.8116 73 3325 4775## male:Chinstrap      3938.971 362.1376 34 3250 4800## male:Gentoo         5484.836 313.1586 61 4750 6300## ## Alpha : 0.05 ; DF erreur : 327 ## Valeur critique de l'étendue studentisée : 4.054126 ## ## Les traitements avec la même lettre ne sont pas significativement différents.## ##                  body_mass_g groups## male:Gentoo         5484.836      a## female:Gentoo       4679.741      b## male:Adelie         4043.493      c## male:Chinstrap      3938.971      c## female:Chinstrap    3527.206      d## female:Adelie       3368.836      d

Si vous avez de nombreux groupes à comparer, les représenter graphiquement peut être plus facile à interpréter :

# définir les marges de l'axe pour que les étiquettes ne soient pas coupéespar(mar = c(4.1, 13.5, 4.1, 2.1))

# créer un intervalle de confiance pour chaque comparaisonplot(TukeyHSD(mod, which = "sex:species"),  las = 2 # faire pivoter les graduations de l'axe des x)

À partir des résultats et du graphique ci-dessus, nous concluons que toutes les combinaisons de sexe et d’espèce sont significativement différentes, sauf entre les femelles Chinstrap et Adélie (valeur p = 0,138) et les mâles Chinstrap et Adélie (valeur p = 0,581).

Ces résultats, qui sont d’ailleurs en accord avec les boîtes à moustaches montrées ci-dessus et qui seront confirmés avec les visualisations ci-dessous, concluent l’ANOVA à deux facteurs en R.

Visualisations

Si vous souhaitez visualiser les résultats d’une manière différente de ce qui a déjà été présenté dans les analyses préliminaires, voici quelques idées de graphiques utiles.

Tout d’abord, avec la moyenne et l’erreur standard de la moyenne par sous-groupe en utilisant la fonction allEffects() du package {effects} :

# méthode 1library(effects)

plot(allEffects(mod))

Ou en utilisant le package {ggpubr} :

# méthode 2library(ggpubr)

ggline(subset(dat, !is.na(sex)), # supprimer le niveau NA pour le sexe  x = "species",  y = "body_mass_g",  color = "sex",  add = c("mean_se") # ajouter la moyenne et l'erreur standard) +  labs(y = "Moyenne de la masse corporelle (g)")

En outre, en utilisant {Rmisc} et {ggplot2} :

library(Rmisc)

# calculer la moyenne et l'erreur standard de la moyenne par sous-groupesummary_stat <- summarySE(dat,  measurevar = "body_mass_g",  groupvars = c("species", "sex"))# tracer la moyenne et l'erreur standard de la moyenneggplot(  subset(summary_stat, !is.na(sex)), # supprimer le niveau NA pour le sexe  aes(x = species, y = body_mass_g, colour = sex)) +  geom_errorbar(aes(ymin = body_mass_g - se, ymax = body_mass_g + se), # ajouter des barres d'erreur    width = 0.1 # largeur des barres d'erreur  ) +  geom_point() +  labs(y = "Moyenne de la masse corporelle (g)")

Deuxièmement, si vous préférez ne tracer que la moyenne par sous-groupe :

with(  dat,  interaction.plot(species, sex, body_mass_g))

Enfin, pour ceux d’entre vous qui connaissent GraphPad, vous êtes probablement familiers avec la représentation des moyennes et des barres d’erreur comme suit :

# tracer la moyenne et l'erreur standard de la moyenne sous forme d'histogrammesggplot(  subset(summary_stat, !is.na(sex)), # supprimer le niveau NA pour le sexe  aes(x = species, y = body_mass_g, fill = sex)) +  geom_bar(position = position_dodge(), stat = "identity") +  geom_errorbar(aes(ymin = body_mass_g - se, ymax = body_mass_g + se), # ajouter des barres d'erreur    width = 0.25, # largeur des barres d'erreur    position = position_dodge(.9)  ) +  labs(y = "Moyenne de la masse corporelle (g)")

Conclusion

Dans cet article, nous avons commencé par rappeler les différents tests existants pour comparer une variable quantitative entre des groupes. Nous nous sommes ensuite concentrés sur l’analyse de variance à deux facteurs, en commençant par son objectif et ses hypothèses jusqu’à son implémentation dans R, ainsi que ses interprétations et certaines visualisations. Nous avons également brièvement mentionné ses hypothèses sous-jacentes et un test post-hoc pour comparer tous les sous-groupes.

Tout cela a été illustré avec l’ensemble de données penguins disponible dans le package {palmerpenguins}.

Merci de votre lecture.

J’espère que cet article vous aidera à réaliser une analyse de variance à deux facteurs avec vos données.

Comme toujours, si vous avez une question ou une suggestion concernant le sujet abordé dans cet article, veuillez l’ajouter en commentaire pour que d’autres lecteurs puissent en bénéficier.

1. En théorie, une analyse de variance à un facteur peut également être utilisée pour comparer 2 groupes, et pas seulement 3 ou plus. Néanmoins, en pratique, il est souvent le cas qu’un test t de Student est effectué pour comparer 2 groupes, et une analyse de variance à un facteur pour comparer 3 groupes ou plus. Les conclusions obtenues via un test t de Student pour des échantillons indépendants et une analyse de variance à un facteur avec 2 groupes seront similaires. ↩︎
2. Notez que si l’hypothèse de normalité n’est pas respectée, de nombreuses transformations peuvent être appliquées pour l’améliorer, la plus courante étant la transformation logarithmique (fonction log() en R). ↩︎
3. Notez que le test de Bartlett est également approprié pour tester l’hypothèse d’égalité des variances. ↩︎
4. Un modèle additif suppose que les 2 variables explicatives sont indépendantes ; elles n’interagissent pas entre elles. ↩︎
5. Où mod est le nom de votre modèle enregistré. ↩︎
6. Ici, nous utilisons la correction de Benjamini & Hochberg (1995), mais vous pouvez choisir entre plusieurs méthodes. Voir ?p.adjust pour plus de détails. ↩︎

Articles connexes

• ANOVA en R
• Test t de Student en R et à la main : comment comparer deux groupes dans différents scénarios ?
• Test de Wilcoxon à un échantillon en R
• Test de Kruskal-Wallis, ou la version non-paramétrique de l’ANOVA
• Test d’hypothèse à la main

Initialement publié sur https://statsandr.com le 19 juin 2023.

We will continue to update IPGirl; if you have any questions or suggestions, please contact us!