ANOVA à un facteur avec R
Dhafer Malouche
3ième année ESSAI, 2015-2016
Objectif
L’analyse de la variance (ANOVA) est une méthode qui permet d’étudier la variantion de la moyenne \(\mu\) du phénoméne étudié \(Y\) (variable quantitative) selon l’influence d’un ou de plusieurs facteurs d’expérience qualitatifs (traitements…).
Dans le cas où la moyenne n’est influencée que par un seul facteur (noté \(X\)), il s’agit d’une analyse de la variance à un facteur (one way ANOVA). Un facteur est une variable qualitative présentant un nombre restreint de modalités.
Le nombre de modalités (c’est-a-dire de niveaux) du facteur \(X\) sera noté \(I\). On suppose que \(Y\) suit une loi normale \({\mathcal N}(\mu_i, \sigma^2)\) sur chaque sous-population \(i\) définie par les modalités de \(X.\)
L’objectif est de tester les l’égalité des moyennes de ces \(I\) populations, à savoir tester l’hypothése nulle : \[H_0\,:\, \mu_1=\ldots=\mu_I,\;\;\mbox{ vs}\;\;H_1\;:\;\exists \,i\not j \mu_i\not=\mu_j.\]
Les données
Pour chaque type de population \(i\) (ou modalité \(i\) de \(X\) ou groupe \(i\)), on dispose d’un échantillon de \(n_i\) observations de la variable quantitative \(Y\) : \[y_{i,1},\ldots,y_{i,n_i}\]
Exemple
On considère une base de données sur des scores de préference de deux variétés de huiles d’olive. Ces scores sont donnés par des panels d’expert.
> sensoHuile <- read.csv("~/Documents/ESSAI_enseig_pfe/MesCourStatAvecR/FormationR_EcolePrintemps/sensoHuile.csv")
> head(sensoHuile)
Echant Product Variety locality System Jury Fusty Moldy Viney
1 2 cheml,bouz,sp cheml bouzC sp A 0 0 0
2 2 cheml,bouz,sp cheml bouzC sp B 0 0 0
3 2 cheml,bouz,sp cheml bouzC sp C 0 0 0
4 2 cheml,bouz,sp cheml bouzC sp D 0 0 0
5 2 cheml,bouz,sp cheml bouzC sp E 0 0 0
6 2 cheml,bouz,sp cheml bouzC sp F 0 0 0
Muddy Métali Rancid Fruitty Bitter Pungent
1 0 0 0 2.0 1.23 0.50
2 0 0 0 2.5 1.55 0.70
3 0 0 0 3.0 1.50 0.55
4 0 0 0 2.8 1.30 0.80
5 0 0 0 3.0 1.46 0.50
6 0 0 0 4.5 2.12 1.50
> attach(sensoHuile)On veut déterminer s’il y a une différence significative dans le goût fruité entre les différentes variétés.
- Une représentation graphique s’impose
- Une matrice résumant les statistiques par variété
> library(doBy)
Loading required package: survival
> cdata <- summaryBy(Fruitty ~ Variety, data=sensoHuile, FUN=c(length,mean,sd))
> colnames(cdata)=c("Variety","N","Mean","sd")
> cdata=format(cdata,digits=3)
> cdata
Variety N Mean sd
1 arbeq 144 2.89 1.34
2 cheml 136 3.22 1.07- Le graphe des boxplots
- Test d’ANOVA
> mon.aov <- aov(Fruitty~Variety)
> summary(mon.aov)
Df Sum Sq Mean Sq F value Pr(>F)
Variety 1 7.5 7.483 5.052 0.0254 *
Residuals 278 411.8 1.481
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- En effet ce tableau est une application du theéoreeème plus geéneéral.
Théorème On consideère un modeèle de reégression \(\mathcal{M}\,:\, y=X\beta +\epsilon\) qui contient \(p\) variables. Nous souhaitons tester la validité d’un sous-modèle (ou modèle emboîté) où un ou plusieurs coefficients sont nuls. Le plan d’expérience privé de ces variables sera noté \(X_0\) ouù \(X_0\) est une matrice \(n\times p_0\) et et le sous-modèle sera \[\mathcal{M}_0\,:\, y=X_0\beta_0 +\epsilon\] Donc nous voulons tester les deux hypotheèses \[\left\{\begin{array}{ll} \mbox{H}_0\;:\; E(y)\in \mathcal{I}(X_0)&\mbox{modèle réduit}\\ \mbox{H}_1\;:\; E(y)\in \mathcal{I}(X) & \mbox{modèle complet} \end{array}\right.\] Pour tester ces deux hypothèses, nous utilisons la statistique de test \(F\) ci-dessous qui possède comme loi sous \(\mbox{H}_0\) : \[F=\displaystyle\frac{||\widehat{Y}-\widehat{Y}_0||^2/(p-p_0)}{||Y-\widehat{Y}||^2/(n-p)}\sim\mbox{F}(p-p_0,n-p)\] Remarquons que $||-_0||=-_0|| où \[\mbox{SCR}=||Y-\widehat{Y}||^2\]
Dans le cas d’une ANOVA d’un facteur à \(I\) modalités, Le modèle réduit est \[y_i=\mu+\epsilon_i,\;\;i=1,\ldots,n\] Donc \(p_0=1\), et on montre \[\widehat{\mu}=\overline{y}\;\;\mbox{et}\;\; \widehat{\sigma}^2_0=\displaystyle\frac{1}{n-1}\sum_{i=1}^I\sum_{j=1}^{n_i}(y_{i,j}-\overline{y})^2,\] \[||\widehat{Y}-\widehat{Y}_0||^2=\sum_{i=1}^In_i(\overline{y}_i-\overline{y})^2\] et \[||Y-\widehat{Y}||^2=\sum_{i=1}^I\sum_{j=1}^{n_i}(y_{i,j}-\overline{y})^2\]
Retrouvons le calcul du tableau dans 3.
> m1=lm(Fruitty~Variety,data=sensoHuile)
> m0=lm(Fruitty~1,data=sensoHuile)
> y0hat=predict(m0,sensoHuile)
> y0hat[1:10]
1 2 3 4 5 6 7 8
3.053107 3.053107 3.053107 3.053107 3.053107 3.053107 3.053107 3.053107
9 10
3.053107 3.053107
> mean(sensoHuile$Fruitty)
[1] 3.053107
> yhat=predict(m1,sensoHuile)
> sum((yhat-y0hat)^2)
[1] 7.482919
> sum((sensoHuile$Fruitty-yhat)^2)
[1] 411.7783
> p0=1
> n=nrow(sensoHuile)
> p=2
> sum((sensoHuile$Fruitty-yhat)^2)/(n-p)
[1] 1.481217
> sum((yhat-y0hat)^2)/p0
[1] 7.482919
> F=(sum((yhat-y0hat)^2)/p0)/(sum((sensoHuile$Fruitty-yhat)^2)/(n-p));F
[1] 5.051873
> 1-pf(F,p0,n-p,lower.tail = T)
[1] 0.02538342Validation des hypothéses : Egalité de la variance
> bartlett.test(Fruitty~Variety)
Bartlett test of homogeneity of variances
data: Fruitty by Variety
Bartlett's K-squared = 6.8875, df = 1, p-value = 0.00868Deux solutions sont possibles
- Transformation de la variable `
Fruitty
> LFruitty=log(Fruitty)
> bartlett.test(LFruitty~Variety)
Bartlett test of homogeneity of variances
data: LFruitty by Variety
Bartlett's K-squared = 11.383, df = 1, p-value = 0.0007412Ou on peut utiliser un test (le test de Levene) plus robuste (moins sensible à la condition de normalité)
> library(car)
> leveneTest(LFruitty,Variety)
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 13.226 0.0003289 ***
278
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Utiliser un autre test qui tient compte de l’inégalité des variances (test de Welch).
> oneway.test(Fruitty~Variety, var.equal = F)
One-way analysis of means
data: Fruitty and Variety
F = 5.0519, num df = 1, denom df = 278, p-value = 0.02538Comparaison entre deux facteurs.
On va s’intéresser maintenant aux deux facteurs System et Variety
> library(doBy)
> int.conf=function(x){
+ z=t.test(x)
+ return(z$conf.int)
+ }
> cdata <- summaryBy(Fruitty ~ Variety+System, data=sensoHuile, FUN=c(length,mean,sd,int.conf))
> colnames(cdata)=c("Variety","System","N","Mean","sd","Lower","Upper")
> cdata1=format(cdata,digits=3)
> cdata1
Variety System N Mean sd Lower Upper
1 arbeq 2p 48 2.92 1.269 2.55 3.28
2 arbeq 3p 48 3.79 1.017 3.50 4.09
3 arbeq sp 48 1.98 1.065 1.67 2.28
4 cheml 2p 48 3.72 0.977 3.44 4.00
5 cheml 3p 48 3.71 0.527 3.56 3.86
6 cheml sp 40 2.04 0.669 1.82 2.25Un graphique des intervalles de confiances
> p1<-ggplot(cdata, aes(x=System, y=Mean)) +
+ geom_errorbar(aes(ymin=Lower, ymax=Upper), width=.2) +
+ geom_line(size=8) +
+ geom_point(size=8)+facet_wrap(~Variety)
> p1
ou
> pd <- position_dodge(0.4)
> p1<-ggplot(cdata, aes(x=System, y=Mean,group=Variety,color=Variety)) +
+ geom_errorbar(aes(y=Mean,ymin=Lower, ymax=Upper), width=.2,position=pd) +
+ geom_point(size=4,position=pd)
> p1
ymax not defined: adjusting position using y instead
L’ANOVA à deux facteurs
> mon.aov <- aov(Fruitty~Variety+System)
> summary(mon.aov)
Df Sum Sq Mean Sq F value Pr(>F)
Variety 1 7.48 7.48 7.825 0.00551 **
System 2 147.86 73.93 77.315 < 2e-16 ***
Residuals 276 263.92 0.96
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Comparaison deux à deux
On veut effectuer toutes les comparaisons deux à deux en proposant plusieurs méthodes de correction du risque afin de tenir compte du problème de la multiplicité des tests.
> TukeyHSD(mon.aov)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Fruitty ~ Variety + System)
$Variety
diff lwr upr p adj
cheml-arbeq 0.3270874 0.09690877 0.5572661 0.005514
$System
diff lwr upr p adj
3p-2p 0.4315625 0.09896904 0.7641560 0.0068994
sp-2p -1.3001608 -1.64022920 -0.9600924 0.0000000
sp-3p -1.7317233 -2.07179170 -1.3916549 0.0000000Comparaison deux à deux (par Variety)
- Deux ANOVA
> mon.aov1 <- aov(Fruitty~System,data=sensoHuile[Variety=='cheml',])
> summary(mon.aov1)
Df Sum Sq Mean Sq F value Pr(>F)
System 2 79.42 39.71 70.03 <2e-16 ***
Residuals 133 75.41 0.57
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> mon.aov2 <- aov(Fruitty~System,data=sensoHuile[Variety=='arbeq',])
> summary(mon.aov2)
Df Sum Sq Mean Sq F value Pr(>F)
System 2 79.24 39.62 31.44 5.13e-12 ***
Residuals 141 177.71 1.26
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Comparaison deux à deux
> TukeyHSD(mon.aov1)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Fruitty ~ System, data = sensoHuile[Variety == "cheml", ])
$System
diff lwr upr p adj
3p-2p -0.012500 -0.376813 0.351813 0.9963604
sp-2p -1.683333 -2.065428 -1.301239 0.0000000
sp-3p -1.670833 -2.052928 -1.288739 0.0000000
> TukeyHSD(mon.aov2)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Fruitty ~ System, data = sensoHuile[Variety == "arbeq", ])
$System
diff lwr upr p adj
3p-2p 0.8756250 0.3328071 1.4184429 0.0005793
sp-2p -0.9410417 -1.4838596 -0.3982238 0.0001992
sp-3p -1.8166667 -2.3594846 -1.2738488 0.0000000- Un graphique pour finir.
> x1=TukeyHSD(mon.aov1)$System
> x2=TukeyHSD(mon.aov2)$System
> cdata=rbind.data.frame(x1,x2)
> cdata=cbind.data.frame(c(rep('Chemlali',3),rep('Arbequina',3)),rownames(cdata),cdata)
> colnames(cdata)[1:2]=c('Variety','System')
> cdata$System[4:6]=cdata$System[1:3]
> cdata
Variety System diff lwr upr p adj
3p-2p Chemlali 3p-2p -0.0125000 -0.3768130 0.3518130 9.963604e-01
sp-2p Chemlali sp-2p -1.6833333 -2.0654280 -1.3012386 0.000000e+00
sp-3p Chemlali sp-3p -1.6708333 -2.0529280 -1.2887386 0.000000e+00
3p-2p1 Arbequina 3p-2p 0.8756250 0.3328071 1.4184429 5.792867e-04
sp-2p1 Arbequina sp-2p -0.9410417 -1.4838596 -0.3982238 1.992048e-04
sp-3p1 Arbequina sp-3p -1.8166667 -2.3594846 -1.2738488 1.828981e-12
> p1<-ggplot(cdata, aes(x=System, y=diff)) +
+ geom_errorbar(aes(ymin=lwr, ymax=upr), width=.2) +
+ geom_line(size=8) +
+ geom_point(size=8)+facet_wrap(~Variety)+geom_hline(y=0,type=2,colour='gray')
> p1
geom_path: Each group consist of only one observation. Do you need to adjust the group aesthetic?
geom_path: Each group consist of only one observation. Do you need to adjust the group aesthetic?
Utilisant la commande pairwise.t.test()
La fonction pairwise.t.test() permet d’effectuer toutes les comparaisons deux é deux en proposant plusieurs méthodes de correction du risque \(\alpha\) afin de tenir compte du probléme de la multiplicité des tests
> pair1=pairwise.t.test(Fruitty[Variety=='cheml'],System[Variety=='cheml'],p.adjust="bonf",pool.sd = F)
> pair1
Pairwise comparisons using t tests with pooled SD
data: Fruitty[Variety == "cheml"] and System[Variety == "cheml"]
2p 3p
3p 1 -
sp <2e-16 <2e-16
P value adjustment method: bonferroni
> pair2=pairwise.t.test(Fruitty[Variety=='arbeq'],System[Variety=='arbeq'],p.adjust="bonf",pool.sd = F)
> pair2
Pairwise comparisons using t tests with pooled SD
data: Fruitty[Variety == "arbeq"] and System[Variety == "arbeq"]
2p 3p
3p 6e-04 -
sp 2e-04 1.8e-12
P value adjustment method: bonferroni Le package agricolae
> library(agricolae)
> HSD.test(mon.aov1,"System", group=TRUE,console = T)
Study: mon.aov1 ~ "System"
HSD Test for Fruitty
Mean Square Error: 0.5669893
System, means
Fruitty std r Min Max
2p 3.720833 0.9773952 48 1.7 5.1
3p 3.708333 0.5274722 48 2.5 4.8
sp 2.037500 0.6685950 40 1.1 4.5
alpha: 0.05 ; Df Error: 133
Critical Value of Studentized Range: 3.352029
Harmonic Mean of Cell Sizes 45
Honestly Significant Difference: 0.3762608
Means with the same letter are not significantly different.
Groups, Treatments and means
a 2p 3.721
a 3p 3.708
b sp 2.038
> HSD.test(mon.aov2,"System", group=TRUE,console = T)
Study: mon.aov2 ~ "System"
HSD Test for Fruitty
Mean Square Error: 1.260349
System, means
Fruitty std r Min Max
2p 2.916042 1.269298 48 1 5.0
3p 3.791667 1.017262 48 2 5.3
sp 1.975000 1.065414 48 1 4.5
alpha: 0.05 ; Df Error: 141
Critical Value of Studentized Range: 3.349881
Honestly Significant Difference: 0.5428179
Means with the same letter are not significantly different.
Groups, Treatments and means
a 3p 3.792
b 2p 2.916
c sp 1.975 ou encore
> library(multcomp)
Loading required package: mvtnorm
Loading required package: TH.data
> tuk1 <- glht(mon.aov1, linfct = mcp(System = "Tukey"))
> summary(tuk1)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: aov(formula = Fruitty ~ System, data = sensoHuile[Variety ==
"cheml", ])
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
3p - 2p == 0 -0.0125 0.1537 -0.081 0.996
sp - 2p == 0 -1.6833 0.1612 -10.442 <1e-06 ***
sp - 3p == 0 -1.6708 0.1612 -10.365 <1e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)
> tuk1.cld <- cld(tuk1)
> tuk2 <- glht(mon.aov2, linfct = mcp(System = "Tukey"))
> summary(tuk2)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: aov(formula = Fruitty ~ System, data = sensoHuile[Variety ==
"arbeq", ])
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
3p - 2p == 0 0.8756 0.2292 3.821 0.000598 ***
sp - 2p == 0 -0.9410 0.2292 -4.106 0.000199 ***
sp - 3p == 0 -1.8167 0.2292 -7.927 < 1e-04 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)
> tuk1.cld <- cld(tuk1)
> tuk2.cld <- cld(tuk2)
> opar <- par(mai=c(1,1,1.5,1),mfrow = c(1,2))
> cols1=as.numeric(factor(tuk1.cld$mcletters$Letters))+1
> cols2=as.numeric(factor(tuk2.cld$mcletters$Letters))+1
> plot(tuk1.cld,col=cols1,xlab="Chemlali")
> plot(tuk2.cld,col=cols2,xlab="Arbequina")
> par(opar)