One-Way ANOVA

By way of a broad overview, the approach we will take to this problem is to partition the total variability into two parts: one attributable to random variability within each group/condition/treatment, and one attributable to the variability between the groups themselves. The ANOVA compares the ratio of the variability between groups to the variability within the groups—if the variability between groups is larger than the random variability within groups, then the grouping variable is said to have a significant effect on the dependent variable. Put another way, if most of the total variability (e.g., in plant height) is due to the different groups, then the grouping factor (e.g., fertilizer) is affecting the outcome more at a greater-than-chance level.

Let’s adopt some useful notation (for observations, grand mean, group means, etc) before we return to our plant-growth scenario.

For subjects (individuals, plants, students, “units”), we use the subscript \(i; \ (i=1,2,3,...,i,...,n_k)\) where \(n_k\) is the total number of observations in a given group \(k\).
For groups (conditions, treatments, etc.), we therefore use the subscript \(k; \ (k=1,2,3,...,k,...,K)\) where \(K\) is the total number of groups.

\(x_{ik}\) : score for subject i in group k
\(\bar{x}_{\bullet k}\) : mean of group k
\(\bar{x}_{\bullet \bullet}\): total (or “grand”) mean

In our example data, the number of groups is \(K=4\) and the number of observations in each group is the same (\(n_1=n_2=...n_K=5\)). The total number of observations across groups is given by \(N=\sum_{k=1}^Kn_k=n_1+n_2+...n_K\)

Visualizing Sums of Squares

These formulas may be daunting at first, and it can be very helpful to see an illustration of what each one represents. So without further ado,

This plot is really important! It allows you to visualize the dispersion of the data in each of the three SS calculations. On the left, we have the Total deviations \((x_{ik}-\bar{x}_{\bullet \bullet})\): black dotted lines represent the deviations of each observation (large dots) from the grand mean (black bar); in the middle, we have Within-Group deviations \((x_{ik}-\bar{x}_{\bullet k})\): the deviations of each observation from their respective group mean (colored lines); on the right, we have Between-Group deviations \((x_{\bullet k}-\bar{x}_{\bullet \bullet})\): the deviations of each group mean from the grand mean times the number of observations in each group. To calculate the Total, Within-Group, and Between-Group Sums of Squares, you simply square each of the deviations in their respective plots above and then add them up, giving you three different sums.

Sums of Squares Plot

OK great, so we have our sums of squares: How do we perform the ANOVA with them? ANOVA uses an F-test for a ratio of two variances. Recall that we can think of our sums of squares as the numerator of a variance… what about the denominator? It will be helpful to consult a generic ANOVA output table:

Degrees of Freedom for \(SS_W\) and \(SS_B\)

The degrees of freedom for \(SS_W\) is easy to remember too: just look at the components of the sum in the formula and note that the observations are constrained by the 4 group means, \(\bar{x}_{\bullet k}\). Thus, 20 observations go in, but once we know 4 of the observations in a given group, we automatically know the fifth one (because we have the group mean, and \(x_5=n_k\bar{x}_{\bullet k}-(x_1+x_2+x_3+x_4)\)). Thus, the degrees of freedom is given by \(N-K = 20-4 = 16\), the total number of observations minus the total number of groups.

Lastly, to find the degrees of freedom for \(SS_B\) we look at the components of the sum-of-squares calculation: here the unique observations are the group means themselves, and they are constrained by the grand mean, such that \[ n_5\bar{x}_{\bullet 5}=N\bar{x}_{\bullet \bullet}-(n_1\bar{x}_{\bullet 1}+n_2\bar{x}_{\bullet 2}+n_3\bar{x}_{\bullet 3}+n_4\bar{x}_{\bullet 4}) \] Thus, the degrees of freedom for \(SS_B\) are found by taking the number of group means and subtracting 1 for the grand mean: \(K-1=4-1=3\).

The F-Test

We saw above that the Mean Squares are found by dividing the Sums of Squares by their DF, \(MS=\frac{SS}{df}\). Let’s fill in the ANOVA table now with all of our calculations.

Our F-statistic of interest is the ratio of the between-group variability (\(MS_B\)) to the within-group variability (\(MS_W\)), \(F=\frac{MS_B}{MS_W}\). (Note that any Mean Square statistic follows a chi-squared distribution, so the F-distribution gives the distribution of the ratio of two chi-squared random variables.) To see if our statistic is unusually large under the null hypothesis of no group differences (that is, to see if the between-group variance is larger than expected relative to the random within-group noise), we look at the probability of observing an F-statistic at least as large as this one under the null hypothesis. Note that the F distribution depends on two parameters: the numerator degrees of freedom (\(K-1=3\)) and the denominator degrees of freedom (\(N-K=16\)), so our test statistic follows the distribution \(F_{(3,16)}\):

F.stat<-(100/12)
numer.df<-3
denom.df<-16

x1<-seq(F.stat,10,len=100)
x2<-seq(qf(.95,numer.df,denom.df),10,len=100)
y1<-df(x1,numer.df,denom.df)
y2<-df(x2,numer.df,denom.df)
{curve(df(x,numer.df,denom.df),main="",xlim=c(0,10),xlab="F")
polygon(c(x2[1],x2,x2[100]),c(0,y2,0),col="black",border=NA)
polygon(c(x1[1],x1,x1[100]),c(0,y1,0),col="grey",border=NA)
abline(v=F.stat,lty=2)}

#p-value
1-pf(F.stat,numer.df,denom.df)

## [1] 0.0014506

We see in the plot above that the probability of observing an F-statistic of \(8.\bar{3}\) (dotted line) or larger under the null hypothesis is very unlikely, \(p>F=.001\). The area in black under the \(F_{(3,16)}\) marks of 5% of the probability distribution. Anything observed F-statistic in the black region would reject the null hypothesis (that group means do no differ) at the traditional significance level \(\alpha=.05\). Another way of stating the null hypothesis is this: all groups are random samples of the same population, and thus that all treatments result in the same effect. Rejecting the null hypothesis means we have non-negligible evidence that different treatments result in different effects.

Doing a One-Way ANOVA in R

Performing a one-way ANOVA in R is very straightforward. Notice that R gives exactly the same results as our hand calculations. Note also that “Mean Squared Residuals” (or “Mean Squared Error”, MSE) refers to the same thing as “Mean Squared Within”—it is whatever variability is left over after accounting for the effect of the groups.

anova1<-aov(Height~Group,data=dat.long)
summary(anova1)

##             Df Sum Sq Mean Sq F value  Pr(>F)   
## Group        3    100   33.33   8.333 0.00145 **
## Residuals   16     64    4.00                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Independence

This one is pretty easy: a given observation should have no influence on any of the others. There are plenty of situations when this simplifying assumption doesn’t hold (imagine if you measure the same person twice!) and we will deal with these situations in the next post. For now, one plant’s height should give you no information about another plant’s height beyond knowing what fertilizer treatment they were exposed to. We achieved this by randomizing which plants got which treatment: thus, on average, the only thing different about each group should be their treatment exposure.

Normality (of Errors)

One easy way to check normality is to look at the distribution of data in each group: if these distributions are normal, then the errors (or residuals) will be normal. Fine, but what’s a residual? We will show in the next post that ANOVA is a special case of regression, in which the concept of residuals plays a crucial role. ANOVA analyses tend to gloss over residuals, but let’s take a minute to compute them.

A residual is what’s left over after we account for the effect of our independent variables. In our example, once we account for the fertilizer treatment group mean \(\bar{x}_{\bullet k}\), the residual is what’s left over. ANOVA assumes that these Within-Group deviations \((x_{ik}-\bar{x}_{\bullet k})\) are normally distributed. Here’s our data again:

So, going down columns, our residuals \((x_{ik}-\bar{x}_{\bullet k})\) are \((6-8),(7-8),(8-8),...,(10-10),(11-10),(13-10)\).

We can get these straight from our ANOVA object in R; typically, they are examined using a quantile-quantile plot (theoretical vs. observed data quantiles should be similar) or tested formally with the Shapiro-Wilk test (which tends to be overly harsh).

round(anova1$residuals,4)

##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 
## -2 -1  0  1  2 -3 -2  0  2  3 -2  0  0  0  2 -3 -1  0  1  3

{qqnorm(anova1$residuals); qqline(anova1$residuals)}

shapiro.test(anova1$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  anova1$residuals
## W = 0.94361, p-value = 0.2804

In the S-W test, the null hypothesis is that the data are normally distributed: here, we fail to reject that null hypothesis, confirming what we see in our Q-Q plot (the quantiles of a normal distribution look similar to the quantiles of our observed data).

Equal Variance (Homoskedasticity)

Fortunately, ANOVA is very robust to non-normality (provided the groups are roughly symmetrical and, if skewed, skewed the same way) and relatively robust to unequal variances (provided the sample size in each group is about the same). Looking at our box plots, we see that the spreads look a little different across the four groups. However, two standard procedures for testing the null hypothesis that group variances are equal (Bartlett’s test and Levene’s test) both fail to reject this assumption.

library(car)

## Loading required package: carData

boxplot(dat.long$Height~dat.long$Group)

bartlett.test(dat.long$Height~dat.long$Group)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  dat.long$Height by dat.long$Group
## Bartlett's K-squared = 1.6465, df = 3, p-value = 0.6489

leveneTest(dat.long$Height~dat.long$Group)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  1.0256 0.4075
##       16

Tukey’s Test

One way around this is to use an approach variously refered to as Tukey’s test, Tukey’s HSD (honestly significant difference), Tukey’s paired comparisons, etc.

For example, to test whether Group A and Group B differed in Height, we would calculate q \[ q=\frac{|\bar{x}_A-\bar{x}_B|}{\sqrt{MS_W/n}}=\frac{|8-14|}{\sqrt{4/5}}=2\sqrt5=4.472 \] We look this value up in table for the studentized-range distribution to find the p-value (or probability of observing a statistic at least as extreme under the null hypothesis), which can be found in R by ptukey(q,nmeans=K,df=N-K,lower.tail=F)

ptukey(4.472,nmeans=4,df=16,lower.tail=F)

## [1] 0.02779534

We can easily get this for all pairwise comparisons with the TukeyHSD() function, which takes ANOVA-type model as its argument:

TukeyHSD(anova1)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Height ~ Group, data = dat.long)
## 
## $Group
##     diff        lwr       upr     p adj
## B-A    4  0.3810644  7.618936 0.0277901
## C-A   -2 -5.6189356  1.618936 0.4161607
## D-A    2 -1.6189356  5.618936 0.4161607
## C-B   -6 -9.6189356 -2.381064 0.0011369
## D-B   -2 -5.6189356  1.618936 0.4161607
## D-C    4  0.3810644  7.618936 0.0277901

Thus, we see that groups A and B are significantly different, C and B are significantly different, and D and C are significantly different while holding the familywise error rate for all six comparisons at \(\alpha_{fw}=.05\), because unlike t-tests, the Tukey test applies simultaneously to all pairwise comparisons.

If you find yourself with unequal sample sizes across groups, you simply use the harmonic mean of the group sizes \(n_k\), \(\frac{K}{\sum_{k=1}^K \frac{1}{n_k}}=\frac{K}{\frac{1}{n_1}+\frac{1}{n_2}+...\frac{1}{n_K}}\), instead of \(n\) in your calculation of q.

Bonferroni T-Tests

Another approach is to say OK, we are doing 6 tests, so lets set the familywise error rate to \(\alpha_{fw}=\frac{\alpha}{5}=\frac{.05}{6}=.008\bar{3}\). Thus, for any one of our paired comparisons to be different, we would need to observe a \(p<.008\bar{3}\). This method is much more conservative (that is, it tends to over-correct, making it more likely that we make a Type II error by failing to reject a false \(H_0\)). We can see this by running all 6 t-tests like so:

# A vs. B
t.test(dat.long[dat.long$Group=='A',]$Height,dat.long[dat.long$Group=='B',]$Height,var.equal = T)$p.value
## [1] 0.01756221
#
# A vs. C
t.test(dat.long[dat.long$Group=='A',]$Height,dat.long[dat.long$Group=='C',]$Height,var.equal = T)$p.value
## [1] 0.06806525
#
# A vs. D
t.test(dat.long[dat.long$Group=='A',]$Height,dat.long[dat.long$Group=='D',]$Height,var.equal = T)$p.value
## [1] 0.1411133
#
# B vs. C
t.test(dat.long[dat.long$Group=='B',]$Height,dat.long[dat.long$Group=='C',]$Height,var.equal = T)$p.value
## [1] 0.00175135
#
# B vs. D
t.test(dat.long[dat.long$Group=='B',]$Height,dat.long[dat.long$Group=='D',]$Height,var.equal = T)$p.value
## [1] 0.2237477
#
# C vs. D
t.test(dat.long[dat.long$Group=='C',]$Height,dat.long[dat.long$Group=='D',]$Height,var.equal = T)$p.value
## [1] 0.009632847

Notice here that only B vs. C is found to be significantly different at \(\alpha=.008\bar{3}\) whereas using Tukey’s test, A vs. B and C vs. D were significant as well.

Stepdown: Bonferroni-Holm

We will look at one last correction that is a less conservative variation on bonferroni. Here, we run our six t-tests like we did above, and we sort our p-values from least to greatest:

bf_holm<-data.frame(Comparison=c("AB","AC","AD","BC","BD","CD"),Pvalue=c(
t.test(dat.long[dat.long$Group=='A',]$Height,dat.long[dat.long$Group=='B',]$Height,var.equal = T)$p.value,
t.test(dat.long[dat.long$Group=='A',]$Height,dat.long[dat.long$Group=='C',]$Height,var.equal = T)$p.value,
t.test(dat.long[dat.long$Group=='A',]$Height,dat.long[dat.long$Group=='D',]$Height,var.equal = T)$p.value,
t.test(dat.long[dat.long$Group=='B',]$Height,dat.long[dat.long$Group=='C',]$Height,var.equal = T)$p.value,
t.test(dat.long[dat.long$Group=='B',]$Height,dat.long[dat.long$Group=='D',]$Height,var.equal = T)$p.value,
t.test(dat.long[dat.long$Group=='C',]$Height,dat.long[dat.long$Group=='D',]$Height,var.equal = T)$p.value))

bf_holm[order(bf_holm$Pvalue),]

##   Comparison      Pvalue
## 4         BC 0.001751350
## 6         CD 0.009632847
## 1         AB 0.017562209
## 2         AC 0.068065248
## 3         AD 0.141113281
## 5         BD 0.223747695

Now then, the Bonferroni-Holm procedure is as follows

• Rank comparisons from smallest p-value to largest (\(c=1,2,...N_c\))
• Start with the comparison with the smallest p-value (i.e., \(c=1\))
• Set the \(\alpha_c\) for each comparison \(c\) at \(\frac{.05}{(N_c+1)-c}\).
• If \(p_1 \le \alpha_1\), move on to \(c=2\)
• STOP when \(p_c \not\le \alpha_c\) (no more tests)

For example, our smallest p-value is \(BC\) and we have \(\alpha_1=\frac{.05}{(6+1)-1}=\frac{.05}{6}=.008\bar{3}\). So for the first comparison, the p-value is the same as Bonferroni. Since \(.00175 \le .008\bar{3}\), we say that groups \(B\) and \(C\) are significantly different, and we move on to our the comparison with the second-smallest p-value, \(CD\). Here we get \(\alpha_2=\frac{.05}{(6+1)-2}=\frac{.05}{5}=.01\), and since \(p_2=.00963<.01\), groups \(C\) and \(D\) are significantly different. For the comparison with the third-smallest p-value, \(AB\), we have \(\alpha_3=\frac{.05}{(6+1)-3}=\frac{.05}{4}=.0125\), but since \(p_3=.0176 \not\le .0125\), groups \(A\) and \(B\) are not significantly different, and we stop.

Notice that using this procedure, 2 comparisons were found to be significantly different (compared to 1 using Bonferonni and 3 using Tukey).

Also note that the Bonferonni and Bonferroni-Holm procedures and general and can be used for any set of hypothesis tests that one is conducting.

kruskal.test(Height~Group,data=dat.long)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Height by Group
## Kruskal-Wallis chi-squared = 12.291, df = 3, p-value = 0.006451