Moyennes marginales
Très récemment, je me suis interrogée sur l’estimation des moyennes marginales, principalement dans le contexte de l’ANOVA à 2 facteurs et celui de l’analyse de covariance.
Comme je sais que je ne suis pas la seule à me poser des questions au sujet de ces moyennes marginales, j’ai écrit cet article pour partager ce que j’ai compris.
Table des matières
Contexte
Je me suis intéressée à l’estimation des moyennes marginales parce que lorsque je réalise des ANOVA à deux facteurs, et que je souhaite ensuite comparer les moyennes des modalités de l’un des facteurs, j’utilise, à présent, la fonction emmeans() (du package emmeans) qui emploie des moyennes marginales pour réaliser ces comparaisons.
De même, lorsque j’utilise une ANCOVA, et que je souhaite ensuite réaliser des tests post-hoc pour comparer les moyennes des différents groupes, j’emploie encore la fonction emmeans() du package emmeans qui a toujours recours à des moyennes marginales.
En outre, très récemment, lorsque j’ai exploré le package modelbased (qui appartient au meta package easystats, j’ai découvert la fonction estimate_means() qui permet d’estimer les moyennes marginales à partir de nombreux modèles de régression.
Il était donc temps de se pencher d’un peu plus près sur cette notion.
Qu’est ce qu’une moyenne marginale ?
Dans la page d’aide de la fonction emmeans() il est dit :
Estimated marginal means or EMMs (sometimes called least-squares means) are predictions from a linear model over a reference grid; or marginal averages thereof.
Et dans le package modelbased :
Marginal means are basically means extracted from a statistical model, and represent average of response variable for each level of predictor variable
À ce stade, je retiens que les moyennes marginales reposent sur des modèles et des prédictions de ces modèles, mais ce n’est pas super limpide !
Je suis donc passé à des choses concrètes en essayant sur une ANOVA à 2 facteurs et une ANCOVA.
Moyennes marginales dans le cadre de l’ANOVA à 2 facteurs
Data pigs
Pour illustrer cet exemple, j’ai choisi d’employer les données “pigs” contenues dans le package emmeans. Ces données correspondent à un plan factoriel 2*2 avec des effectifs déséquilibrés.
Elles comportent 29 observations et 3 variables :
- source : Source de protéines dans l’alimentation (facteur à 3 niveaux : farine de poisson, farine de soja, lait écrémé en poudre) ; c’est le premier facteur
- percent: Pourcentage de protéines dans l’alimentation (numérique avec 4 valeurs : 9, 12, 15 et 18); c’est le second facteur
- conc : Concentration de leucine plasmatique libre, en mcg/ml; c’est la réponse.
library(emmeans)
# passage de la variable percent en factor
pigs$percent <- as.factor(pigs$percent)
head(pigs)
## source percent conc
## 1 fish 9 27.8
## 2 fish 9 23.7
## 3 fish 12 31.5
## 4 fish 12 28.5
## 5 fish 12 32.8
## 6 fish 15 34.0
# nombre d'observation par croisement des modalités
table(pigs$percent, pigs$source)
##
## fish soy skim
## 9 2 3 3
## 12 3 3 3
## 15 2 3 2
## 18 3 1 1 Les données sont bien déséquilibrées.
Voici une visualisation des moyennes observées :
Ajustement du modèle ANOVA à 2 facteurs avec interaction
pigs.lm1 <- lm(conc~source*percent,
contrasts=list(source=contr.sum, percent=contr.sum),
data=pigs)
library(car)
Anova(pigs.lm1, type=3)
## Anova Table (Type III tests)
##
## Response: conc
## Sum Sq Df F value Pr(>F)
## (Intercept) 37176 1 1671.5239 < 2.2e-16 ***
## source 1253 2 28.1801 4.003e-06 ***
## percent 508 3 7.6193 0.00193 **
## source:percent 215 6 1.6074 0.20547
## Residuals 378 17
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Test post hoc pour la variable percent
Si je veux, subséquemment à l’ANOVA, comparer les moyennes relatives aux différents pourcentages (ici je considère que j’ai le droit de le faire puisque l’interaction est non significative), alors j’emploie la fonction emmeans(), comme ceci :
emmeans(pigs.lm1, specs=pairwise~percent)
## NOTE: Results may be misleading due to involvement in interactions
## $emmeans
## percent emmean SE df lower.CL upper.CL
## 9 31.9 1.70 17 28.3 35.5
## 12 38.0 1.57 17 34.7 41.3
## 15 40.3 1.82 17 36.4 44.1
## 18 45.0 2.40 17 39.9 50.1
##
## Results are averaged over the levels of: source
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## percent9 - percent12 -6.08 2.31 17 -2.629 0.0754
## percent9 - percent15 -8.35 2.49 17 -3.359 0.0177
## percent9 - percent18 -13.08 2.94 17 -4.449 0.0018
## percent12 - percent15 -2.27 2.40 17 -0.944 0.7820
## percent12 - percent18 -7.00 2.87 17 -2.439 0.1071
## percent15 - percent18 -4.73 3.01 17 -1.572 0.4193
##
## Results are averaged over the levels of: source
## P value adjustment: tukey method for comparing a family of 4 estimates Les moyennes marginales sont rapportées dans la première partie des résultats. Elles sont égales à :
- 31.9 pour la modalité 9 percent
- 38.0 pour la modalité 12 percent
- 40.3 pour la modalité 15 percent
- 45.0 pour la modalité 18 percent
Dans la seconde partie, on peut voir les comparaisons de ces moyennes deux à deux.
Intuitivement, on comprend que ces moyennes sont des niveaux moyens de concentration, en prenant en considération toutes les sources, mais comment est fait le calcul ?
En première intention, on pourrait penser que les moyennes marginales sont simplement les moyennes de l’ensemble des données relatives à une modalité (9, 12, 15 et 18). Par exemple que 45.0 est la moyenne de toutes les données observées avec la condition pourcentage = 18. Vérifions :
library(tidyverse)
pigs %>%
filter(percent=="18") %>%
summarise(moy=mean(conc))
## moy
## 1 39.94 Et ben non, ce n’est pas ça !
En seconde intention (et parce qu’on est bien forcé de réfléchir à présent, 😉) on peut faire l’hypothèse que les moyennes marginales sont les moyennes des moyennes observées pour chaque source. Cela à d’ailleurs plus de sens avec le terme “marginal”.
Et c’est bien ça !
On peut le recalculer à la main avec R, comme ça :
pigs %>%
group_by(percent, source) %>%
summarise(moy=mean(conc)) %>%
ungroup() %>%
group_by(percent) %>%
summarise(mm=mean(moy))
## `summarise()` has grouped output by 'percent'. You can override using the
## `.groups` argument.
## # A tibble: 4 × 2
## percent mm
## <fct> <dbl>
## 1 9 31.9
## 2 12 38.0
## 3 15 40.3
## 4 18 45.0 Remarque: si seules les estimations des moyennes marginales vous intéresse, vous pouvez plus facilement employer la fonction estimate_means() du package modelbased, comme ceci :
library(modelbased)
estimate_means(pigs.lm1, at="percent")
## NOTE: Results may be misleading due to involvement in interactions
## Estimated Marginal Means
##
## percent | Mean | SE | 95% CI
## ---------------------------------------
## 9 | 31.93 | 1.70 | [28.35, 35.51]
## 12 | 38.01 | 1.57 | [34.69, 41.33]
## 15 | 40.28 | 1.82 | [36.45, 44.11]
## 18 | 45.01 | 2.40 | [39.94, 50.08]
##
## Marginal means estimated at percent Bon à ce stade, je comprends comment ça marche pour l’ANOVA à 2 facteurs, mais le lien avec les modèles et les prédictions est encore un peu flou. Mais ça va s’éclaircir avec l’ANCOVA !
Et là, je vous rapporte l’exemple que j’ai trouvé dans la vignette du package modelbased parce que je le trouve très bon !
Moyennes marginales et ANCOVA
Commençons par un modèle ANOVA à 1 facteur
Nous nous intéressons ici à la largeur des sépales en fonction de l’espèce. Nous pouvons calculer les moyennes observées comme ceci :
iris %>%
group_by(Species) %>%
summarise(moy=mean(Sepal.Width))
## # A tibble: 3 × 2
## Species moy
## <fct> <dbl>
## 1 setosa 3.43
## 2 versicolor 2.77
## 3 virginica 2.97 Le modèle ANOVA nous permet de de modéliser la largeur du sépale en fonction de l’espèce :
iris.lm1 <- lm(Sepal.Width~Species, data=iris)
Anova(iris.lm1)
## Anova Table (Type II tests)
##
## Response: Sepal.Width
## Sum Sq Df F value Pr(>F)
## Species 11.345 2 49.16 < 2.2e-16 ***
## Residuals 16.962 147
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Ici, comme vous pouvez le voir, les moyennes marginales estimées à partir de ce modèle sont simplement les moyennes observées.
means <- estimate_means(iris.lm1)
## We selected `at = c("Species")`.
means
## Estimated Marginal Means
##
## Species | Mean | SE | 95% CI
## ---------------------------------------
## setosa | 3.43 | 0.05 | [3.33, 3.52]
## versicolor | 2.77 | 0.05 | [2.68, 2.86]
## virginica | 2.97 | 0.05 | [2.88, 3.07]
##
## Marginal means estimated at Species Nous pouvons les représenter comme ceci :
library(see)
ggplot(iris, aes(x = Species, y = Sepal.Width, fill = Species)) +
geom_violin() +
geom_jitter2(width = 0.05, alpha = 0.5) +
geom_line(data = means, aes(y = Mean, group = 1), size = 1) +
geom_pointrange(
data = means,
aes(y = Mean, ymin = CI_low, ymax = CI_high),
size = 1,
color = "white"
) +
theme_modern() 
Modele ancova
Comme je l’expliquais dans cet article, l’ANCOVA peut être vue sous l’angle d’une comparaison de moyennes ajustées sur une co-variable numérique.
Ici, il s’agit ici d’estimer les moyennes de la largeur du sépale (Sepal.Width) des 3 espèces d’iris, ajustées sur la valeur de la largeur du pétale (Petal.Width); autrement dit en prenant en compte l’effet de la largeur du Sepal.
Pour cela nous pouvons ajuster le modèle ANCOVA suivant :
iris.lm2<- lm(Sepal.Width ~ Species * Petal.Width, data = iris)
Anova(iris.lm2)
## Anova Table (Type II tests)
##
## Response: Sepal.Width
## Sum Sq Df F value Pr(>F)
## Species 11.3059 2 62.7151 < 2.2e-16 ***
## Petal.Width 3.7554 1 41.6638 1.551e-09 ***
## Species:Petal.Width 0.2269 2 1.2585 0.2872
## Residuals 12.9797 144
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Si je calcule les moyennes marginales d’après ce modèle ANCOVA, je vais voir qu’elles ne seront pas les mêmes que précédemment, puisqu’il s’agit ici des moyennes ajustées sur la largeur des pétales.
means_complex <- estimate_means(iris.lm2, at="Species")
means_complex
means_complex
## Estimated Marginal Means
##
## Species | Mean | SE | 95% CI
## ---------------------------------------
## setosa | 4.23 | 0.39 | [3.45, 5.00]
## versicolor | 2.64 | 0.05 | [2.54, 2.74]
## virginica | 2.45 | 0.14 | [2.18, 2.72]
##
## Marginal means estimated at Species Nous pouvons visualiser l’évolution des moyennes marginales, entre le modèle ANOVA à un facteur, et le modèle ANCOVA, comme ceci :
ggplot(iris, aes(x = Species, y = Sepal.Width, fill = Species)) +
geom_violin() +
geom_jitter2(width = 0.05, alpha = 0.5) +
geom_line(data = means, aes(y = Mean, group = 1), size = 1, alpha = 0.25) +
geom_pointrange(
data = means,
aes(y = Mean, ymin = CI_low, ymax = CI_high),
size = 1,
color = "white"
) +
geom_line(data = means_complex, aes(y = Mean, group = 1), size = 1) +
geom_pointrange(
data = means_complex,
aes(y = Mean, ymin = CI_low, ymax = CI_high),
size = 1,
color = "orange"
) +
theme_modern() 
En blanc les moyennes marginales estimées issues de l’ANOVA à 1 facteur (qui sont aussi les moyennes observées), et en jaune les moyennes marginales estimées issues de l’ANCOVA, (qui sont des moyennes ajustées sur la co-variable Petal.Width).
Si j’emploie, à présent, la fonction emmeans() pour réaliser mes tests post-hoc, ceux-ci vont bien être basés sur les moyennes marginales (qui sont, bien évidement, identiques à celles estimées par la fonction estimate_means()) :
emmeans(iris.lm2, specs=pairwise~Species)
## NOTE: Results may be misleading due to involvement in interactions
## $emmeans
## Species emmean SE df lower.CL upper.CL
## setosa 4.23 0.3903 144 3.45 5.00
## versicolor 2.64 0.0506 144 2.54 2.74
## virginica 2.45 0.1359 144 2.18 2.72
##
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## setosa - versicolor 1.590 0.394 144 4.039 0.0003
## setosa - virginica 1.774 0.413 144 4.293 0.0001
## versicolor - virginica 0.184 0.145 144 1.272 0.4131
##
## P value adjustment: tukey method for comparing a family of 3 estimates Le plot précédent permet de mieux comprendre la différence entre une moyenne marginale et une moyenne observée, mais concrètement, comment sont-elles caculées ici ?
J’ai trouvé un peu d’info dans la vignette “Basics of estimated marginal means “du package emmeans.
Voici la traduction :
“Les moyennes marginales estimées sont basées sur un modèle et non directement sur des données. Leur base est ce que nous appelons la grille de référence pour un modèle donné. Pour obtenir la grille de référence, considérez tous les prédicteurs du modèle. Voici les règles par défaut pour construire la grille de référence:
- Pour chaque prédicteur qui est un facteur, utilisez ses niveaux (supprimez ceux qui ne sont pas utilisés)
- Pour chaque prédicteur numérique (covariable), utilisez sa moyenne.
Donc d’après ce que je comprends, pour calculer les moyennes marginales dans une ANCOVA, des prédictions sont faites en fixant la covariable à son niveau moyen, puis les prédictions sont moyennées.
On essaye !
Je crée un dataset “iris2” contenant les variables Sepal.Width, Species et Petal.Width, mais je fixe Petal Width à sa moyenne( 1.199)
new_iris <- iris %>%
select(Sepal.Width,Species,Petal.Width) %>%
mutate(Petal.Width=mean(Petal.Width))
head(new_iris)
## Sepal.Width Species Petal.Width
## 1 3.5 setosa 1.199333
## 2 3.0 setosa 1.199333
## 3 3.2 setosa 1.199333
## 4 3.1 setosa 1.199333
## 5 3.6 setosa 1.199333
## 6 3.9 setosa 1.199333 Ensuite, je fais des prédictions de Sepal.Width pour chacun de 150 fleurs du fichier, en considérant son espèce et que la valeur de Petal Width égale 1.199.
new_iris$predSW <- predict(iris.lm2, new_iris)
head(new_iris)
## Sepal.Width Species Petal.Width predSW
## 1 3.5 setosa 1.199333 4.226123
## 2 3.0 setosa 1.199333 4.226123
## 3 3.2 setosa 1.199333 4.226123
## 4 3.1 setosa 1.199333 4.226123
## 5 3.6 setosa 1.199333 4.226123
## 6 3.9 setosa 1.199333 4.226123 Puis je calcule les moyennes de ces prédictions de sepal.Width, pour chaque espèce :
new_iris %>%
group_by(Species) %>%
summarise(moy=mean(predSW))
## # A tibble: 3 × 2
## Species moy
## <fct> <dbl>
## 1 setosa 4.23
## 2 versicolor 2.64
## 3 virginica 2.45 Et c’est bien ça, je retombe sur les moyennes marginales !
Si on résume un peu les choses :
- les moyennes marginales sont estimées en réalisant des prédictions selon une grille de valeurs pour la co-variable ou le second facteur, puis ces prédictions sont moyennées
- dans le cadre de l’ANOVA à un facteur, les moyennes marginales sont égales aux moyennes observées puisque les prédictions de l’ANOVA à un facteur sont les moyennes observées (la prédiction de largeur de sépale pour l’espèce setosa, c’est la moyenne des largeurs de fleurs de l’espèce setosa)
- dans le cadre de l’ANOVA à 2 facteurs avec interaction, j’ai pu aussi retrouver les moyennes marginales sans passer par les prédictions du modèle, mais là encore parce que les prédictions sont égales aux moyennes observées (la prédiction de la source fish pour percent = 9, est la moyenne observée pour cette condition)
J’espère que cet article vous permettra, à vous aussi, de mieux appréhender cette notion de moyennes marginales, et de bien comprendre comment elles sont calculées.
Vous souhaitez soutenir mon travail ?
Si vous souhaitez soutenir mon travail, vous pouvez faire un don libre sur la page tipeee du blog :
2 réponses
Bonjour
Merveilleuse analyse et modélisation.
C’est super.
Mais est ce possible d’avoir le PDF ?
Bonjour,
Merci pour vote retour sur cet article.
Celui-ci est uniquement consultable, il n’est pas possible d’obtenir le pdf;
Bonne continuation