One-way Analysis of VarianceIn this section, we will discuss how can we assess whether the population means of a continuous variable are different for different groups.

One-way ANOVA AssumptionsLet \[Y_{ij} = \mu_i + \varepsilon_{ij}\] denote the value of the response variable for the #j#th replicate randomly-assigned to treatment group #i# from the population. Thus #Y_{ij}# has mean #\mu_i# and randomly deviates from this mean by #\varepsilon_{ij}#.

As stated previously, we assume that the error \[\varepsilon_{ij} \sim N\Big(0,\sigma^2\Big)\] in this model, where #\sigma^2 > 0# does not depend on #i#.

This is the same as assuming: \[Y_{ij} \sim N\Big(\mu_i,\sigma^2 \Big)\] Consequently, for #i = 1,...,I#, the variance of #\bar{Y}_i# is #\cfrac{\sigma^2}{n_i}#.

The assumption that #\sigma^2# is the same for all the treatment groups is called constant variance or homoscedasticity.

This assumption can be checked by checking the residuals #Y_{ij} - \bar{Y}_i# for #i = 1,...,I# and #j = 1,...,n_i# for each #i#. This is because the residuals should represent a random sample from the #N\Big(0,\sigma^2\Big)# if the assumptions are true.

A scatterplot of these residuals against the corresponding fitted values #\bar{Y}_1,...\bar{Y}_I# should show approximately the same spread for each fitted value, and no extreme outliers, as seen in the following plot:

Furthermore, for each #i# a Normal quantile-quantile plot should also be used to confirm the assumption of normality of #Y_{ij}# for each level #i#. The following plot casts doubt on this assumption for group 1 and group 2:

Using the above plots to assess the validity of these two assumptions is a qualitative approach. We can formally assess the validity of the two assumptions using various hypothesis tests. The Shapiro-Wilk test is a commonly-used test for the normality assumption, while Bartlett's test is commonly used to test the constant variance assumption. For both tests we do not want to reject the null hypothesis.

One-Way ANOVA in R#\mathrm{R}# can be used to efficiently conduct analysis of variance (ANOVA). We demonstrate this using the following data, which you can copy and paste into #\mathrm{R}#:

> Y = c(24.6,22.8,27.5,21.3,24.0,33.4,25.2,19.6,29.4,20.8,32.7,23.3,18.2,28.7,24.0)

> Group = c(2,1,3,1,2,3,2,1,3,1,3,2,1,3,2)

> Data = data.frame(Y, Group)

Here #\mathtt{Y}# is the response variable, and #\mathtt{Group}# is the categorical factor, which has 3 levels.

There are #15# subjects, and the #k#th position in the two vectors gives the group level and the response variable measurement for the #k#th subject, for #k = 1,...,15#. We merge these together into a data frame which we aptly name #\mathtt{Data}#. (If you read data into #\mathrm{R}# from a file using the #\mathtt{read.table()}# function, it will be loaded as a data frame. You can give it any name you like, though Data usually makes the most sense.)

We have a one-factor experiment, with #\mathtt{Y}# as the response variable and #\mathtt{Group}# as the factor. We can create a model in #\mathrm{R}# using the #\mathtt{aov()}# function and then inspect the results using the #\mathtt{summary()}# function as follows:

> Model = aov(Y ~ factor(Group), data=Data)

> summary(Model) 

We stored the output from the #\mathtt{aov()}# function in an object that we named #\mathtt{Model}#. You can choose any name you like.

Note that we have to convert the numeric vector #\mathtt{“Group”}# into the factor class in #\mathrm{R}# using the #\mathtt{factor()}# function if we want the #\mathtt{aov() }# function to work properly. Even if #\mathtt{Group}# was already in the factor class, using the #\mathtt{factor()}# function would do no harm, so to be on the safe side it is a good idea to always do this.

Note also that the #\mathtt{aov()}# function allows us to specify the data frame in which to find the variables, so that we don't have to type #\mathtt{Data$Y}# and #\mathtt{Data$Group}#. If the variables are not nested inside a data frame, but are in the #\mathrm{R}# workspace, then we would not use the #\mathtt{“data=”}# setting.

The output of the #\mathtt{summary()}# function is an ANOVA table (which we can also obtain using #\mathtt{anova(Model)}#), which allows us to see the treatment sum of squares (in the first row, after #\mathtt{“factor(Group)”}#) and the error sum of squares (in the second row, after #\mathtt{“Residuals”}#), preceded by their respective degrees of freedom and followed by their mean squares, the F-statistic and the P-value for the F-test of the null hypothesis #H_0: \mu_1 = \mu_2 = \mu_3# against the alternative hypothesis #H_1: \mu_i \neq \mu_j# for at least one pair #i# and #j#:

Signif.codes: 0'***' 0.001'**' 0.01'*' 0.05'.' 0.1'' 1

The #\mathtt{“Signif.codes”}# let you know how small the P-value of the F-test is, in case you can't figure that out yourself.

Here the given P-value translates as #0.00000838#, which is between #0# and #0.001#, as signified by #3# stars. So we would reject the null hypothesis and conclude that there are at least two population means which differ out of the three groups. However, this conclusion is only valid if the assumptions are satisfied.

To qualitatively check the normality of #\mathtt{Y}# in each group, we can use a single Q-Q-plot of the model residuals:

> qqnorm(Model$residuals, pch=20)

> qqline(Model$residuals)

We get:

The points are close enough to the line that we would not reject the assumption of normality.

There is also a hypothesis test, called the Shapiro-Wilk test, available in #\mathtt{R}#, whereby you can formally test #H_0#: the residuals come from a normal distribution, against #H_1#: the residuals do not come from a normal distribution. In #\mathtt{R}#, use:

> shapiro.test(Model$residuals)

A small P-value implies that the residuals do not come from a normal distribution, which in turn implies that the dependent variable is not normally-distributed in each group of the factor. This violates the assumption and makes any inference based on that assumption invalid.

To qualitatively check the assumption of constant variance (homoscedasticity), we can make a scatterplot of the model residuals against the fitted values (i.e., the group sample means):

> plot(Model$fitted.values, Model$residuals, pch=20, xlab="Fitted values", ylab="Residuals")

The spread of points in the middle group (group 2) is smaller than the spread of points in the other two groups, which indicates we may have a problem with the constant variance assumption. But this is hard to address visually.

There is also a hypothesis test, called Barlett's test, available in #\mathtt{R}#, whereby you can formally test #H_0#: the variance is the same in all groups of the factor, against #H_1#: the variance is not the same in all groups of the factor. In #\mathtt{R}#, this test can be performed for the above data as follows:

> bartlett.test(Y~factor(Group), data=Data) 

A small P-value implies that the constant variance assumption is violated, and thereby invalidates any inference based on this assumption.

For the data, if you need the group sample means, along with the overall mean (the "grand mean"), use the #\mathtt{model.tables()}# function, with the setting #\mathtt{type=“means”}#:

> model.tables(Model, type="means")

However, if you want confidence intervals for the population means, you have to first fit a new model with a small adjustment: you have to put '#\mathtt{0 +}#' after the '#\mathtt{\sim}#' symbol.

Then you can use the #\mathtt{confint()}# function on this adjusted model to obtain confidence intervals for the population means in each group, at any specified confidence level (default is 95%):

> Model2 = aov(Y ~ 0 + factor(Group), data=Data)

> confint(Model2, level=0.99)

We get the output:

We can see that the #99\%# CI for the mean of group 3 has no overlap with the #99\%# CIs for the means of the other two groups, so it seems very clear that group 3 has a larger mean than groups 1 and 2.