Carrés de type III et carrés de type II dans l'analyse de variance

carrés de type III ou de type II

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)”

(Husson, F., Cornillon, P. A., Guyader, A., Jégou, N., Josse, J., Klutchnikoff, N., … & Thieurmel, B. (2018). R pour la statistique et la science des données. Presses universitaires de Rennes.)

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.sum et en utilisant la fonction Anova() avec l’argument type =3, ou encore la fonction drop1().

  • préciser les contrastes contr.sum pendant la modélisation, par l’argument contrasts des fonctions aov() ou lm(), et en employant ensuite la fonction Anova avec l’argument type =3, ou la fonction drop1()
  • en utilisant directement la fonction AovSum()du package FactoMineR

 

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  

Poursuivez votre lecture

Image par lisa runnels de Pixabay

Laisser un commentaire

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

Bonjour !

vous venez souvent ?

Identifiez-vous pour avoir accès à toutes les fontionnalités !

Aide mémoire off'R ;)

Enregistrez vous pour recevoir gratuitement mes fiches “aide mémoire” (ou cheat sheets) qui vous permettront de réaliser facilement les principales analyses biostatistiques avec le logiciel R et pour être informés des mises à jour du site.