ANOVA

• Analyse de variance ou ANOVA
• One-way analysis of variance (ANOVA)
• One-factor ANOVA

sur la base des cours de L. Reboul et du “R book”.

ANOVA

Une extension de la comparaison de moyennes quand il y a plus de deux groupes. Les groupes sont des facteurs (factor), nous allons tester si la moyenne des différents groupes est la même ou pas. L’hypothèse Ho est qu’il y a au moins un groupe dont la moyenne n’est pas égale à celle des autres groupes.

ANOVA : hypothèses

• les observations sont indépendantes
• pour chaque groupe (factor), les données suivent une loi Normale
• pour chaque groupe, la variance est la même

ANOVA et R

df <- data.frame(
  myGp = rep(1:3, times = 30),
  myData = rnorm(90)
)
print(df)

##    myGp       myData
## 1     1  0.231054022
## 2     2 -0.289179134
## 3     3 -1.280162893
## 4     1  1.072713244
## 5     2  0.215148253
## 6     3 -1.041524613
## 7     1 -1.592904833
## 8     2  0.237602054
## 9     3  0.507547292
## 10    1 -0.436476788
## 11    2 -0.827732999
## 12    3  1.816606232
## 13    1  0.496288170
## 14    2 -0.473338790
## 15    3  0.096940305
## 16    1 -0.565519591
## 17    2  1.183138913
## 18    3 -0.157524688
## 19    1 -0.639482046
## 20    2 -0.999999966
## 21    3 -0.457214343
## 22    1 -0.015784216
## 23    2  0.554783877
## 24    3  1.044781271
## 25    1  0.551952859
## 26    2  0.336166229
## 27    3 -2.143860102
## 28    1  0.068475877
## 29    2 -1.019164615
## 30    3  0.351129907
## 31    1 -1.093372937
## 32    2  0.367960655
## 33    3  0.613626848
## 34    1  0.803314853
## 35    2 -0.351708575
## 36    3  0.644295080
## 37    1  1.629525175
## 38    2 -1.276627060
## 39    3  1.010594324
## 40    1 -0.478130690
## 41    2 -1.137574471
## 42    3 -1.850821962
## 43    1 -1.636169441
## 44    2  0.015931003
## 45    3  0.004911134
## 46    1 -0.520214353
## 47    2 -1.810873680
## 48    3 -1.727836718
## 49    1  1.102038725
## 50    2 -1.421477370
## 51    3  0.725185746
## 52    1  0.133766537
## 53    2  0.695776695
## 54    3 -0.063354107
## 55    1  0.290947665
## 56    2 -0.102632573
## 57    3  1.868037279
## 58    1 -2.225499922
## 59    2 -0.700208123
## 60    3  1.020031993
## 61    1 -0.539738946
## 62    2  2.735205972
## 63    3  0.957983197
## 64    1 -0.089488483
## 65    2  0.523134904
## 66    3  0.519632447
## 67    1 -1.081405571
## 68    2 -0.598408912
## 69    3  0.764627767
## 70    1 -1.337090605
## 71    2  1.170173433
## 72    3 -0.669016413
## 73    1  0.395627675
## 74    2  0.930320959
## 75    3 -0.634532996
## 76    1  0.515346091
## 77    2 -0.206944827
## 78    3  0.892421647
## 79    1 -0.342180830
## 80    2 -0.216201463
## 81    3  1.691006745
## 82    1 -0.412201398
## 83    2  0.092082861
## 84    3  0.204286020
## 85    1 -0.314259735
## 86    2  0.680878209
## 87    3 -1.607875073
## 88    1 -0.376465894
## 89    2 -0.072720816
## 90    3 -0.109153163

ANOVA et R

str(df)

## 'data.frame':    90 obs. of  2 variables:
##  $ myGp  : int  1 2 3 1 2 3 1 2 3 1 ...
##  $ myData: num  0.231 -0.289 -1.28 1.073 0.215 ...

Le vecteur qui détermine le groupe doit être de type factor !

• Eviter les groupes 1, 2, 3
• TOUJOURS vérifier le type de données avec str()

ANOVA et R

df$myGp <- as.factor(df$myGp)
str(df)

## 'data.frame':    90 obs. of  2 variables:
##  $ myGp  : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...
##  $ myData: num  0.231 -0.289 -1.28 1.073 0.215 ...

ANOVA

boxplot(df$myData ~ df$myGp, 
  col = levels(df$myGp))

ANOVA

boxplot(df$myData ~ df$myGp, 
  col = as.numeric(levels(df$myGp)) + 1)

ANOVA

myCol <- c("olivedrab1", "salmon", "slateblue4")
boxplot(df$myData ~ df$myGp, 
  col = myCol[as.numeric(levels(df$myGp))])

ANOVA

myCol <- c("olivedrab1", "salmon", "slateblue4")
boxplot(df$myData ~ df$myGp, 
  col = myCol[as.numeric(levels(df$myGp))], 
  names = paste0("trt_",levels(df$myGp)))

ANOVA

• Variable explicative qualitative (type factor) : “1”, “3”, ou “2”
• Variable à expliquer quantitative continue : tirage au hasard dans une loi Normale de paramètres μ=0 et σ=1.

Ici nous avons 1 facteur qui prend comme niveaux (ou modalités) “1”, “2”, “3” (p=3), et pour chaque niveau un échantillon de taille 30 (plan équilibré n1=...=np).

ANOVA

tapply(df$myData, INDEX = df$myGp, FUN = length)

##  1  2  3 
## 30 30 30

tapply(df$myData, INDEX = df$myGp, FUN = summary)

## $`1`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -2.2255 -0.5591 -0.3282 -0.2135  0.3695  1.6295 
## 
## $`2`
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -1.81087 -0.67476 -0.08768 -0.05888  0.48434  2.73521 
## 
## $`3`
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -2.14386 -0.59020  0.27771  0.09969  0.86047  1.86804

ANOVA

Nous allons tester si la moyenne des différents groupes est la même ou pas = mettre en évidence une différence de valeur moyenne selon le niveau du facteur.

Yj=μ+αj+ϵj,j=1,...,p

H0:α1=...=αj=...=αP

ANOVA

On suppose :

yij=μ+αj+ϵij,i=1,...,nj,j=1,...,p

ANOVA

Dans notre modèle, on ne connait pas μ, et α1,...,αp, ce qui fait p+1 inconnus. Comme nous n’avons que p groupes, on impose que :

∑pj=1njαj=0

compensation entre groupes

ANOVA

On minimise (moindres carrés) :

∑i∑j(yij−αj−μ)2

ANOVA

Somme des carrés des variations dans le niveau j

Dj=∑(yij−¯yj)2

Dj <- sapply(levels(df$myGp), function(j){
  sum((df$myData[df$myGp == j] - 
    mean(df$myData[df$myGp == j]))^2)
})
print(Dj)

##        1        2        3 
## 21.99745 25.23369 34.23868

ANOVA

Somme des carrés des variations intra-niveaux

SSintra=∑jDj=(n−p)S2

SSintra <- sum(Dj)
print(SSintra)

## [1] 81.46982

ANOVA

Somme des carrés des variations inter-niveaux (SSE)

SSinter=∑jnj(¯yj−¯y)2

SSinter <- sum(sapply(levels(df$myGp), function(j){
  length(df$myData[df$myGp == j]) * 
    (mean(df$myData[df$myGp == j]) - 
      mean(df$myData))^2
}))
print(SSinter)

## [1] 1.471524

ANOVA

Somme des carrés des variations totales (SST)

SST=∑(yij−¯y)2

SST <- sum((df$myData - mean(df$myData))^2)
print(SST)

## [1] 82.94134

ANOVA

SST=SSinter+SSintra

print(SSinter + SSintra)

## [1] 82.94134

print(SST)

## [1] 82.94134

ANOVA

Indice de qualité : pourcentage de variation expliqué par le modèle

η2=SSinter/SST

eta2 <- SSinter / SST
print(eta2)

## [1] 0.01774174

ANOVA

Indice de Fisher

F=SSinter/(p−1)SSintra/n−p

Fi <- (SSinter/(3-1)) / (SSintra/(90-3))
print(Fi)

## [1] 0.7857056

ANOVA en images… (the R book)

plot(df$myData, type = 'n', axes = FALSE, 
  xlab = "", ylab = "")
axis(1)
axis(2)
df <- df[order(df$myGp),]
rownames(df) <- NULL
trash <- lapply(levels(df$myGp), function(p){
  points(
    x = as.numeric(rownames(df[df$myGp == p,])), 
    y = df$myData[df$myGp == p], 
    col = as.numeric(p), pch = 16)
})

ANOVA en images…

ANOVA en images…

plot(df$myData, type = 'n', axes = FALSE, 
  xlab = "", ylab = "")
axis(1)
axis(2)
df <- df[order(df$myGp),]
rownames(df) <- NULL
trash <- lapply(levels(df$myGp), function(p){
  points(
    x = as.numeric(rownames(df[df$myGp == p,])), 
    y = df$myData[df$myGp == p], 
    col = as.numeric(p), pch = 16)
})
abline(h = mean(df$myData))
segments(
  x0 = 1:nrow(df), 
  y0 = rep(mean(df$myData), nrow(df)), 
  x1 = 1:nrow(df), 
  y1 = df$myData)

ANOVA en images…

ANOVA en images…

Représentation de la variance totale.

Somme des carrés des variations totales (SST)

SST=∑(yij−¯y)2

SST <- sum((df$myData - mean(df$myData))^2)
print(SST)

## [1] 82.94134

ANOVA en images…

plot(df$myData, type = 'n', axes = FALSE, 
  xlab = "", ylab = "")
axis(1)
axis(2)
df <- df[order(df$myGp),]
rownames(df) <- NULL
trash <- lapply(levels(df$myGp), function(p){
  points(
    x = as.numeric(rownames(df[df$myGp == p,])), 
    y = df$myData[df$myGp == p], 
    col = as.numeric(p), pch = 16)
})
trash <- lapply(levels(df$myGp), function(p){
  points(
    x = as.numeric(rownames(df[df$myGp == p,])), 
    y = df$myData[df$myGp == p], 
    col = as.numeric(p), pch = 16)
  segments(
    x0 = as.numeric(rownames(df[df$myGp == p,])),
    y0 = mean(df$myData[df$myGp == p]),
    x1 = as.numeric(rownames(df[df$myGp == p,])),
    y1 = df$myData[df$myGp == p], col = p
  )
  segments(
    x0 = as.numeric(rownames(df[df$myGp == p,]))[1],
    x1 = tail(as.numeric(rownames(df[df$myGp == p,])),1),
    y0 = mean(df$myData[df$myGp == p]),
    y1 = mean(df$myData[df$myGp == p]), col = p
  )
})

ANOVA en images…

ANOVA en images…

Somme des carrés des variations intra-niveaux

SSintra=∑jDj=(n−p)S2

Dj=∑(yij−¯yj)2

Dj <- sapply(levels(df$myGp), function(j){
  sum((df$myData[df$myGp == j] - 
    mean(df$myData[df$myGp == j]))^2)
})
SSintra <- sum(Dj)
print(SSintra)

## [1] 81.46982

ANOVA en images…

Si SST=SSintra, pas de différence entre moyennes.

Calcul de SSinter=SST−SSintra

=> quel seuil pour différence significqtive, test de Fisher sur les variances.

ANOVA

Indice de Fisher

F=SSinter/(p−1)SSintra/n−p

Fi <- (SSinter/(3-1)) / (SSintra/(90-3))
print(Fi)

## [1] 0.7857056

# seuil critique 5% (F(p-1;n-p))
qf(p=0.95, df1=3-1, df2=90-3)

## [1] 3.101296

ANOVA

Validation des hypothèses

• échantillonnage aléatoire
• variances égales
• indépendance des erreurs
• distribution normale des erreurs

ANOVA ; variances

tapply(df$myData, INDEX = df$myGp, FUN = var)

##         1         2         3 
## 0.7585327 0.8701273 1.1806440

# Flinger-Killen : homogeneité des variances
fligner.test(df$myData ~ df$myGp)

## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  df$myData by df$myGp
## Fligner-Killeen:med chi-squared = 1.836, df = 2, p-value = 0.3993

ANOVA

summary(aov(df$myData ~ df$myGp))

##             Df Sum Sq Mean Sq F value Pr(>F)
## df$myGp      2   1.47  0.7358   0.786  0.459
## Residuals   87  81.47  0.9364

ANOVA ; hypothèses

par(mfrow = c(2, 2), ask = FALSE)
plot(aov(df$myData ~ df$myGp))

ANOVA ; hypothèses

graph 1 : variances égales

graph 2 : distribution normale des erreurs

grpah 3/4 : indépendance des erreurs

ANOVA 2f

• ANOVA à 2 facteurs
• Two-way ANOVA

ANOVA 2f

df2 <- data.frame(
  myData = rnorm(360),
  myGp1 = rep(c("t1", "t2", "t3"), each = 120),
  myGp2 = rep(c("x1", "x2", "x3", "x4"), times = 90)
)

ANOVA 2f : hypothèses

• Il n’y a pas de différence entre moyennes pour le groupe 1
• Il n’y a pas de différence entre moyennes pour le groupe 2
• Il n’y a pas d’interaction antre les groupes 1 et 2

ANOVA 2f

Cas d’un design équilibré 3x4 : 30 observations par combinaison de facteurs.

table(df2$myGp1, df2$myGp2)

##     
##      x1 x2 x3 x4
##   t1 30 30 30 30
##   t2 30 30 30 30
##   t3 30 30 30 30

ANOVA 2f

boxplot(df2$myData ~ df2$myGp1, 
  col = c(2:4), las = 3)

ANOVA 2f

boxplot(df2$myData ~ df2$myGp1 + df2$myGp2, 
  col = c(2:4), las = 3)
abline(v = c(3.5, 6.5, 9.5))

ANOVA 2f

boxplot(df2$myData ~ df2$myGp2, 
  col = c(2:5), las = 3)

ANOVA 2f

boxplot(df2$myData ~ df2$myGp2 + df2$myGp1, 
  col = c(2:5), las = 3)
abline(v = c(4.5, 8.5))

ANOVA 2f

interaction.plot(
  x.factor = df2$myGp1, 
  trace.factor = df2$myGp2, 
  response = df2$myData, 
  fun = mean, 
  type = "b", legend = TRUE, 
  xlab = "Treatment", ylab="Results",
  pch = c(1, 19), col = 1:4)

ANOVA 2f

ANOVA 2f (additif : pas d’interaction)

res2 <- aov(myData ~ myGp1 + myGp2, data = df2)
summary(res2)

##              Df Sum Sq Mean Sq F value Pr(>F)
## myGp1         2    0.8  0.3806   0.338  0.713
## myGp2         3    3.4  1.1412   1.014  0.386
## Residuals   354  398.4  1.1254

ANOVA 2f (test interaction)

res3 <- aov(myData ~ myGp1 * myGp2, data = df2)
summary(res3)

##              Df Sum Sq Mean Sq F value Pr(>F)
## myGp1         2    0.8  0.3806   0.338  0.714
## myGp2         3    3.4  1.1412   1.012  0.387
## myGp1:myGp2   6    5.9  0.9843   0.873  0.515
## Residuals   348  392.5  1.1278

ANOVA 2f (stats descriptives)

model.tables(res2, type = "means", se = TRUE)

## Tables of means
## Grand mean
##             
## -0.01544385 
## 
##  myGp1 
## myGp1
##       t1       t2       t3 
## -0.02459 -0.06664  0.04489 
## 
##  myGp2 
## myGp2
##       x1       x2       x3       x4 
## -0.18015  0.03959  0.06999  0.00879 
## 
## Standard errors for differences of means
##          myGp1  myGp2
##         0.1370 0.1581
## replic.    120     90

ANOVA 2f + Tukey Honest Significant Differences

Comparaisons mutliples 2 à 2.

TukeyHSD(res2)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = myData ~ myGp1 + myGp2, data = df2)
## 
## $myGp1
##              diff        lwr       upr     p adj
## t2-t1 -0.04204985 -0.3643853 0.2802856 0.9493698
## t3-t1  0.06947382 -0.2528616 0.3918092 0.8678456
## t3-t2  0.11152367 -0.2108118 0.4338591 0.6943992
## 
## $myGp2
##              diff        lwr       upr     p adj
## x2-x1  0.21973424 -0.1884690 0.6279375 0.5067607
## x3-x1  0.25013522 -0.1580680 0.6583385 0.3903686
## x4-x1  0.18893809 -0.2192652 0.5971413 0.6305677
## x3-x2  0.03040098 -0.3778023 0.4386042 0.9974783
## x4-x2 -0.03079615 -0.4389994 0.3774071 0.9973796
## x4-x3 -0.06119713 -0.4694004 0.3470061 0.9802826

ANOVA 2f

df3 <- data.frame(
  myData = c(
    rnorm(120), 
    rnorm(120, mean = 2), 
    rnorm(120, mean = 5)),
  myGp1 = rep(c("t1", "t2", "t3"), each = 120),
  myGp2 = rep(c("x1", "x2", "x3", "x4"), times = 90)
)

ANOVA 2f

res4 <- aov(myData ~ myGp1 + myGp2, data = df3)
summary(res4)

##              Df Sum Sq Mean Sq F value Pr(>F)    
## myGp1         2 1621.0   810.5 783.344 <2e-16 ***
## myGp2         3    4.2     1.4   1.339  0.262    
## Residuals   354  366.3     1.0                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANOVA 2f

boxplot(df3$myData ~ df3$myGp2 + df3$myGp1, 
  col = c(2:5), las = 3)
abline(v = c(4.5, 8.5))

ANOVA 2f

TukeyHSD(res4)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = myData ~ myGp1 + myGp2, data = df3)
## 
## $myGp1
##           diff      lwr      upr p adj
## t2-t1 1.951714 1.642644 2.260783     0
## t3-t1 5.147817 4.838748 5.456886     0
## t3-t2 3.196103 2.887034 3.505173     0
## 
## $myGp2
##              diff        lwr       upr     p adj
## x2-x1  0.15256917 -0.2388340 0.5439724 0.7459037
## x3-x1 -0.10751430 -0.4989175 0.2838889 0.8934836
## x4-x1  0.14156625 -0.2498369 0.5329694 0.7868022
## x3-x2 -0.26008347 -0.6514867 0.1313197 0.3172379
## x4-x2 -0.01100292 -0.4024061 0.3804003 0.9998627
## x4-x3  0.24908056 -0.1423226 0.6404837 0.3560447

ANOVA 2f

par(mfrow = c(1, 2))
plot(TukeyHSD(res4))

ANOVA (design non équilibré)

… quand il n’y a pas le même nombre d’observations pour toutes les modalités !

ANOVA

df4b <- data.frame(
  myData = c(
    rnorm(240), 
    rnorm(240, mean = 2), 
    rnorm(240, mean = 5)),
  myGp1 = rep(c("t1", "t2", "t3"), each = 240),
  myGp2 = rep(c("x1", "x2", "x3", "x4"), times = 180)
)
sampleRow <- sample(1:nrow(df4b), size = 701)
df5 <- df4b[sampleRow, ]

ANOVA 2f

table(df5$myGp1, df5$myGp2)

##     
##      x1 x2 x3 x4
##   t1 59 58 59 58
##   t2 59 60 57 59
##   t3 57 59 57 59

ANOVA et les types I, II, et III

En bref : Il y a différentes façons de calculer les sommes des carrés (SS). Il existe trois approches que l’on appelle I, II et III. si le design est non équilibré, il faut utiliser le type III alors que R utilise le type I par default. Le package car a une fonction prête à l’emploi pour faire cela.

ANOVA : fonction Anova du package “car”

pkgCheck <- function(x){ 
    if (!require(x, character.only = TRUE)){
        install.packages(x, dependencies = TRUE)
        if(!require(x, character.only = TRUE)) {
            stop()
        }
    }
}
pkgCheck("car")

## Loading required package: car

## Warning: package 'car' was built under R version 3.6.1

## Loading required package: carData

ANOVA et les types I, II, et III

myAnova <- aov(df5$myData ~ df5$myGp1 * df5$myGp2)
Anova(myAnova, type = "III")

## Anova Table (Type III tests)
## 
## Response: df5$myData
##                     Sum Sq  Df  F value Pr(>F)    
## (Intercept)           0.05   1   0.0442 0.8335    
## df5$myGp1           745.25   2 330.0907 <2e-16 ***
## df5$myGp2             0.23   3   0.0691 0.9764    
## df5$myGp1:df5$myGp2   3.21   6   0.4740 0.8279    
## Residuals           777.78 689                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANOVA et Tukey

dfx <- data.frame(
  myData = c(
    rnorm(60), 
    rnorm(60, mean = 1), 
    rnorm(60, mean = 2), 
    rnorm(60, mean = 3),
    rnorm(60, mean = 4), 
    rnorm(60, mean = 5)),
  myGp1 = rep(c(
    paste0("t0", 1:9), 
    paste0("t", 10:12)), each = 30)
)

ANOVA et Tukey

boxplot(dfx$myData ~ dfx$myGp1)

ANOVA et Tukey

res6 <- aov(dfx$myData ~ dfx$myGp1)
summary(res6)

##              Df Sum Sq Mean Sq F value Pr(>F)    
## dfx$myGp1    11 1030.1   93.65    92.6 <2e-16 ***
## Residuals   348  351.9    1.01                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANOVA et Tukey

TukeyHSD(res6)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = dfx$myData ~ dfx$myGp1)
## 
## $`dfx$myGp1`
##                diff          lwr       upr     p adj
## t02-t01 -0.16848654 -1.022889955 0.6859169 0.9999628
## t03-t01  0.88054403  0.026140621 1.7349474 0.0367186
## t04-t01  0.79325256 -0.061150852 1.6476560 0.0976784
## t05-t01  1.64335011  0.788946697 2.4977535 0.0000000
## t06-t01  1.73215061  0.877747199 2.5865540 0.0000000
## t07-t01  2.84750503  1.993101614 3.7019084 0.0000000
## t08-t01  2.80697791  1.952574494 3.6613813 0.0000000
## t09-t01  3.90891647  3.054513053 4.7633199 0.0000000
## t10-t01  3.66854689  2.814143479 4.5229503 0.0000000
## t11-t01  4.85723744  4.002834028 5.7116409 0.0000000
## t12-t01  4.82186201  3.967458596 5.6762654 0.0000000
## t03-t02  1.04903058  0.194627164 1.9034340 0.0037493
## t04-t02  0.96173910  0.107335690 1.8161425 0.0129900
## t05-t02  1.81183665  0.957433239 2.6662401 0.0000000
## t06-t02  1.90063715  1.046233742 2.7550406 0.0000000
## t07-t02  3.01599157  2.161588156 3.8703950 0.0000000
## t08-t02  2.97546445  2.121061037 3.8298679 0.0000000
## t09-t02  4.07740301  3.222999596 4.9318064 0.0000000
## t10-t02  3.83703343  2.982630021 4.6914368 0.0000000
## t11-t02  5.02572398  4.171320570 5.8801274 0.0000000
## t12-t02  4.99034855  4.135945138 5.8447520 0.0000000
## t04-t03 -0.08729147 -0.941694886 0.7671119 1.0000000
## t05-t03  0.76280608 -0.091597337 1.6172095 0.1324546
## t06-t03  0.85160658 -0.002796835 1.7060100 0.0516390
## t07-t03  1.96696099  1.112557580 2.8213644 0.0000000
## t08-t03  1.92643387  1.072030460 2.7808373 0.0000000
## t09-t03  3.02837243  2.173969020 3.8827758 0.0000000
## t10-t03  2.78800286  1.933599445 3.6424063 0.0000000
## t11-t03  3.97669341  3.122289994 4.8310968 0.0000000
## t12-t03  3.94131797  3.086914562 4.7957214 0.0000000
## t05-t04  0.85009755 -0.004305864 1.7045010 0.0525422
## t06-t04  0.93889805  0.084494639 1.7933015 0.0176078
## t07-t04  2.05425247  1.199849053 2.9086559 0.0000000
## t08-t04  2.01372535  1.159321934 2.8681288 0.0000000
## t09-t04  3.11566391  2.261260493 3.9700673 0.0000000
## t10-t04  2.87529433  2.020890918 3.7296977 0.0000000
## t11-t04  4.06398488  3.209581467 4.9183883 0.0000000
## t12-t04  4.02860945  3.174206035 4.8830129 0.0000000
## t06-t05  0.08880050 -0.765602910 0.9432039 1.0000000
## t07-t05  1.20415492  0.349751504 2.0585583 0.0003092
## t08-t05  1.16362780  0.309224385 2.0180312 0.0006136
## t09-t05  2.26556636  1.411162944 3.1199698 0.0000000
## t10-t05  2.02519678  1.170793369 2.8796002 0.0000000
## t11-t05  3.21388733  2.359483918 4.0682907 0.0000000
## t12-t05  3.17851190  2.324108486 4.0329153 0.0000000
## t07-t06  1.11535441  0.260951002 1.9697578 0.0013465
## t08-t06  1.07482730  0.220423883 1.9292307 0.0025375
## t09-t06  2.17676585  1.322362442 3.0311693 0.0000000
## t10-t06  1.93639628  1.081992867 2.7907997 0.0000000
## t11-t06  3.12508683  2.270683416 3.9794902 0.0000000
## t12-t06  3.08971140  2.235307984 3.9441148 0.0000000
## t08-t07 -0.04052712 -0.894930532 0.8138763 1.0000000
## t09-t07  1.06141144  0.207008027 1.9158149 0.0031127
## t10-t07  0.82104187 -0.033361547 1.6754453 0.0727277
## t11-t07  2.00973241  1.155329001 2.8641358 0.0000000
## t12-t07  1.97435698  1.119953569 2.8287604 0.0000000
## t09-t08  1.10193856  0.247535147 1.9563420 0.0016653
## t10-t08  0.86156898  0.007165572 1.7159724 0.0460024
## t11-t08  2.05025953  1.195856121 2.9046629 0.0000000
## t12-t08  2.01488410  1.160480689 2.8692875 0.0000000
## t10-t09 -0.24036957 -1.094772987 0.6140338 0.9988276
## t11-t09  0.94832097  0.093917562 1.8027244 0.0155484
## t12-t09  0.91294554  0.058542130 1.7673490 0.0246018
## t11-t10  1.18869055  0.334287136 2.0430940 0.0004027
## t12-t10  1.15331512  0.298911704 2.0077185 0.0007278
## t12-t11 -0.03537543 -0.889778845 0.8190280 1.0000000

ANOVA et Tukey

# adapted from https://www.r-graph-gallery.com/84-tukey-test.html
pkgCheck <- function(x){ 
    if (!require(x, character.only = TRUE)){
        install.packages(x, dependencies = TRUE)
        if(!require(x, character.only = TRUE)) {
            stop()
        }
    }
}
pkgCheck("multcompView")
tuk <- TukeyHSD(res6)
generate_label_df <- function(TUKEY, variable){
  Tukey.levels <- TUKEY[[variable]][,4]
  Tukey.labels <- data.frame(multcompLetters(Tukey.levels)['Letters'])
  Tukey.labels$treatment=rownames(Tukey.labels)
  Tukey.labels=Tukey.labels[order(Tukey.labels$treatment) , ]
  return(Tukey.labels)
}
LABELS <- generate_label_df(TUKEY = tuk, "dfx$myGp1")
myCol <- sample(colors(), size = 50)
boxA <- boxplot(dfx$myData ~ dfx$myGp1, 
  ylim=c(min(dfx$myData), 1.1*max(dfx$myData)), 
  col = myCol[as.numeric(LABELS[,1])])
over <- 0.1 * max(boxA$stats[nrow(boxA$stats),])
text(
  x = c(1:nlevels(dfx$myGp1)), 
  y = boxA$stats[nrow(boxA$stats),] + over, 
  labels = LABELS[,1], 
  col = 1)

ANOVA et Tukey

## Loading required package: multcompView

ANOVA et Tukey

# adapted from https://www.r-graph-gallery.com/84-tukey-test.html
tuk <- TukeyHSD(res6)
generate_label_df <- function(TUKEY, variable){
  Tukey.levels <- TUKEY[[variable]][,4]
  Tukey.labels <- data.frame(multcompLetters(Tukey.levels)['Letters'])
  Tukey.labels$treatment=rownames(Tukey.labels)
  Tukey.labels=Tukey.labels[order(Tukey.labels$treatment) , ]
  return(Tukey.labels)
}
LABELS <- generate_label_df(TUKEY = tuk, "dfx$myGp1")
myCol <- sample(colors(), size = 50)
boxA <- boxplot(dfx$myData ~ dfx$myGp1, 
  ylim=c(min(dfx$myData), 1.1*max(dfx$myData)), 
  col = myCol[as.numeric(LABELS[,1])], 
  axes = FALSE, xlab = "", ylab = "")
grid()
axis(1, at = 1:12, labels = levels(dfx$myGp1))
axis(2)
over <- 0.1 * max(boxA$stats[nrow(boxA$stats),])
text(
  x = c(1:nlevels(dfx$myGp1)), 
  y = boxA$stats[nrow(boxA$stats),] + over, 
  labels = LABELS[,1], 
  col = 1)

ANOVA et Tukey

Lorsque la variable à expliquer ne suit pas une loi Normale, nous pouvons transformer les données. Les transformations utilisent le plus souvent le log ou la puissance.

InsecrSprays : The counts of insects in agricultural experimental units treated with different insecticides.

data(InsectSprays)
head(InsectSprays)

##   count spray
## 1    10     A
## 2     7     A
## 3    20     A
## 4    14     A
## 5    14     A
## 6    12     A

mod01 <- aov(sqrt(InsectSprays$count) ~ InsectSprays$spray)
summary(mod01)

##                    Df Sum Sq Mean Sq F value Pr(>F)    
## InsectSprays$spray  5  88.44  17.688    44.8 <2e-16 ***
## Residuals          66  26.06   0.395                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

shapiro.test(resid(mod01)) # normalité résidus

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(mod01)
## W = 0.98721, p-value = 0.6814

bartlett.test( # homogénéité de la variance
  sqrt(InsectSprays$count) ~ InsectSprays$spray) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  sqrt(InsectSprays$count) by InsectSprays$spray
## Bartlett's K-squared = 3.7525, df = 5, p-value = 0.5856

pkgCheck("lmtest")

## Loading required package: lmtest

## Warning: package 'lmtest' was built under R version 3.6.1

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

dwtest(mod01) # indépendance des résidus

## 
##  Durbin-Watson test
## 
## data:  mod01
## DW = 1.8804, p-value = 0.1325
## alternative hypothesis: true autocorrelation is greater than 0

Autres cas

warpbreaks : This data set gives the number of warp breaks per loom, where a loom corresponds to a fixed length of yarn.

data(warpbreaks)
head(warpbreaks)

##   breaks wool tension
## 1     26    A       L
## 2     30    A       L
## 3     54    A       L
## 4     25    A       L
## 5     70    A       L
## 6     52    A       L

mod02 <- aov(warpbreaks$breaks ∼ warpbreaks$wool * warpbreaks$tension)
summary(mod02)

##                                    Df Sum Sq Mean Sq F value   Pr(>F)    
## warpbreaks$wool                     1    451   450.7   3.765 0.058213 .  
## warpbreaks$tension                  2   2034  1017.1   8.498 0.000693 ***
## warpbreaks$wool:warpbreaks$tension  2   1003   501.4   4.189 0.021044 *  
## Residuals                          48   5745   119.7                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

bartlett.test(warpbreaks$breaks, warpbreaks$wool)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  warpbreaks$breaks and warpbreaks$wool
## Bartlett's K-squared = 6.9314, df = 1, p-value = 0.00847

bartlett.test(warpbreaks$breaks, warpbreaks$tension)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  warpbreaks$breaks and warpbreaks$tension
## Bartlett's K-squared = 9.6571, df = 2, p-value = 0.007998

mod03 <- aov(log(warpbreaks$breaks) ∼ 
  warpbreaks$wool * warpbreaks$tension)
summary(mod03)

##                                    Df Sum Sq Mean Sq F value  Pr(>F)   
## warpbreaks$wool                     1  0.313  0.3125   2.234 0.14151   
## warpbreaks$tension                  2  2.176  1.0881   7.779 0.00118 **
## warpbreaks$wool:warpbreaks$tension  2  0.913  0.4566   3.264 0.04686 * 
## Residuals                          48  6.714  0.1399                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

shapiro.test(resid(mod03))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(mod03)
## W = 0.97292, p-value = 0.2583

bartlett.test(log(warpbreaks$breaks), warpbreaks$wool)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  log(warpbreaks$breaks) and warpbreaks$wool
## Bartlett's K-squared = 2.3674, df = 1, p-value = 0.1239

bartlett.test(log(warpbreaks$breaks), warpbreaks$tension)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  log(warpbreaks$breaks) and warpbreaks$tension
## Bartlett's K-squared = 0.94898, df = 2, p-value = 0.6222