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.
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.
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.
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 :
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 :
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 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 à :
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 !
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() 
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:
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 :
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.
Si vous souhaitez soutenir mon travail, vous pouvez faire un don libre sur la page tipeee du blog :
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.
4 Responses
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
Bonjour,
Je pensais être la seule à me casser la tête à trouver les valeurs estimées de emmeans. Je n’arrive toujours pas à les trouver. J’ai une anova à deux facteurs sans interaction. Donc je n’ai pas de répétitions par couple de modalité. Mes données sont déséquilibrées( Il y a des couples de modalités où je n’aurai pas d’observations). Si on s’interesse à percent par exemple. Alors la moyenne marginale dans ce cas (sans repetition et des couples d’effectifs 0) est la moyenne des valeurs (de source) par modalité de percent. Je ne tombe pas sur les résultats de emmeans. J’ai essayé aussi avec la fonction predict en suivant les étapes de la vignette du package modelbased. Toujours pas les mêmes résultats entre le predict et les moyennes estimées de emmeans. Si vous avez une idée…. Je vous remercie bien
Bonne journée
En fait ça marche avec predict. Il a dû prendre un résultat d’avant. Merci en tout cas pour vos cours, toujours très clair. Bonne journée