The Assumptions of One-Way ANOVA

1. The response variable is normally distributed in each group
2. The response variable has equal variance across the groups
3. The data are drawn on the basis of random sampling

The Null and the Alternative Hypotheses

H0: There is no difference between the group means (μ1=μ2=μ3)
HA: There is a difference between the group means (at least one μi≠μj)

The Testing Sequence

As usual we will test the assumptions of (a) normality in each group (with the shapiro.test) and then (b) equality of variance across all groups (with Levene’s test for homogeneity of variances).

However, since ANOVA is fairly robust to violations of these assumptions we won’t abandon it unless we absolutely have to. In a practical research situation you will have to go by standards accepted by your field/sub-field so if few folks opt for a non-parametric test then you should be fine doing so as well. However, it never hurts to have run a non-parametric test just in case someone criticizes your work or asks why you ran the ANOVA when the assumptions were clearly violated. My suggestion would be to point to the results of the non-parametric test which, hopefully, lead to the same conclusion as that generated by ANOVA. Remember …

1. If normality is not rejected but we have unequal variances then Welch’s one-way ANOVA is a good choice.

W∗=∑wj(ˉyj−ˆμ)2/(k−1)1+[2(k−2)/(k2−1)]∑hj

where wj=njs2j; ˆμ=∑wjyjW; W=∑wj; hj=(1−wj/W2)(nj−1), and; f=(k2−1)3∑hj

W∗∼F(df1=k−1,df2=f)

2. If both normality and equal variance are rejected and no transformation works, then we may have to settle for the weaker test – the Kruskal-Wallis test.

H=[12N(N+1)×∑Tc2nc]−3(N+1)

where N is the total number of scores in the study; Tc is the rank total for each group; nc is the number of units in each group

H∼χ2with df=k−1

Welch’s One-Way ANOVA

Welch’s one-way ANOVA is used when normality holds but homogeneity of variance does not hold. One may be tempted to attempt a transformation in the hope that this will ensure equal variances but this seldom works.

I’ll generate some fake data investigating the effects of grip-types on cross-country skiing performance.

score = c(164.2, 162.4, 163.2, 169.7, 171.2, 172.9, 140.1, 164.7, 172.9, 151.3)
grips = c(rep(1, 3), rep(2, 3), rep(3, 4))
grips = factor(grips, levels = c(1, 2, 3), labels = c("Classic", "Integrated", "Modern"))
ski = cbind.data.frame(score, grips)

Let us do a quick check for normality and equal variances.

library(car)
with(ski, boxplot(score ~ grips, horizontal = TRUE))

with(ski, tapply(score, grips, FUN = shapiro.test))

## $Classic
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.9959, p-value = 0.8777
## 
## 
## $Integrated
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.9987, p-value = 0.9311
## 
## 
## $Modern
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.97406, p-value = 0.8665

with(ski, leveneTest(score, grips, center = median))

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value   Pr(>F)   
## group  2  9.5746 0.009925 **
##        7                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Quite clearly Shapiro-Wilk indicates that normality holds but Levene’s test soundly rejects the NULL of equal variances. This would be the ideal situation for Welch’s one-way ANOVA.

oneway.test(score ~ grips, data = ski, var.equal = FALSE)

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  score and grips
## F = 24.975, num df = 2.0000, denom df = 4.1069, p-value = 0.005029

The p-value is 0.005 and so we can reject H0; the data indicate that mean scores do differ for at least one grip-type.

This then begs the question: How do we determine which groups are really different? The Games-Howell test is commonly used in these situations because it is designed for situations where (i) the sample sizes differ across the groups, and/or (ii) the variances are unequal. The test uses the sample size and variance of each group being compared in a specific pairwise comparison.

library(userfriendlyscience)
posthocTGH(y = ski$score, x = ski$grips)

##            n means variances
## Classic    3   163      0.81
## Integrated 3   171      2.56
## Modern     4   157    209.98
## 
##                    diff ci.lo ci.hi    t  df   p
## Integrated-Classic    8   3.7    12 7.54 3.2 .01
## Modern-Classic       -6 -36.2    24 0.83 3.0 .71
## Modern-Integrated   -14 -43.9    16 1.92 3.1 .27

The comparisons indicate a significant difference only for Classic vs. Integrated, with a p-value of 0.0084.

Non-Parametric Tests: The Kruskal-Wallis

The amount of oxygen obtained in each breath at high altitude can be as low as one-third that obtained at sea level. Studies have begun to examine whether indigenous people living at high elevations have physiological attributes that compensate for the reduced availability of oxygen. A reasonable expectation is that they should have more hemoglobin, the molecule that binds and transports oxygen in the blood. To test this, researchers sampled blood from males in three high-altitude human populations: the high Andes, high-elevation Ethiopia, and Tibet, along with a sea-level population from the USA (Bell et al. 2002). These data are given below:

http://whitlockschluter.zoology.ubc.ca/wp-content/data/chapter02/chap02e3cHumanHemoglobinElevation.csv

Do mean hemoglobin levels differ between these groups?

We start with the usual checks for normality.

hemo = read.csv("http://whitlockschluter.zoology.ubc.ca/wp-content/data/chapter02/chap02e3cHumanHemoglobinElevation.csv")
summary(hemo)

##              id         hemoglobin       population  
##  Andean.male1 :   1   Min.   :10.40   Andes   :  71  
##  Andean.male10:   1   1st Qu.:14.80   Ethiopia: 128  
##  Andean.male11:   1   Median :15.45   Tibet   :  59  
##  Andean.male12:   1   Mean   :15.56   USA     :1704  
##  Andean.male13:   1   3rd Qu.:16.15                  
##  Andean.male14:   1   Max.   :25.06                  
##  (Other)      :1956

boxplot(hemo$hemoglobin ~ hemo$population, horizontal = TRUE)

library(lattice)
qqmath(~ hemoglobin | population, data = hemo,
       prepanel = prepanel.qqmathline,
       panel = function(x, ...) {
          panel.qqmathline(x, ...)
          panel.qqmath(x, ...)
       })

library(moments)
tapply(hemo$hemoglobin, hemo$population, FUN = skewness)

##      Andes   Ethiopia      Tibet        USA 
##  0.4230040  0.0437093  0.6690375 -0.2542965

tapply(hemo$hemoglobin, hemo$population, FUN = shapiro.test)

## $Andes
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.97937, p-value = 0.293
## 
## 
## $Ethiopia
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.995, p-value = 0.9359
## 
## 
## $Tibet
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.9551, p-value = 0.02922
## 
## 
## $USA
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.99235, p-value = 9.567e-08

library(car)
leveneTest(hemo$hemoglobin, hemo$population, center = median)

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value    Pr(>F)    
## group    3  30.114 < 2.2e-16 ***
##       1958                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Neither normality nor homogeneity of variances holds across the four groups. Might a transformation work?

hemo$log.hemo = log(hemo$hemoglobin)

boxplot(hemo$log.hemo ~ hemo$population, horizontal = TRUE)

library(lattice)
qqmath(~ log.hemo | population, data = hemo,
       prepanel = prepanel.qqmathline,
       panel = function(x, ...) {
          panel.qqmathline(x, ...)
          panel.qqmath(x, ...)
       })

library(moments)
tapply(hemo$log.hemo, hemo$population, FUN = skewness)

##       Andes    Ethiopia       Tibet         USA 
##  0.06522668 -0.15719251  0.35690592 -0.55585979

tapply(hemo$log.hemo, hemo$population, FUN = shapiro.test)

## $Andes
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.9873, p-value = 0.6952
## 
## 
## $Ethiopia
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.99394, p-value = 0.8618
## 
## 
## $Tibet
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.97191, p-value = 0.1888
## 
## 
## $USA
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.98154, p-value = 5.211e-14

library(car)
leveneTest(hemo$log.hemo, hemo$population, center = median)

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value    Pr(>F)    
## group    3  15.136 9.706e-10 ***
##       1958                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

hemo$sqrt.hemo = sqrt(hemo$hemoglobin)

boxplot(hemo$sqrt.hemo ~ hemo$population, horizontal = TRUE)

library(lattice)
qqmath(~ sqrt.hemo | population, data = hemo,
       prepanel = prepanel.qqmathline,
       panel = function(x, ...) {
          panel.qqmathline(x, ...)
          panel.qqmath(x, ...)
       })

library(moments)
tapply(hemo$sqrt.hemo, hemo$population, FUN = skewness)

##       Andes    Ethiopia       Tibet         USA 
##  0.24327995 -0.05585916  0.51125721 -0.39930875

tapply(hemo$sqrt.hemo, hemo$population, FUN = shapiro.test)

## $Andes
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.98481, p-value = 0.5481
## 
## 
## $Ethiopia
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.99511, p-value = 0.9418
## 
## 
## $Tibet
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.96472, p-value = 0.08487
## 
## 
## $USA
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.98789, p-value = 9.55e-11

library(car)
leveneTest(hemo$sqrt.hemo, hemo$population, center = median)

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value    Pr(>F)    
## group    3  21.402 1.222e-13 ***
##       1958                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

hemo$sq.hemo = hemo$hemoglobin^2

boxplot(hemo$sq.hemo ~ hemo$population, horizontal = TRUE)

library(lattice)
qqmath(~ sq.hemo | population, data = hemo,
       prepanel = prepanel.qqmathline,
       panel = function(x, ...) {
          panel.qqmathline(x, ...)
          panel.qqmath(x, ...)
       })

library(moments)
tapply(hemo$sq.hemo, hemo$population, FUN = skewness)

##      Andes   Ethiopia      Tibet        USA 
## 0.78917353 0.23857087 0.99256569 0.01031461

tapply(hemo$sq.hemo, hemo$population, FUN = shapiro.test)

## $Andes
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.95991, p-value = 0.02324
## 
## 
## $Ethiopia
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.99109, p-value = 0.5886
## 
## 
## $Tibet
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92881, p-value = 0.001958
## 
## 
## $USA
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.99612, p-value = 0.0002425

library(car)
leveneTest(hemo$sq.hemo, hemo$population, center = median)

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value    Pr(>F)    
## group    3  56.011 < 2.2e-16 ***
##       1958                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

kruskal.test(hemo$hemoglobin ~ hemo$population)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  hemo$hemoglobin by hemo$population
## Kruskal-Wallis chi-squared = 199.56, df = 3, p-value < 2.2e-16