ANOVA à un facteur : partie 2 – la pratique
Table des matières
Pour faire suite au précédent article sur l’ANOVA à un facteur qui était plutôt théorique, je vais vous montrer, dans cet article comment réaliser, en pratique, une ANOVA à un facteur avec R.
Pour cela, je vais utiliser le jeu de données cholesterol du package multcomp. Si vous voulez reproduire les exemples, il faudra donc importer et charger ce package.
Le jeu de données “cholestérol” provient d’une étude clinique visant a évaluer l’effet de cinq traitements sur la réduction du cholesterol.
La modalité “1time” correspond à 20 mg d’une molécule une fois par jour, la modalité “2times correspond” à 10 mg deux fois par jour, la modalité “4times” correspond à 5mg quatre fois par jour, les modalités “drugD” et “drugE” correspondent à l’utilisation de 2 molécules qui servent de groupe “contrôle”.
Pour plus d’information sur ce jeu de données, utilisez la syntaxe?cholesterol.
library(multcomp)
data(cholesterol)
head(cholesterol)
## trt response
## 1 1time 3.8612
## 2 1time 10.3868
## 3 1time 5.9059
## 4 1time 3.0609
## 5 1time 7.7204
## 6 1time 2.7139 Exploration visuelle des données
En premier lieu, il est toujours utile de représenter les données pour se faire une première idée. Pour cette étape, j’ai l’habitude d’utiliser des boxplots, en ajoutant par-dessus les points. Cela permet de se rendre compte du nombre de données par modalité, ainsi que de leur distribution (plutôt normal, plutôt asymétrique ?), et de la présence éventuelle d’outliers (valeurs extrêmes). Ces deux derniers points pouvant biaiser les résultats de l’analyse.
library(ggplot2)
ggplot(cholesterol, aes(y=response, x=trt,colour=trt ,fill=trt))+
geom_boxplot(alpha=0.5, outlier.alpha=0)+
geom_jitter(width=0.25)+
theme_classic() 
L’argument outlier.alpha=0 dans la fonction geom_boxplot permet de ne pas représenter deux fois un point outlier (une fois avec la fonction geom_boxplot, un fois avec la fonction geom_jitter). La fonction geom_jitter permet de représenter les points sans qu’ils se chevauchent. L’argument width=0.25 permet de gerer l’écartement des points. Ici, par exemple, les données de la modalité “2times” sont un peu asymétriques, avec la présence de plusieurs outliers (points dépassant les traits verticaux).
Calculer les moyennes et leurs intervalles de confiance
Il peut également être intéressant de calculer les moyennes et leur intervalle de confiance. Une façon rapide de le faire est d’employer la fonction ci.mean du package Publish. Les estimations des intervalles de confiances sont basées sur une distribution t.
library(Publish)
ci.mean(response~trt,data=cholesterol)
## trt mean CI-95%
## 1time 5.78 [3.72;7.84]
## 2times 9.22 [6.73;11.72]
## 4times 12.37 [10.28;14.47]
## drugD 15.36 [12.89;17.83]
## drugE 20.95 [18.55;23.34] Il est également possible d’estimer les intervalles de confiance par approche bootstrap, à l’aide du package slipper, comme expliqué dans ici, mais cela prend peu plus de temps.
Réalisation de l'ANOVA
Sous R, il y a deux commandes principales pour réaliser l’ANOVA : aov et lm. J’utilise toujours lm, mais c’est surtout une question d’habitude.
lm1 <- lm(response~trt, data=cholesterol)
aov1 <- aov(response~trt, data=cholesterol) Visualisation des résultats
Les résultats de l’ANOVA sont généralement représentés sous la forme d’une table de variance.
Lorsque l’ANOVA a été réalisée avec la fonction lm, je vous recommande d’utiliser la fonction Anova du package car, plutôt que la commande anova du package de base.
Dans le cadre d’une ANOVA à un facteur, la commande anova n’est pas problématique, mais elle le devient lorsqu’il s’agit d’une ANOVA à 2 facteurs, car elle fournit des résultats séquentiels, qui ont pour conséquence de dépendre de l’ordre d’apparition des facteurs dans la fonction lm, et qui ne sont pas, au final, ceux de la table ANOVA attendue.
library(car)
Anova(lm1)
## Anova Table (Type II tests)
##
## Response: response
## Sum Sq Df F value Pr(>F)
## trt 1351.37 4 32.433 9.819e-13 ***
## Residuals 468.75 45
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Si vous utilisez la fonction aov(), vous pouvez simplement utiliser la fonction summary() pour visualiser la table de résultat:
summary(aov1)
## Df Sum Sq Mean Sq F value Pr(>F)
## trt 4 1351.4 337.8 32.43 9.82e-13 ***
## Residuals 45 468.8 10.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Vérification des hypothèses
Comme évoqué dans le premier article sur l’ANOVA à un facteur le résultat de l’ANOVA à un facteur est valide (on peut avoir confiance dans ce résultat), si 3 hypothèses sont vérifiées :
- Les résidus sont indépendants :
- Les résidus suivent une loi normale de moyenne 0
- Les résidus relatifs aux différentes modalités sont homogènes (ils ont globalement la même dispersion), autrement dit leur variance est constante.
Indépendance des résidus
L’indépendance des résidus signifie que les résidus ne doivent pas être corrélés entre eux. Par exemple, il ne doit pas avoir de lien entre un résidu et celui de la donnée suivante, ou précédente. On voit cela, lorsque des données sont répétées sur des sujets identiques. On parle alors d’autocorrélation des résidus. De la même façon, les résidus ne doivent pas être corrélés au facteur étudié.
L’absence d’autocorrélation se valide par l’étude du plan expérimental : pour réaliser une ANOVA à un facteur, il ne doit pas avoir de données répétées. Si c’est le cas, il faut utiliser un autre type de modèle, comme un modèle linéaire à effet mixte.
Il est aussi possible de faire un test statistique, le test de Dubin-Watson. L’hypothèse nulle de ce test spécifie un coefficient d’autocorrélation = 0, alors que l’hypothèse
alternative spécifie un coefficient d’autocorrélation différent de 0 . On conclut donc à l’absence d’autocorrélation lorsque la pvalue du test est supérieure à 0.05.
library(car)
durbinWatsonTest(lm1)
## lag Autocorrelation D-W Statistic p-value
## 1 -0.1265478 2.216953 0.824
## Alternative hypothesis: rho != 0 Ici la p-value est < 0.05, l’hypothèse H0 n’est donc pas rejetée, et on conclut à l’absence d’auto-corrélation.
Normalité des résidus
Pour vérifier cette hypothèse, on utilise généralement un QQplot et/ou un test de normalité comme le test de Shapiro-Wilk.
plot(lm1,2) 
Ici les points sont bien répartis le long de la ligne, cela signifie que les résidus sont distribués selon une loi normale. Le fait que les points soient centrés sur 0 (sur l’axe des y), montre que leur moyenne est égale à 0.
L’hypothèse nulle du test de normalité de Shapiro Wilk spécifie que les résidus suivent une loi normale, alors que son hypothèse alternative spécifie qu’ils suivent une autre distribution quelconque. Pour accepter la normalité des résidus, il est donc nécessaire d’obtenir une p-value > 0.05.
shapiro.test(residuals(lm1))
##
## Shapiro-Wilk normality test
##
## data: residuals(lm1)
## W = 0.98864, p-value = 0.9094 Homogénéité des variances
L’hypothèse d’homogénéité des variances, c’est-à-dire l’hypothèse que les résidus ont une variance constante, peut s’évaluer graphiquement et/ou à l’aide d’un test statistique.
La méthode graphique consiste à représenter les résidus standardisés en fonction des valeurs prédites (les moyennes des différents traitements).
plot(lm1,3) 
On voit ici que les dispersions des résidus (leurs écartements verticaux) relatives à chaque modalité de traitement sont globalement identiques, l’hypothèse d’homogénéité des résidus est acceptée.
On peut également utiliser le test de Bartlett, le test de Levene, ou encore le test de Fligner-Killeen. Leurs hypothèses nulles spécifient que les variances des différents groupes sont globalement identiques. A l’inverse, leurs hypothèses alternatives spécifient qu’au moins 2 variances (les variances de 2 modalités) sont différentes.
Pour accepter l’hypothèse d’homogénéité des résidus, il est donc nécessaire d’obtenir une pvalue >0.05.
bartlett.test(residuals(lm1)~cholesterol$trt)
##
## Bartlett test of homogeneity of variances
##
## data: residuals(lm1) by cholesterol$trt
## Bartlett's K-squared = 0.57975, df = 4, p-value = 0.9653
library(car)
leveneTest(residuals(lm1)~cholesterol$trt)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 4 0.0755 0.9893
## 45
fligner.test(residuals(lm1)~cholesterol$trt)
##
## Fligner-Killeen test of homogeneity of variances
##
## data: residuals(lm1) by cholesterol$trt
## Fligner-Killeen:med chi-squared = 0.74273, df = 4, p-value = 0.946 Ici, la p-avlue est largement supérieure à 0.05, l’hypothèse d’homogénéité des résidus est donc acceptée.
Astuce
Pour visualiser tous les plots du diagnostic de régression en une fois, il est possible d’utiliser les commande suivantes :
par(mfrow=c(2,2)) # permet de séparer la fenêtre graphique en 4 parties (2 lignes, et 2 colonnes)
plot(lm1) 
Pour repasser à une fenêtre graphique sans division, utilisez par(mfrow=c(1,1))
Interprétation des résultats
Les hypothèses étant validées, les résultats peuvent être interprétés.
Anova(lm1)
## Anova Table (Type II tests)
##
## Response: response
## Sum Sq Df F value Pr(>F)
## trt 1351.37 4 32.433 9.819e-13 ***
## Residuals 468.75 45
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 La pvalue de l’ANOVA est < 0.05, ce résultat indique alors que les niveaux de réduction du cholestérol des 5 traitements sont globalement différents. Se pose ensuite évidemment la question de l’identification des moyennes différentes entre elles, ou différentes du contrôle. Pour cela on parle de tests post-hoc ou de comparaisons multiples.
Comparaisons multiples
Le problème de l'augmentation du risque alpha global
Les comparaisons multiples consistent à comparer deux à deux plusieurs moyennes, voire toutes les moyennes entre elles. Le problème de ce type de procédure est que la multiplication des tests de comparaisons entraîne une augmentation du risque de se tromper en mettant en évidence des différences significatives qui ne sont dues qu’au hasard. D’un point de vue statistique, on dit que le risque alpha global augmente. En général, le risque alpha, pour un seul test, est fixé à 5% . Lorsque k tests sont réalisés, le risque alpha global devient:
$ \alpha_{global} = 1- (1- \alpha)^k $
Dans cette situation, il est nécessaire d’utiliser une procédure d’ajustement des p-values obtenues lors des comparaisons pour que les résultats soient en quelque sorte recalibrés pour un risque alpha global de 5%. Les principales procédures sont celles de “holm”, “bonferroni”, etc..
D’un point de vu pratique, le package multcomp permet de réaliser facilement ces comparaisons multiples en appliquant automatiquement une procédure d’ajustement.
Comparaisons au témoin : le test de Dunnett
Dans le jeu de données cholestrol, d’après la page d’aide (?cholesterol), il existe deux modalités contrôle, les modalités “drugD”” et “drugE”. Imaginons que le seul groupe contrôle soit drugE, par exemple.
Pour comparer les moyennes de tous les autres groupes à la moyenne du contrôle, on utilise un test de Dunnett. Pour commencer, il faut définir la modalité “drugE” comme modalité de référence. Pour cela, on utilise la fonction relevel
Par défaut, la modalité de référence, est la première dans l’ordre alphabétique. On peut la visualiser comme cela :
levels(cholesterol$trt)
## [1] "1time" "2times" "4times" "drugD" "drugE" Remarque : pour plus de détails sur la manipulation des variables catégorielles, vous pouvez consulter cet article.
Ici, la modalité, ou groupe, de référence est “1time”. Pour passer la remplacer par “drugE”, on utilise cette commande :
cholesterol$trt <- relevel(cholesterol$trt, ref="drugE")
# verification
levels(cholesterol$trt)
## [1] "drugE" "1time" "2times" "4times" "drugD" La modalité drugE est bien passé en première place.
Ensuite, il est nécessaire de refaire tourner l’ANOVA à un facteur :
lm2 <- lm(response~trt, data=cholesterol) Evidemment, le résultat de l’ANOVA n’est pas modifié :
Anova(lm2)
## Anova Table (Type II tests)
##
## Response: response
## Sum Sq Df F value Pr(>F)
## trt 1351.37 4 32.433 9.819e-13 ***
## Residuals 468.75 45
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Pour réaliser toutes les comparaisons deux à deux avec le groupe contôle “drugE”, il suffit d’utiliser les commandes suivantes :
mc_dunnett <- glht(lm2, linfct=mcp(trt="Dunnett"))
summary(mc_dunnett)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Dunnett Contrasts
##
##
## Fit: lm(formula = response ~ trt, data = cholesterol)
##
## Linear Hypotheses: ## Estimate Std. Error t value Pr(>|t|)
## 1time - drugE == 0 -15.166 1.443 -10.507 < 1e-04 ***
## 2times - drugE == 0 -11.723 1.443 -8.122 < 1e-04 ***
## 4times - drugE == 0 -8.573 1.443 -5.939 < 1e-04 ***
## drugD - drugE == 0 -5.586 1.443 -3.870 0.00131 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method) Chaque ligne correspond à une comparaison. Ici, toutes les moyennes sont différentes du groupe contrôle “drugE”. L’avantage de cette procédure est qu’elle fournit, pour chaque comparaison, une p-value ajustée qui va prendre en compte le nombre de comparaisons pour conserver, au final, un risque alpha global de 5%.
On peut également obtenir une représentation graphique des comparaisons:
plot(mc_dunnett) 
Si les étiquettes des comparaisons sont tronquées, il suffit de modifier les marges , comme ceci :
par(mar=c(3,7,3,3))
plot(mc_dunnett) 
NB : Par défaut, les tests réalisés sont des tests bilatéraux, mais cela peut être changé au sein de la fonction glht, ,à l’aide de l’argument alternative = "less" par exemple
Toutes les comparaisons deux à deux : le test de Tukey
Parfois, ce qu’on souhaite, c’est comparer toutes les moyennes deux à deux. Pour cela, on utilise le test de Tukey :
mc_tukey <- glht(lm2, linfct=mcp(trt="Tukey")) summary(mc_tukey)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = response ~ trt, data = cholesterol)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1time - drugE == 0 -15.166 1.443 -10.507 < 1e-04 ***
## 2times - drugE == 0 -11.723 1.443 -8.122 < 1e-04 ***
## 4times - drugE == 0 -8.573 1.443 -5.939 < 1e-04 ***
## drugD - drugE == 0 -5.586 1.443 -3.870 0.003063 **
## 2times - 1time == 0 3.443 1.443 2.385 0.137993
## 4times - 1time == 0 6.593 1.443 4.568 0.000344 ***
## drugD - 1time == 0 9.579 1.443 6.637 < 1e-04 ***
## 4times - 2times == 0 3.150 1.443 2.182 0.205053
## drugD - 2times == 0 6.136 1.443 4.251 0.000926 ***
## drugD - 4times == 0 2.986 1.443 2.069 0.251337
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method) On peut également visualiser les comparaisons :
par(mar=c(3,7,3,3))
plot(mc_tukey) 
Le package multcomp, contient également une fonction cld qui permet, dans le cadre du test de Tukey, d’indiquer par des lettres la significativité des comparaisons. Lorsque deux modalités partagent une même lettre, cela signifie que leurs différences ne sont pas significativement différentes. A l’inverse, lorsque deux modalités ne partagent pas de lettres en commun, alors cela signifie que leurs moyennes sont significativement différentes.
tuk.cld <- cld(mc_tukey)
tuk.cld
## drugE 1time 2times 4times drugD
## "d" "a" "ab" "bc" "c" Comme expliqué dans un de mes précédents articles, on peut alors utiliser ces lettres pour les ajouter sur un graph réalisé avec ggplot2.
# permet de mettre les lettres dans un data frame, ensuite utilisé dans ggplot
letters <- tuk.cld$mcletters$Letters
myletters_df <- data.frame(trt=levels(cholesterol$trt),letters=letters)
myletters_df
## trt letters
## drugE drugE d
## 1time 1time a
## 2times 2times ab
## 4times 4times bc
## drugD drugD c
ggplot(cholesterol, aes(x=trt, y=response, colour=trt, fill=trt))+
geom_boxplot(outlier.alpha = 0, alpha=0.25)+
geom_jitter(width=0.25)+
stat_summary(fun.y=mean, colour="black", geom="point", shape=18, size=3) +
theme_classic()+
theme(legend.position="none")+
theme(axis.text.x = element_text(angle=30, hjust=1, vjust=1))+
geom_text(data = myletters_df, aes(label = letters, y = 30 ), colour="black", size=5) 
NB : les losanges noirs représentent les moyennes
J’espère que cet article vous sera utile pour réaliser vos ANOVA à un facteur avec R. Si quelque chose vous freine ou vous pose problème lorsque vous faites ce type d’analyse, indiquez le moi en commentaire, et j’essaierais de vous apporter une réponse.
Dans un prochain article, je vous montrerai quelles sont les solutions envisageables (transformation, ANOVA non-paramétrique, ANOVA basée sur les tests de permutation) lorsque les hypothèses de validité de l’ANOVA paramétrique classique, ne sont pas remplies.
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 
Bonjour, je constate une erreur à la section “5.1 Indépendance des résidus”.
Vous écrivez “Ici la p-value est 0.05 d’où le non-rejet de H0.
Bien cordialement.
Merci pour les explications claires et détaillées.
Au paragraphe 5.1, dans l’interprétation du résultat du test de Durbin-Watson, ne serait ce pas :
” Ici la p-value est > 0.05″ ?
Bonjour,
J’explore votre article “ANOVA à un facteur : partie 2” en essayant de reproduire vos syntaxes de console R.
Le chargement du package multcomp par la commande library (multcomp) me conduit au message d’erreur : Error in library(multcomp) : aucun package nommé ‘multcomp’ n’est trouvé.
Pouvez-vous m’aider ?
Merci d’avance
Bonjour Antoine,
il faut installer le package multcomp avant de le charger dans la session R. Vous pouvez le faire à partir de l’outil de R Studio (fenêtre en bas à droite, onglet package, puis install), ou sinon employer la commande:
install.packages(“multcomp”)
Bonne continuation
Bonjour,
Merci pour votre réponse ; je débute avec R…
Je suis en interface RGui 3.6.0, dans la console R, et je n’ai pas de fenêtre RStudio…
Donc, j’ai cherché multcomp dans la liste des packages :
https://cran.r-project.org/web/packages/
Puis j’ai téléchargé multicomp_1.4-10.zip (pour Windows).
De retour dans R (que j’ai “Executé en tant qu’administrateur”), menu packages / Install packages from local files ; et pointer sur le .zip téléchargé.
Ça a l’air de marcher, mais, pour des questions de dépendances, il faut installer (préalablement) trois autres packages, mvtnorm, MASS et TH.data
Malgré tout, lorsque je fais ?cholesterol la réponse est : No documentation for ‘cholesterol’ in specified packages…
Probablement je me suis totalement fourvoyé…
Merci pour votre aide
Antoine.
Il faut installer Rstudio, cela sera beaucoup plus simple !
Vous trouverez un article pour vous guider ici : https://wp.me/p93iR1-4O
Bonne continuation.
Bonjour,
Merci beaucoup pour tous ces tutos qui me permettent de progresser très rapidement dans l’utilisation de R.
Je trouve que la présentation graphique de l’anova est vraiment très bien et je trouve particulièrement intéressant l’affichage des groupes de traitement sur le graphe. J’aimerai réaliser cette même étude (anova et tests Tuckey) et présenter ce même graphe avec le moyennes des traitements classées par ordre croissant (ou décroissant). Je n’arrive pas à voir à quel niveau faire ce classement, avant l’anova, ou après ?
Merci d’avance pour votre réponse
Olivier
Bonjour Olivier,
il faut le faire avant l’anova, en modifiant les levels de la variable catégoriel, pour les mettre dans l’ordre souhaité.
Comme ceci :
> levels(iris$Species)
[1] “setosa” “versicolor” “virginica”
> library(forcats)
> iris$Species <- fct_relevel(iris$Species, c("versicolor", "virginica", "setosa")) > levels(iris$Species)
[1] “versicolor” “virginica” “setosa”
Bonne continuation.
merci beaucoup pour cette réponse hyper rapide, je vais essayer très vite.
Bonne journée
Bonjour Claire, merci pour l’article.
Je voudrai une aide pour réaliser le SNK.test après l’anova aux mesures répétés
Merci
Bonjour,
J’adore vos articles et je suis vos publications ( y compris twitter ). Je les trouve hyper clairs et détaillés. Bravo !
A la fin de l’article, j’ai du mal à conclure finalement sur le tests d’indépendance de Tukey.
Que concluriez vous en tant que Data Scientiste ?
On voit que DrugE est bien isolé
2time pas fondamentalement différent de 1time et 4time qui eux sont différents 2 à 2 … etc ..
Mais que peut on en tirer comme décision ?
Par ailleurs que signifie la p-value du test de Tukey à chaque ligne ? Ex sur la ligne 1 évalue < 0,05 donc on rejette l'hypothèse nulle comme quoi la 1time et DrugE ont des moyennes semblables ?? Est ce bien cela ?
Merci d'avance
Bonjour Sylvain,
J’en conclus que les deux molécules contrôles (E et D) sont différentes, ce qui est, a priori, problématique. Ensuite que donner 20 mg en 4 prises semble être la meilleure solution car pas de différence mise en évidence avec le contrôle drug D et meilleur que les 20 mg en 2 prises même si la différence n’est pas significative. Mais conclure seulement avec des stats, ce n’est pas une bonne chose à faire, à mon avis. Parce que peut être que 20 mg en 4 prises, dans la réalité c’est compliqué, et que 20 mg en deux prises est peut être un meilleur compromis entre l’efficacité et la compliance au traitement. Je ne sais pas si cela vous aide dans votre réflexion…
Bonjour Claire,
Vos articles sont vraiment très aidants, merci infiniment !
Alice
Bonjour,
Merci pour cet article très clair et utile !
Savez-vous comment changer l’ordre des lettre de significativité pour que sur le graphique, le a soit à gauche, puis vienne le b, etc. ?
Merci !
Bonjour,
non, je ne sais pas, désolée.
Bonne continuation.
Bonjour,
Merci pour votre blog très didactique, une bible !
J’aimerais indiquer par des lettres la significativité des comparaisons avec un témoin (test mc_dunnet), savez-vous comment faire ?
Cordialement
Bonjour,
Vous pouvez regarder les deux articles suivants :
Comment indiquer les résultats des comparaisons sur un graph ?
p-values ajustées : affichez les facilement sur un graph !
Ajouter les pvalues sur un ggplot, manuellement
En espérant que cela vous aide.
Bonne continuation.
Bonjour Claire,
Merci pour votre blog, il m’est d’une grande aide !
J’aimerai générer des lettres après le test de Tukey avec une p-value < 0.01 (et non <0.05)
Es que cela est possible ? Si oui comment faire ?
Merci beaucoup !
Cordialement,
Lucie
Bonjour Lucie,
a priori oui c’est possible avec l’argument level de la fonction cld : tuk.cld <- cld(mc_tukey, level=0.01). Pour trouver cela, j'ai simplement regardé l'aide de la fonction cld. Bonne continuation
Bonjour,
On est plusieurs a faire une remarques sur le 5.1 et la “p-value est < 0.05″ on comprend que ça devrait plutôt être ” p-value est > 0.05″ mais vous ne faites pas de réponse …
Merci pour votre éclairage
Fabien
Bonjour Fabien,
voilà, c’est fait, j’ai modifié la phrase. A présent vous pourrez lire : “Ici la p-value est > 0.05, l’hypothèse H0 n’est donc pas rejetée, et on conclut à l’absence d’auto-corrélation.”
Bonne continuation.
Super ! Merci à vous 🙂
Bonne continuation et merci pour votre partage de toutes ces techniques R
Bonjour,
Quand j’essaie de faire le Test de Tukey sur mes données j’obtiens ceci :
Coefficients:
(Intercept) ParcellePERMA ParcellePYLONE
87.5 145.8 145.8
Alors que j’aimerais avoir :
PERMAPARCELLE 4 PERMAPYLONE PYLONEPARCELLE4
En fait je ne comprend pas pourquoi il prend le titre de la colonne pour faire les coefficients ?
Voici mon jeu de donnée :
Parcelle Effectif
PERMA 325
PERMA 150
PERMA 175
PERMA 175
PERMA 325
PERMA 250
PYLONE 250
PYLONE 175
PYLONE 125
PYLONE 200
PYLONE 450
PYLONE 200
PARCELLE 4 25
PARCELLE 4 0
PARCELLE 4 75
PARCELLE 4 125
PARCELLE 4 150
PARCELLE 4 150
Bonjour,
Tout d’abord, vos effectifs étant du comptage, il faut utiliser un glm et pas une anova simple.
Vous obtenez des coefficients parce que vous avez utilisé la fonction summary. R utilise le titre de la colonne parce que ….c’est son choix ! Les coeffs correspondent à l’écrat de myenne de la modalité avec celle qui sert de référence (première dans l’odre alphabétique, ici parcelle4). Vous pouvez ensuite faire des comparaisons multiples avec le package emmeans. Les pavalues à considérer sont dans la partie $contrasts
parcelles <- read.csv2("data/parcelles.csv") my_anova <- lm(Effectif~Parcelle, data=parcelles) library(car) ## Loading required package: carData Anova(my_anova) ## Anova Table (Type II tests) ## ## Response: Effectif ## Sum Sq Df F value Pr(>F)## Parcelle 72358 2 3.4709 0.0577 .
## Residuals 156354 15
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(my_anova)
##
## Call:
## lm(formula = Effectif ~ Parcelle, data = parcelles)
##
## Residuals:
## Min 1Q Median 3Q Max
## -191.667 -61.458 -6.667 57.708 243.333
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 87.50 41.68 2.099 0.0531 .
## ParcellePERMA 145.83 58.95 2.474 0.0258 *
## ParcellePYLONE 119.17 58.95 2.022 0.0614 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 102.1 on 15 degrees of freedom
## Multiple R-squared: 0.3164, Adjusted R-squared: 0.2252
## F-statistic: 3.471 on 2 and 15 DF, p-value: 0.0577
library(emmeans)
emmeans(my_anova, specs=pairwise~Parcelle)
## $emmeans
## Parcelle emmean SE df lower.CL upper.CL
## PARCELLE 4 87.5 41.7 15 -1.34 176
## PERMA 233.3 41.7 15 144.49 322
## PYLONE 206.7 41.7 15 117.83 296
##
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## PARCELLE 4 - PERMA -145.8 58.9 15 -2.474 0.0630
## PARCELLE 4 - PYLONE -119.2 58.9 15 -2.022 0.1412
## PERMA - PYLONE 26.7 58.9 15 0.452 0.8941
##
## P value adjustment: tukey method for comparing a family of 3 estimates
Avec un glm (https://delladata.fr/tutoriel-glm-sur-donnees-de-comptage-regression-de-poisson-avec-r/)
my_glm <- glm(Effectif~Parcelle, data=parcelles, family=poisson) summary(my_glm) ## ## Call: ## glm(formula = Effectif ~ Parcelle, family = poisson, data = parcelles) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -17.4540 -5.3800 -0.4663 5.1869 14.6171 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.47164 0.04364 102.46 <2e-16 ***
## ParcellePERMA 0.98083 0.05118 19.17 <2e-16 ***
## ParcellePYLONE 0.85947 0.05207 16.51 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1483.1 on 17 degrees of freedom
## Residual deviance: 1023.0 on 15 degrees of freedom
## AIC: 1144.9
##
## Number of Fisher Scoring iterations: 5
my_glm_qp <- glm(Effectif~Parcelle, data=parcelles, family=quasipoisson) Anova(my_glm_qp, test="F") ## Analysis of Deviance Table (Type II tests) ## ## Response: Effectif ## Error estimate based on Pearson residuals ## ## Sum Sq Df F value Pr(>F)
## Parcelle 460.11 2 3.9327 0.04235 *
## Residuals 877.48 15
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(my_glm_qp, specs=pairwise~Parcelle)
## $emmeans
## Parcelle emmean SE df asymp.LCL asymp.UCL
## PARCELLE 4 4.47 0.334 Inf 3.82 5.13
## PERMA 5.45 0.204 Inf 5.05 5.85
## PYLONE 5.33 0.217 Inf 4.91 5.76
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## PARCELLE 4 - PERMA -0.981 0.391 Inf -2.506 0.0327
## PARCELLE 4 - PYLONE -0.859 0.398 Inf -2.158 0.0785
## PERMA - PYLONE 0.121 0.298 Inf 0.407 0.9128
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates
Bonjour Claire,
Merci pour votre blog, tout est vraiment clair étape par étape.
J’aimerais générer une analyse Anova avec une p-value < 0.10 (et non <0.05)
Est ce que cela est possible ?
Merci,
Cordialement,
Alexandre
Bonjour Alexandre,
il suffit de comparer les pvalues obtenues à 0.1 plutôt que 0.05.
Bonne continuation
Merci pour votre réponse.
Du coup ce n’est pas possible par exemple de réaliser un test post-hoc de Tukey et d’avoir ce seuil de 0.1 comme référence afin de pouvoir insérer directement les résultats dans le graphique ?
Bonjour Claire,
Merci pour cet article que je trouve intéressant. J’avais une question mineure à vous poser.
Qu’est ce qui peut expliquer, à la sortie de cld, que certaines modalités ont 2 lettres tandis que d’autres n’en ont qu’une seule? Comment peut-on interpréter la différence de nombre de lettre?
tuk.cld <- cld(mc_tukey)
tuk.cld
## drugE 1time 2times 4times drugD
## "d" "a" "ab" "bc" "c"
Bonjour Lasana
La convention veut que si deux moyennes partagent une même lettre alors elles ne sont pas significativement différentes, et au contraire, si deux moyennes ne partagent pas une même lettre alors elles sont significativement différentes.
cf :
Bonne continuation
Bonjour Claire,
Merci pour tous ces tutoriels très bien fait.
J’ai une question concernant l’indépendance des résidus. Vous dites :”pour réaliser une ANOVA à un facteur, il ne doit pas avoir de données répétées”. Je ne comprend pas ce que vous entendez par la… Pourriez vous me donner un exemple?
Merci d’avance !
Bonjour,
On peut rencontrer des données répétées lorsqu’on observe quelque chose sur les mêmes individus.
Par exemple si vous observez une note (entre 0 et 20) obtenue en math, français et anglais, sur les mêmes individus (vous avez les 3 notes pour chaque individu), et bien, vous ne pourrez pas comparer les moyennes globales des notes en maths, français et anglais avec une anova simple. Il faudra employer une anova sur mesures répétées, car il y a un effet individu (une part de variabilité est causée par les individus).
J’espère que cela vous aide.
Bonne continuation
Merci Claire pour ta réponse,
Si je comprend bien, en prenant l’exemple de suivis marins, si j’échantillonne une fois chaque année la même station pour y compter le taux de recouvrement corallien, je ne pourrais pas comparer les moyennes du taux de recouvrement entre les années avec une ANOVA car c’est la même station?
Merci beaucoup !
Je précise: en effectuant un comptage du benthos sur un transect avec avec un échantillonnage tous les 0,5m.
Avec pour type de résultat: Hard coral Sand Rock Etc..
20 2 15 …
Bonjour Claire ! Je suis buté à petit problème, après avoir fait le test de comparaison multiples comment peut-on identifier la plus petite différence significative (ppds)?
Bonjour,
je pense qu’il faut regarder la colonne estimate qui est la différence entre les 2 groupes comparés;
Bonen continuation.
Bonjour,
Je viens de commencer avec R et souhaite comparer de le recouvrement de mélanges semenciers entre eux. J’ai réussi à suivre votre tutoriel mais j’ai un problème quand j’arrive pour le test de Tukey (ou de Dunett)… Avec l’erreur qui suis: Error in mcp2matrix(model, linfct = linfct) :
Pourtant comme vous pouvez le voir j’ai converti les caractères en facteurs…
Variable(s) ‘Melange’ of class ‘character’ is/are not contained as a factor in ‘model’.
> str(rtot)
‘data.frame’: 16 obs. of 2 variables:
$ Melange : chr “TV2” “TV2” “TV2” “K” …
$ Recouvrement: num 16 25 15 5 5 …
> rtot str(rtot)
‘data.frame’: 16 obs. of 2 variables:
$ Melange : Factor w/ 5 levels “CSR”,”K”,”RR”,..: 5 5 5 2 2 2 2 1 1 1 …
$ Recouvrement: num 16 25 15 5 5 …
> mc_dunnett <- glht(lm1, linfct=mcp(Melange="Dunnett"))
Error in mcp2matrix(model, linfct = linfct) :
Variable(s) ‘Melange’ of class ‘character’ is/are not contained as a factor in ‘model’.
Merci d'avance pour votre aide 🙂
Bonjour,
Est ce que vous avez transformé Melange en factor avant de faire l’ANOVA ?
Je suspecte que non et que c’est un problème pour la fonction glht()
Bonne continuation
Effectivement le problème venais de l’exécution de la transformation après avoir fait lm1.
Merci pour votre aide et votre réactivité.
Bonjour Claire,
Merc infiniment,
cependant, je sais face a une difficulté décrite ci-dessous
> mc_dunnett <- glht(lm2, linfct=mcp(Engrais="Dunnett"))
Error in glht(lm2, linfct = mcp(Engrais = "Dunnett")) :
could not find function "glht"
j'a cru que le problème &tait lié à la fonction glht.
donc jai essayé de l'installer mais le problème persiste.
merci d'avance, cordiallement.
Bonjour,
Vérifiez que vous avez bien installé le package multcomp.
Sinon, regardez du côté du package emmeans.
Bonne continuation.
Bonjour,
A priori le lien vers l’article de la partie 1 ne semble pas fonctionner (https://delladata.fr/anova-a-un-facteur-partie-1)
A bravo pour tous les articles, je suis un fidèle de blog depuis des années !