11.2 Modèle sans interactions

La version la plus simple consiste à considérer simplement deux facteurs successivement, c’est-à-dire que la variance est décomposée d’abord selon le premier facteur, et ensuite selon le second.

\[y_{ijk} = \mu + \tau1_j + \tau2_k + \epsilon_i \mathrm{\ avec\ } \epsilon \sim N(0, \sigma)\]

avec \(\tau1_j\) correspondant à l’écart de la moyenne générale \(µ\) à la moyenne selon la j^ème population pour la variable fact1, et \(\tau2_k\) correspondant à l’écart vers le k^ème niveau d’une seconde variable fact2. La formule qui spécifie ce modèle dans R avec les trois variables y, fact1 et fact2 s’écrit :

\[y \sim fact1 + fact2\]

Notez que, quel que soit le niveau considéré pour \(\tau1\), un niveau donné de \(\tau2\) est constant dans l’équation qui décrit ce modèle. Cela signifie que l’on considère que les écarts pour les moyennes de la variable fact2 sont toujours les mêmes depuis les moyennes de fact1. Donc, si une sous-population de fact2 tend à avoir une moyenne, disons, supérieure pour la première sous-population de fact1, elle sera considérée comme ayant les mêmes écarts pour toutes les autres sous-populations de fact1. Evidemment, cette condition n’est pas toujours rencontrée dans la pratique. Le graphique des interactions (Fig. 11.1) permet de visualiser les écarts des moyennes respectives des différentes sous-populations.

read("crabs", package = "MASS", lang = "fr") %>.%
  mutate(., aspect = labelise(
    as.numeric(rear / width),
    "Ratio largeur arrière / max", units = NA)) %>.%
  mutate(., aspect5 = labelise(
    aspect^5,
    "(Ratio largeur arrière /max)^5", units = NA)) %>.%
  select(., species, sex, aspect, aspect5) ->
  crabs2
# Graphique de base pour visualiser les interactions
#chart$base(interaction.plot(crabs2$species, crabs2$sex, crabs2$aspect5))
# Version avec ggplot2
crabs2 %>.%
  group_by(., species, sex) %>.%
  summarise(., aspect5_groups = mean(aspect5)) %>.%
  print(.) %>.% # Tableau des moyennes par groupes
  chart(data = ., aspect5_groups ~ species %col=% sex %group=% sex) +
    geom_line() +
    geom_point()

# `summarise()` has grouped output by 'species'. You can override using the `.groups`
# argument.

# # A tibble: 4 x 3
# # Groups:   species [2]
#   species sex   aspect5_groups
#   <fct>   <fct>          <dbl>
# 1 B       F            0.00727
# 2 B       M            0.00363
# 3 O       F            0.00811
# 4 O       M            0.00427

Figure 11.1: Graphique des interactions entre les variables facteurs (espèce et sexe). Les traits (pratiquement) parallèles indiquent qu’il n’y a pas d’interactions, comme c’est le cas ici.

À vous de jouer !

Au niveau de la description préliminaire des données, nous pourrons utiliser un tableau qui résume la moyenne, l’écart type et le nombre d’observations pour chaque combinaison des deux variables facteurs. Le template de ce code est disponible dans un “snippet” à partir du menu hypothesis tests: means ou .hm, et ensuite two-way ANOVA (description).

crabs2 %>.%
  group_by(., species, sex) %>.%
  summarise(., mean = mean(aspect5), sd = sd(aspect5), count = sum(!is.na(aspect5)))

# `summarise()` has grouped output by 'species'. You can override using the `.groups`
# argument.

# # A tibble: 4 x 5
# # Groups:   species [2]
#   species sex      mean      sd count
#   <fct>   <fct>   <dbl>   <dbl> <int>
# 1 B       F     0.00727 0.00115    50
# 2 B       M     0.00363 0.00124    50
# 3 O       F     0.00811 0.00138    50
# 4 O       M     0.00427 0.00130    50

Pour la visualisation graphique, nous sommes tributaires du nombre d’observations. Avec moins d’une petite dizaine d’observations, nous représenterons des points pour chaque observation et superposerons les moyennes. Lorsque le nombre est plus grand nous pourrons utiliser soit les boites de dispersion, soit le graphique en violon si ce nombre est encore plus grand. Voyons cela sur notre exemple (les “snippets” dans le menu chart: bivariate peuvent être utilisés comme point de départ auquel nous ajoutons la seconde variable facteur pour les facettes multi-graphiques séparée par un |). La Fig. 11.2 montre ce que cela donne si l’on opte pour les boites de dispersion.

chart(data = crabs2, aspect5 ~ species | sex) +
  geom_boxplot()

Figure 11.2: Taille relative de la carapace à l’arrière de crabes L. variegatus (deux variétés et deux sexes), version simple.

La Fig. 11.3 est une version améliorée avec les observations et les moyennes pour chaque sous-groupe ajoutées au graphique selon la même technique que nous avions utilisé pour représenter les données pour l’ANOVA à un facteur.

chart(data = crabs2, aspect5 ~ species | sex) +
  geom_boxplot() +
  geom_jitter(width = 0.05, alpha = 0.5) +
  geom_point(data = group_by(crabs2, species, sex) %>.%
    summarise(., means = mean(aspect5, na.rm = TRUE)),
    f_aes(means ~ species), size = 3, col = "red")

# `summarise()` has grouped output by 'species'. You can override using the `.groups`
# argument.

Figure 11.3: Taille relative de la carapace à l’arrière de crabes L. variegatus (deux variétés et deux sexes), version annotée.

Maintenant que nous avons décrit correctement nos données par rapport à l’analyse que nous souhaitons faire, nous pouvons réaliser notre ANOVA à deux facteurs. Nous devons vérifier l’homoscédasticité, mais le test de Bartlett que nous réalisons revient au même que celui que nous avons fait en décomposant toutes les sous-populations. Comme nous n’avons pas nécessairement ce calcul réalisé (la variable group que nous avions calculée au module 10), nous utilisons la fonction interaction() qui effectue ce calcul pour nous directement dans la formule :

bartlett.test(data = crabs2, aspect5 ~ interaction(species, sex))

# 
#   Bartlett test of homogeneity of variances
# 
# data:  aspect5 by interaction(species, sex)
# Bartlett's K-squared = 1.7948, df = 3, p-value = 0.6161

À vous de jouer !

Si vous comparez avec le test que nous avions fait dans le cas de l’ANOVA à un facteur sur la variable group, vous constaterez qu’il donne exactement le même résultat. Nous continuons avec notre ANOVA. Nous avons un “snippet” pour cela dans le menu hypothesis tests: means à partir de .hm qui se nomme two-way ANOVA (without interactions).

anova(anova. <- lm(data = crabs2, aspect5 ~ species + sex))

# Analysis of Variance Table
# 
# Response: aspect5
#            Df     Sum Sq    Mean Sq F value    Pr(>F)    
# species     1 0.00002753 0.00002753   17.15  5.12e-05 ***
# sex         1 0.00069935 0.00069935  435.66 < 2.2e-16 ***
# Residuals 197 0.00031624 0.00000161                      
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comme la variance est décomposée en trois étapes (selon l’espèce, puis selon le sexe, puis les résidus), nous avons trois lignes dans le tableau de l’ANOVA. Nous effectuons deux tests. Le premier consiste à comparer les carrés moyens (Mean Sq) pour l’espèce par rapport aux résidus. Donc, la valeur F est le ratio de la somme des carrés species divisé par la somme des carrés des résidus, et cette valeur est reportée sur la loi de distribution théorique F pour obtenir une première valeur P (ici 5,12 . 10^-5). De même, le second test qui s’intéresse au sexe calcule la valeur F via la ratio de la somme des carrés pour sex divisé par la somme des carrés des résidus, et la loi F nous permet de calculer une seconde valeur P (ici 2,2 . 10^-16).

Nous interprétons chacun des deux tests séparément. Dans notre cas, nous pouvons dire avec un seuil \(\alpha\) de 5% que nous rejetons \(H_0\) dans les deux cas. Donc le rapport largeur arrière sur largeur max de la carapace est significativement différent au seuil \(\alpha\) de 5% à la fois en fonction de l’espèce (F = 17,15, ddl = 197 et 1, valeur P << 0,001) et du sexe (F = 436, ddl = 197 et 1, valeur P << 0,001)⁵¹.

La suite logique consiste à réaliser des tests “post hoc”. Ils ne sont pas vraiment nécessaires ici puisque nous n’avons que deux niveaux pour chacune des deux variables, mais nous les réalisons quand même pour montrer le code correspondant. Un template est accessible via le “snippet” anova - multiple comparisons [multcomp] du menu .hm. Pensez juste à rajouter le second facteur sex dans les arguments de la fonction mcp().

summary(anovaComp. <- confint(multcomp::glht(anova.,
  linfct = multcomp::mcp(species = "Tukey", sex = "Tukey")))) # Add a second factor if you want

# 
#    Simultaneous Tests for General Linear Hypotheses
# 
# Multiple Comparisons of Means: Tukey Contrasts
# 
# 
# Fit: lm(formula = aspect5 ~ species + sex, data = crabs2)
# 
# Linear Hypotheses:
#                       Estimate Std. Error t value Pr(>|t|)    
# species: O - B == 0  0.0007420  0.0001792   4.141 0.000102 ***
# sex: M - F == 0     -0.0037399  0.0001792 -20.872  < 1e-10 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# (Adjusted p values reported -- single-step method)

.oma <- par(oma = c(0, 5.1, 0, 0)); plot(anovaComp.); par(.oma); rm(.oma)

Ceci confirme que les différences sont significatives au seuil \(\alpha\) de 5%. Il ne nous reste plus qu’à vérifier la distribution des résidus de l’ANOVA pour que notre analyse soit complète (Fig. 11.4).

Figure 11.4: Graphique quantile-quantile des résidus pour l’ANOVA à deux facteurs sans interactions de aspect^5.

# [1]  65 103

Encore une fois, nous voyons que les résidus sont quasiment les mêmes que précédemment, mais cela n’aurait pas été le cas si une interaction existait. La distribution s’éloigne un peu d’une Gaussienne pour les valeurs élevées surtout. Mais comme l’ANOVA est robuste à ce critère, et que l’homoscédasticité a été vérifiée sur la transformation puissance 5 de notre variable, nous pouvons conserver notre analyse moyennant une précaution supplémentaire : vérifier que les valeurs P sont beaucoup plus petites que notre seuil comme sécurité supplémentaires contre les approximations liées à la légère violation de la contrainte de distribution normale des résidus. C’est le cas ici, et nous pouvons donc conclure notre analyse.

Conditions d’application

• échantillon représentatif (par exemple, aléatoire),
• observations indépendantes,
• variable réponse quantitative,
• deux variables explicatives qualitatives à deux niveaux ou plus,
• distribution normale des résidus \(\epsilon_i\),
• homoscédasticité (même variance intragroupes),
• pas d’interactions entre les deux variables explicatives.