Carrés de type III et carrés de type II dans l'analyse de variance
L’analyse de la variance à deux facteurs est une approche statistique qui permet d’évaluer les effets de deux facteurs et de leur interaction sur une réponse quantitative continue. C’est une méthode qui est employée pour analyser les plans factoriels à deux facteurs.
Pour plus d’information sur cette approche, vous pouvez consulter cet article “Introduction à l’ANOVA à 2 facteurs” .
Le principe est de décomposer la variabilité des observations en 4 parties :
- une imputable à l’effet du 1er facteur
- une imputable à l’effet du 2ème facteur
- une imputable à l’effet de l’interaction des deux facteurs
- une non imputable, c’est la part résiduelle.
Pour estimer ces différentes parts, on calcule des sommes des carrés relatives à ces différentes sources (1er facteur, second facteur, interaction et résiduelle)
Pour plus d’informations sur ce calcul de la somme des carrés, vous pouvez consulter cet article : ANOVA à 2 facteurs : principe.
Lorsque le nombre d’observations n’est pas identique pour chaque combinaison des modalités des facteurs, on dit que les données sont déséquilibrées.
Et lorsque les données sont déséquilibrées, il existe différentes façons de calculer la sommes des carrés relatives aux différentes sources. On parle de carrés de type II ou III.
“Lors de la décomposition de la variabilité, ce déséquilibre a un impact sur le calcul des différentes sommes de carrés, et par la suite sur les tests sui sont construits. Il est donc nécessaire de choisir quel type de carrés seront utilisés (type I, type II, ou type III)”
Bien qu’il existe toujours un débat sur le choix du type de carrés (II ou III), l’approche généralement employée pour évaluer les effets des facteurs en présence d’une interaction consiste à calculer des carrés de type III.
En cas d’ajustement d‘un modèle ANOVA sans interaction, les carrés de type III sont égaux aux carrés de type II. La question du choix du type de carrés ne se pose donc pas.
Dans cet article je vais vous montrer comment obtenir ces carrés de types III (il existe différentes possibilités), lorsqu’un modèle complet, avec interaction, est ajusté.
Pour plus d’information sur la réalisation d’ANOVA à deux facteurs avec le logiciel R, vous pouvez consulter cet article ANOVA à 2 facteurs avec R : Tutoriel.
Table des matières
Les data
Pour illustrer cet article, j’ai créé, à partir des données penguins du package palmerpenguins deux jeux de données :
- un déséquilibré,
- et un équilibré.
Données déséquilibrées
Voici le code pour générer le jeu de données déséquilibrées :
#install.packages("palmerpenguins")
#install.packages("tidyverse")
library(palmerpenguins)
library(tidyverse)
penguins07 <- penguins %>%
filter(year==2007)
ind1 <- c(2,3,5,7,49,51,53,56,102,103,105,107)
mydata_DesEq <- penguins07[-ind1,]
Nous nous intéresserons aux effets des variables species et sex sur la variable body_mass_g.
Ces données sont déséquilibrées puisque les effectifs de chaque groupe au croisement des variables species et sex ne sont pas identiques :
table( mydata_DesEq$sex,mydata_DesEq$species)
##
## Adelie Chinstrap Gentoo
## female 17 9 13
## male 22 13 17 Données équilibrées
Voici le code pour créer les données équilibrées :
mydata_Eq <-penguins07[-c(2,3,5,7,13,16,17,19,21,1,6,8,14,15,18,20,22,24,51,53,56,52,54,55,58),]
Ces données sont équilibrées puisque les effectifs de chaque groupe au croisement des variables species et sex sont identiques :
table( mydata_Eq$sex,mydata_Eq$species)
##
## Adelie Chinstrap Gentoo
## female 13 13 13
## male 13 13 13 Analyse de la variance à deux facteurs avec interaction lorsque les données sont déséquilibrées
Obtention des carrés de type II
Nous allons commencer par calculer les carrés de type II, car il s’agit de l’approche par défaut Pour cela, il suffit d’ajuster le modèle (le terme sex*species permet de modéliser l’effet de chaque facteur et de leur interaction), et d’employer la fonction Anova() avec ses arguments par défaut.
mydata_DesEq_type2<- aov(body_mass_g~sex*species, data=mydata_DesEq)
library(car)
Anova(mydata_DesEq_type2)
## Anova Table (Type II tests)
##
## Response: body_mass_g
## Sum Sq Df F value Pr(>F)
## sex 9071259 1 90.651 4.761e-15 ***
## species 40705547 2 203.390 < 2.2e-16 ***
## sex:species 1667930 2 8.334 0.0004954 ***
## Residuals 8505754 85
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 La variable Sum Sq correspond à la somme des carrés (ici de type II).
Ici, l’interaction est significative.
Remarque : il est ensuite possible de réaliser des comparaisons multiples à l’intérieur des modalités d’un des deux facteurs (puisque l’interaction est significative), comme ceci :
library(emmeans)
emmeans(mydata_DesEq_type2, specs=pairwise~species|sex)
## $emmeans
## sex = female:
## species emmean SE df lower.CL upper.CL
## Adelie 3353 76.7 85 3200 3505
## Chinstrap 3639 105.4 85 3429 3849
## Gentoo 4646 87.7 85 4472 4821
##
## sex = male:
## species emmean SE df lower.CL upper.CL
## Adelie 4039 67.4 85 3905 4173
## Chinstrap 3819 87.7 85 3645 3994
## Gentoo 5553 76.7 85 5400 5705
##
## Confidence level used: 0.95
##
## $contrasts
## sex = female:
## contrast estimate SE df t.ratio p.value
## Adelie - Chinstrap -286 130 85 -2.193 0.0782
## Adelie - Gentoo -1293 117 85 -11.096 <.0001
## Chinstrap - Gentoo -1007 137 85 -7.343 <.0001
##
## sex = male:
## contrast estimate SE df t.ratio p.value
## Adelie - Chinstrap 219 111 85 1.983 0.1227
## Adelie - Gentoo -1514 102 85 -14.824 <.0001
## Chinstrap - Gentoo -1734 117 85 -14.875 <.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimates Obtention des carrés de type III
Les carrés de type III peuvent être obtenus de différentes façons. Nous allons commencer par modifier les contrastes des facteurs en amont de la modélisation.
En modifiant les contrastes des variables en amont
Pour illustrer cette étape, je cré une copie des données déséquilibrées dans un data frame nommé mydata_DesEqC` ; C pour “contrasts”.
Ensuite, je modifie les contrastes par défaut (qui sont de type contr.treatment), par des contrastes de type contr.sum :
mydata_DesEqC <- mydata_DesEq
contrasts(mydata_DesEqC$species)
## Chinstrap Gentoo
## Adelie 0 0
## Chinstrap 1 0
## Gentoo 0 1
contrasts(mydata_DesEqC$species) <- "contr.sum"
contrasts(mydata_DesEqC$species)
## [,1] [,2]
## Adelie 1 0
## Chinstrap 0 1
## Gentoo -1 -1
contrasts(mydata_DesEqC$sex)
## male
## female 0
## male 1
contrasts(mydata_DesEqC$sex) <- "contr.sum"
contrasts(mydata_DesEqC$sex)
## [,1]
## female 1
## male -1 Nous allons, à présent, ajuster le modèle complet de l’ANOVA à deux facteurs, avec l’interaction en employant le terme `sex*species` (le caractère * permettant de modéliser l’effet de chaque facteur et de leur interaction)
mydata_DesEq_type3 <- aov(body_mass_g~sex*species, data=mydata_DesEqC)
Pour obtenir les carrés de type III, il est encore nécessaire de le préciser type=3 en argument de la fonction Anova(), comme ceci :
Anova(mydata_DesEq_type3, type=3)
## Anova Table (Type III tests)
##
## Response: body_mass_g
## Sum Sq Df F value Pr(>F)
## (Intercept) 1465784351 1 14647.928 < 2.2e-16 ***
## sex 7342229 1 73.373 4.086e-13 ***
## species 38222607 2 190.984 < 2.2e-16 ***
## sex:species 1667930 2 8.334 0.0004954 ***
## Residuals 8505754 85
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Nous pouvons voir que les sommes des carrés (et par conséquent les F et p values) relatives aux variables sex et species sont un peu différentes (plus faibles) que celles obtenues précédemment avec les carrés de type II.
En revanche les sommes des carrés relatives à l’interaction des variables sex et species sont identiques.
Les comparaisons multiples peuvent ensuite être réalisées, et nous obtenons des résultats identiques à ceux du modèle mydata_DesEq.aov1 (sans modification des contrastes), vu précédemment.
Cela n’est pas étonnant puisqu’il s’agit de tests t, dans lesquels les contrastes employés dans l’ANOVA n’ont pas d’influence.
library(emmeans)
emmeans(mydata_DesEq_type3, specs=pairwise~species|sex)
## $emmeans
## sex = female:
## species emmean SE df lower.CL upper.CL
## Adelie 3353 76.7 85 3200 3505
## Chinstrap 3639 105.4 85 3429 3849
## Gentoo 4646 87.7 85 4472 4821
##
## sex = male:
## species emmean SE df lower.CL upper.CL
## Adelie 4039 67.4 85 3905 4173
## Chinstrap 3819 87.7 85 3645 3994
## Gentoo 5553 76.7 85 5400 5705
##
## Confidence level used: 0.95
##
## $contrasts
## sex = female:
## contrast estimate SE df t.ratio p.value
## Adelie - Chinstrap -286 130 85 -2.193 0.0782
## Adelie - Gentoo -1293 117 85 -11.096 <.0001
## Chinstrap - Gentoo -1007 137 85 -7.343 <.0001
##
## sex = male:
## contrast estimate SE df t.ratio p.value
## Adelie - Chinstrap 219 111 85 1.983 0.1227
## Adelie - Gentoo -1514 102 85 -14.824 <.0001
## Chinstrap - Gentoo -1734 117 85 -14.875 <.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimates Remarque :
Attention, si vous utilisez uniquement l’argument type = 3, dans la fonction Anova() sans modifier au préalable les contrastes, vous n’obtiendrez pas des carrés de type III.
En effet, si on fait cela avec le modèle mydata_DesEq_type2 (celui ajusté sur les données mydata_DesEq, et donc sans modification des contrastes), vous pouvez voir que les résultats sont différents :
Anova(mydata_DesEq_type2, type=3)
## Anova Table (Type III tests)
##
## Response: body_mass_g
## Sum Sq Df F value Pr(>F)
## (Intercept) 191117647 1 1909.884 < 2.2e-16 ***
## sex 4508885 1 45.058 2.026e-09 ***
## species 12841450 2 64.164 < 2.2e-16 ***
## sex:species 1667930 2 8.334 0.0004954 ***
## Residuals 8505754 85
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 En modifiant les contrasts dans le modèle de l’ANOVA
Pour illustrer cette approche nous repartons des données originelles mydata_DesEq.
Les contrastes peuvent être modifiés dans la fonction aov(), en employant l’argument contrasts=list(sex=contr.sum,species=contr.sum). Il est encore nécessaire d’utiliser ensuite la fonction Anova(), avec l’argument type=3.
mydata_DesEq_type3bis <- aov(body_mass_g~sex*species,
contrasts=list(sex=contr.sum,
species=contr.sum),
data=mydata_DesEq)
Anova(mydata_DesEq_type3bis, type=3)
## Anova Table (Type III tests)
##
## Response: body_mass_g
## Sum Sq Df F value Pr(>F)
## (Intercept) 1465784351 1 14647.928 < 2.2e-16 ***
## sex 7342229 1 73.373 4.086e-13 ***
## species 38222607 2 190.984 < 2.2e-16 ***
## sex:species 1667930 2 8.334 0.0004954 ***
## Residuals 8505754 85
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 On retrouve bien les résultats du modèle mydata_DesEq_type3 obtenus précédemment.
En utilisant la fonction AovSum() du package FactoMineR
Le package FactoMineR est surtout connu pour la réalisation d’ACP ( pour plus d’informations vous pouvez consulter cet article “Ressources pour l’ACP“) mais il contient aussi la fonction AovSum() qui permet d’obtenir directement des carrés de type III :
#install.packages("FactoMineR")
library(FactoMineR)
mydata_DesEq_aovsum <- AovSum(body_mass_g~sex*species, data=mydata_DesEq)
mydata_DesEq_aovsum$Ftest
## SS df MS F value Pr(>F)
## sex 7342229 1 7342229 73.373 4.086e-13 ***
## species 38222607 2 19111303 190.984 < 2.2e-16 ***
## sex:species 1667930 2 833965 8.334 0.0004954 ***
## Residuals 8505754 85 100068
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Nous obtenons, là encore, les mêmes résultats que précédemment avec la commande Anova(mydata_DesEq_type3, type=3).
Cette approche par la fonction AovSum() à l’avantage de sa simplicité, car nous n’avons pas à nous préoccuper de modifier les contrastes.
En revanche, il n’est alors plus possible d’employer la fonction emmeans() pour réaliser les comparaisons multiples ! C’est un peu dommage….
emmeans(mydata_DesEq_aovsum, specs=pairwise~species|sex)
## Error in (function (object, at, cov.reduce = mean, cov.keep = get_emm_option("cov.keep"), : Can't handle an object of class "AovSum"
## Use help("models", package = "emmeans") for information on supported models. En utilisant la fonction drop1()
Avant d’utiliser la fonction Anova() et son argument type=3 pour obtenir les carrés de type 3, j’utilisais la fonction drop1() ; elle conduit, là encore, aux mêmes résultats :
drop1(mydata_DesEq_type3bis, .~., test="F")
## Single term deletions
##
## Model:
## body_mass_g ~ sex * species
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 8505754 1053.5
## sex 1 7342229 15847983 1108.2 73.373 4.086e-13 ***
## species 2 38222607 46728361 1204.6 190.984 < 2.2e-16 ***
## sex:species 2 1667930 10173683 1065.8 8.334 0.0004954 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Egalité des carrés de type II et III en absence d’interaction
Il s’agit simplement ici, de montrer que le problème du choix des carrés de type II ou III ne se pose pas lorsque le modèle est ajusté sans interaction.
Pour commencer, un modèle sans interaction est ajusté, et des carrés de type II sont calculés :
mydata_DesEq_mod1<- aov(body_mass_g~sex+species, data=mydata_DesEq)
Anova(mydata_DesEq_mod1)
## Anova Table (Type II tests)
##
## Response: body_mass_g
## Sum Sq Df F value Pr(>F)
## sex 9071259 1 77.573 1.113e-13 ***
## species 40705547 2 174.046 < 2.2e-16 ***
## Residuals 10173683 87
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Puis le même modèle sans interaction est ajusté, mais cette fois en calculant des carrés de type III :
mydata_DesEq_mod2 <- aov(body_mass_g~sex+species,
contrasts=list(sex=contr.sum,
species=contr.sum),
data=mydata_DesEq)
Anova(mydata_DesEq_mod2, type=3)
## Anova Table (Type III tests)
##
## Response: body_mass_g
## Sum Sq Df F value Pr(>F)
## (Intercept) 1466152995 1 12537.771 < 2.2e-16 ***
## sex 9071259 1 77.573 1.113e-13 ***
## species 40705547 2 174.046 < 2.2e-16 ***
## Residuals 10173683 87
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Les résultats montrent qu’en absence d’interaction, les sommes des carrés de type II et de type III sont strictement égales, et par conséquence les F et p values aussi .
Remarque :
L’ajustement du modèle de l’anova à deux facteurs peut également se réaliser avec la fonction lm(), à la place de la fonction aov().
Analyse de la variance à deux facteurs avec des données équilibrées
Il s’agit ici, de vous montrer que le problème du choix des carrés de type II ou III dans un modèle complet, avec interaction, ne se pose pas non plus lorsque les données sont équilibrées, car les résultats sont alors strictement identiques.
Nous employons ici les données mydata_Eq, qui sont parfaitement équilibrées :
table( mydata_Eq$sex,mydata_Eq$species)
##
## Adelie Chinstrap Gentoo
## female 13 13 13
## male 13 13 13 Resultats des carrés de type II
# sans modif des contrastes --> carrés type II
mydata_Eq_type2 <- aov(body_mass_g~sex*species, data=mydata_Eq)
Anova(mydata_Eq_type2)
## Anova Table (Type II tests)
##
## Response: body_mass_g
## Sum Sq Df F value Pr(>F)
## sex 7508205 1 72.8937 1.519e-12 ***
## species 34962756 2 169.7186 < 2.2e-16 ***
## sex:species 1517372 2 7.3657 0.001229 **
## Residuals 7416154 72
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Resultats des carrés de type III
# avec modif des contrastes --> type III
mydata_Eq_type3 <- aov(body_mass_g~sex*species,
contrasts=list(sex=contr.sum,
species=contr.sum),
data=mydata_Eq)
Anova(mydata_Eq_type3, type=3)
## Anova Table (Type III tests)
##
## Response: body_mass_g
## Sum Sq Df F value Pr(>F)
## (Intercept) 1350419263 1 13110.5946 < 2.2e-16 ***
## sex 7508205 1 72.8937 1.519e-12 ***
## species 34962756 2 169.7186 < 2.2e-16 ***
## sex:species 1517372 2 7.3657 0.001229 **
## Residuals 7416154 72
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Nous pouvons observer que les résultats, en termes de sommes des carrés, et ensuite F et p values, sont strictement identiques entre les carrés de type II et ceux de type III.
Nous obtenons encore la même chose avec la fonction AovSum() du package FactoMineR :
#install.packages("FactoMineR")
library(FactoMineR)
mydata_Eq_aovsum <- AovSum(body_mass_g~sex*species, data=mydata_Eq)
mydata_Eq_aovsum$Ftest
## SS df MS F value Pr(>F)
## sex 7508205 1 7508205 72.8937 1.519e-12 ***
## species 34962756 2 17481378 169.7186 < 2.2e-16 ***
## sex:species 1517372 2 758686 7.3657 0.001229 **
## Residuals 7416154 72 103002
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Conclusion
Lorsqu’un modèle d’analyse de la variance comportant un terme d’interaction est ajusté à des données déséquilibrées, il est usuel d’évaluer les effets des facteurs étudiés en calculant des carrés de type III. Ces carrés de type III ne sont pas ceux obtenus par défaut. Trois méthodes permettent d’accéder aux carrés de type III :
modifier les contrastes des facteurs en amont de la modélisation, en utilisant le type
contr.sumet en utilisant la fonctionAnova()avec l’argumenttype =3,ou encore la fonctiondrop1().
- préciser les contrastes
contr.sumpendant la modélisation, par l’argumentcontrastsdes fonctionsaov()oulm(), et en employant ensuite la fonctionAnovaavec l’argumenttype =3, ou la fonctiondrop1() - en utilisant directement la fonction
AovSum()du packageFactoMineR
Cette problématique du calcul des carrés de type III ne se pose pas si :
- le modèle ajusté comporte un terme d’interaction, mais que les données sont équilibrées
- le modèle ajusté ne comporte pas de terme d’interaction (que les données soient équilibrées ou pas).
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. 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
Une réponse
bonjour,
merci pour vos efforts et votre partage:
Si possible un eclaircissement pour le calcul et une interpretation de la somme des carrés de type 3 pour l’intercept