ANOVA and t-test
Michael Sinclair
May 23, 2019
setwd("~/VCU/Lab/R_Tutorial/AnovaTtest")A Brief Review of ANOVA
ANOVA stands for Analysis of Variance. It tests for differences in means among 2 or more groups, and with two-way ANOVA, you can test two factors. Your response variable should technically be continuous (i.e. numbers where decimals are possible) or integers with a wide possible range. Technically you should not use ANOVA with variables that have a very fixed range (e.g. proportions that range from 0-1). Many people do though, so do what you think is best and what you see supported in the literature.
*You may be interested to learn about the different ‘types’ of ANOVA (I, II, and III). They deal with unequal sample sizes differently. R has a different default than many other statistical programs. For a first glance, check out this website by R-bloggers about ANOVA types I/II/III sums of squares: https://www.r-bloggers.com/anova-%E2%80%93-type-iiiiii-ss-explained/.*
One-way ANOVA
Let’s take a look at the relatively simply coded one-way ANOVA first. We will use the iris dataset, a freely available example data set in R that we can call directly (instead of reading in data from a .csv like we did in Getting Started with R).
We will test if sepal width differs among 3 iris species.
Checking Assumptions
The assumptions of ANOVA are 1. Samples are independent from one another (within and among groups). 2. Samples were randomly chosen. 3. Response variable is normally distributed (often people check residual normality too, but I won’t be showing that here). To do this, google qqplots in R. 4. Variances are equal among groups.
You have to determine if your sampling was independent and random. The other assumptions can be tested with the code below.
Testing Normality
First we will test normality using two methods. The first is by looking at the histogram using the hist() function. This is a good way to ‘eyeball’ the data. The second way is to use a statistical test for normality called the Shapiro-Wilk test. There is some disagreement on if this test does a good job, but it is widely used, so I do too (at least for the moment).
attach(iris) #Attachs data to make it easier to call the data later
hist(Sepal.Width)
shapiro.test(Sepal.Width)##
## Shapiro-Wilk normality test
##
## data: Sepal.Width
## W = 0.98492, p-value = 0.1012The histogram is rather normal, but maybe has a bit of right skew. The Shapiro-Wilk test shows a p-value of 0.10, which is greater than 0.05, suggesting that the data is normally distributed. In summary, if the Shapiro-Wilk p-value is greater than 0.05, it suggests the data is normal; if it’s less than 0.05, then it’s nonnormal.
Testing Variances
Next we will test if the variances are the same among the groups using the Levene’s test.
We will call the package car to run the function leveneTest. If you have not installed the package car before, run the line #install.packages… below once and never again. One installation is all you need per computer. To run it, delete the pound sign and hit run. Then put the hashtag back to make it a ‘comment’ and not code.
#install.packages("car") ##Only run this once ever per computer
library(car) #run this every time you open r and need the package
lm.iris1<-lm(Sepal.Width~Species,data=iris) #makes a 'linear model' object which we will use for the next line of code and the ANOVA code
leveneTest(lm.iris1)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.5902 0.5555
## 147The variances are equal, as the p-value is 0.556 (greater than 0.05).
This dataset satisfies our assumptions, so we can run the ANOVA!
Running the one-way ANOVA
anova(lm.iris1)## Analysis of Variance Table
##
## Response: Sepal.Width
## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 11.345 5.6725 49.16 < 2.2e-16 ***
## Residuals 147 16.962 0.1154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is at least one group that is significantly different (p is super small: p< 2.2 x 10^-16).
Let’s also look at the summary of the ANOVA. It gives you the average of the first species, setosa (estimate of the intercept) and how the other groups are different from setosa (see the estimates of Speciesversicolor and Speciesvirginica). I. ersicolor has, on average, sepals that are 0.658 mm narrower that I. setosa; I. virginica is 0.454 mm narrower than I. setosa.
summary(lm.iris1)##
## Call:
## lm(formula = Sepal.Width ~ Species, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.128 -0.228 0.026 0.226 0.972
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.42800 0.04804 71.359 < 2e-16 ***
## Speciesversicolor -0.65800 0.06794 -9.685 < 2e-16 ***
## Speciesvirginica -0.45400 0.06794 -6.683 4.54e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3397 on 147 degrees of freedom
## Multiple R-squared: 0.4008, Adjusted R-squared: 0.3926
## F-statistic: 49.16 on 2 and 147 DF, p-value: < 2.2e-16Post-hoc Multiple Comparisons
Okay, so we know that one group is different, but which one/s? It’s an excellent question; so let’s talk about post-hoc multiple comparison tests. There are lots out there, but today we will just use a Tukey test. Check out other multiple comparison options also, they may be better suited to your analyses.
In the first line of code below, we create another ANOVA object using aov(). This is very similar to using lm() and anova(), like we did above. We use aov() below because it is compatible with the function TukeyHSD() in the second line of code.
iris.aov<-aov(Sepal.Width~Species,data=iris)
TukeyHSD(iris.aov)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Sepal.Width ~ Species, data = iris)
##
## $Species
## diff lwr upr p adj
## versicolor-setosa -0.658 -0.81885528 -0.4971447 0.0000000
## virginica-setosa -0.454 -0.61485528 -0.2931447 0.0000000
## virginica-versicolor 0.204 0.04314472 0.3648553 0.0087802Look at the resulting Tukey table. You should see 4 columns and 3 rows. Each row is a comparison between 2 of the species. Look at the 4th column, p adj. This is the p-value for the comparison; less than 0.05 the two groups are different, greater than 0.05 and they are not significanlty different from one another.
I do not know the adjustment method used for these p-values, but this may be helpful from the R documentation for the function: “This function incorporates an adjustment for sample size that produces sensible intervals for mildly unbalanced designs.”
Two-way ANOVA
Now we will do a two-way ANOVA using the iris dataset again. This time, we will predict if sepal width differs among species and community.
The next code makes the community variable in the iris dataset.
v<-c("high","low")
iris$community<-vChecking Assumptions
Can check normality using the same code above.
Check variances
We will again use a Levene test.
library(car) #run this every time you open r and need the package
lm.iris2<-lm(Sepal.Width~Species*community,data=iris) #makes a 'linear model' object which we will use for the next line of code and the ANOVA code
leveneTest(lm.iris2)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 0.3307 0.8938
## 144The Levene p-value is greater than 0.05, so we are golden like Ponyboy (if you don’t get this reference, read The Outsiders by S.E. Hinton, it’s a fricking classic).
Running the two-way ANOVA
anova(lm.iris2)## Analysis of Variance Table
##
## Response: Sepal.Width
## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 11.3449 5.6725 48.7581 <2e-16 ***
## community 1 0.0067 0.0067 0.0573 0.8111
## Species:community 2 0.2025 0.1013 0.8704 0.4210
## Residuals 144 16.7528 0.1163
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1First of all, let’s look back at where we made lm.iris2. It was Sepal.Width~Speciescommunity. That ’‘tells R that we are testing for an interaction between Species and community. If you didn’t want to test the interaction, use a’+’ instead of the asterisk. Note: you almost always want to test the interaction.*
So now that we have that established, let’s look at the anova table. Let’s start in the Species:community row. The farthest right column (Pr(>F)) shows a p-value of 0.42, the interaction of the 2 variables is not significant. Now in the community row, we see a p-value of 0.81 - community does not have a significant effect on sepal width. Finally, we see that Species has a super small p-value (less than 2 X 10^-16). Species of iris has a significant effect on sepal width.
Post-hoc comparisons
Because the only signicant factor was Species, this post-hoc comparison will be the same as the one we did earlier.
iris2.aov<-aov(Sepal.Width~Species,data=iris)
TukeyHSD(iris2.aov)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Sepal.Width ~ Species, data = iris)
##
## $Species
## diff lwr upr p adj
## versicolor-setosa -0.658 -0.81885528 -0.4971447 0.0000000
## virginica-setosa -0.454 -0.61485528 -0.2931447 0.0000000
## virginica-versicolor 0.204 0.04314472 0.3648553 0.0087802But what if we did find a significant interaction of Species*community? The code below is what we would run. This would compare every combination of Species and community.
iris3.aov<-aov(Sepal.Width~Species*community,data=iris)
TukeyHSD(iris3.aov)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Sepal.Width ~ Species * community, data = iris)
##
## $Species
## diff lwr upr p adj
## versicolor-setosa -0.658 -0.81955127 -0.4964487 0.0000000
## virginica-setosa -0.454 -0.61555127 -0.2924487 0.0000000
## virginica-versicolor 0.204 0.04244873 0.3655513 0.0091385
##
## $community
## diff lwr upr p adj
## low-high -0.01333333 -0.1234264 0.09675978 0.8111496
##
## $`Species:community`
## diff lwr upr p adj
## versicolor:high-setosa:high -0.704 -0.98266057 -0.42533943 0.0000000
## virginica:high-setosa:high -0.544 -0.82266057 -0.26533943 0.0000013
## setosa:low-setosa:high -0.104 -0.38266057 0.17466057 0.8894028
## versicolor:low-setosa:high -0.716 -0.99466057 -0.43733943 0.0000000
## virginica:low-setosa:high -0.468 -0.74666057 -0.18933943 0.0000459
## virginica:high-versicolor:high 0.160 -0.11866057 0.43866057 0.5613861
## setosa:low-versicolor:high 0.600 0.32133943 0.87866057 0.0000001
## versicolor:low-versicolor:high -0.012 -0.29066057 0.26666057 0.9999958
## virginica:low-versicolor:high 0.236 -0.04266057 0.51466057 0.1475984
## setosa:low-virginica:high 0.440 0.16133943 0.71866057 0.0001553
## versicolor:low-virginica:high -0.172 -0.45066057 0.10666057 0.4800283
## virginica:low-virginica:high 0.076 -0.20266057 0.35466057 0.9692150
## versicolor:low-setosa:low -0.612 -0.89066057 -0.33333943 0.0000000
## virginica:low-setosa:low -0.364 -0.64266057 -0.08533943 0.0031579
## virginica:low-versicolor:low 0.248 -0.03066057 0.52666057 0.1113213t-test
A t-test looks for a difference in means between two groups.
First we will create a new dataset that only has two species of iris using the subset() function. This is only to create this example, not something you would need to do in the future.
iris.2<-subset(iris,Species!="virginica") #drops all virginica recordsTesting assumptions
A t-test has pretty much the same assumptions as ANOVA, so I’m going to breeze through them.
Checking normality
attach(iris.2)## The following objects are masked from iris:
##
## Petal.Length, Petal.Width, Sepal.Length, Sepal.Width, Specieshist(Sepal.Width)
shapiro.test(Sepal.Width)##
## Shapiro-Wilk normality test
##
## data: Sepal.Width
## W = 0.98977, p-value = 0.6463The histogram has a bit of right skew, but looks passable. The Shapiro test p-value is above 0.05, so it’s normally distributed!
Testing variances
leveneTest(Sepal.Width,Species,data=iris.2)## Levene's Test for Homogeneity of Variance (center = median: iris.2)
## Df F value Pr(>F)
## group 1 0.591 0.4439
## 98A Levene test gives a p-value of 0.44, above 0.05, so the variances are not different.
Note: if your variances ARE different, you can use a Welch’s t-test, which is friendly to violations of the variance equality assumption. This is the default in R with the function t.test, which I use below. In short, you can violate variance assumptions and use t.test() below.
Second note: we used a little different format for the Levene test this time. They are interchangeable: you could have made the lm (linear model) and ran the Levene’s test on it, like we did with the ANOVAs.
Running the t-test
t.test(Sepal.Width~Species,data=iris.2)##
## Welch Two Sample t-test
##
## data: Sepal.Width by Species
## t = 9.455, df = 94.698, p-value = 2.484e-15
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.5198348 0.7961652
## sample estimates:
## mean in group setosa mean in group versicolor
## 3.428 2.770The resulting t-test table is more easily interpreted than most tables in R. It gives you a t-value, degrees of freedom, p-value, and means of both groups. Because the p-value is less than 0.05, sepal width is significantly different between the two species.
Plotting t-test and ANOVA results
Please see the R plotting tutorial. This includes a tutorial for visualizing results of one-way ANOVA, two-way isn’t really covered unfortunately, but the package ggplot and google might help you out.
Now With Your Data: ANOVA (one and two-way) and t-test
When reading in your own data, set the working directory using setwd() or by the directions in Getting Started with R. Use read.csv() to read in your data.
And here is all of the code you will need, replace the variables with your own!
#ANOVA
#Testing Normality
attach(iris) #Attachs data to make it easier to call the data later
hist(Sepal.Width)
shapiro.test(Sepal.Width)
#If your data is not normal, check this out for transformations: http://fmwww.bc.edu/repec/bocode/t/transint.html
#If none of these work, check out the tutorial on nonparametric (Kruskal-Wallis and Wilcoxon) tests.
#Testing Variances among Groups
#If your data is normal, you can test variances.
library(car) #run this every time you open r and need the package
lm.iris1<-lm(Sepal.Width~Species,data=iris) #makes a 'linear model' object which we will use for the next line of code and the ANOVA code
#Remember, to make a 2-way ANOVA, you just change to lm.iris1<-lm(Sepal.Width~Species*Community,data=iris)
leveneTest(lm.iris1)
#The ANOVA
anova(lm.iris1)
summary(lm.iris1)
#Multiple comparisons: Tukey Test
iris.mc<-aov(Sepal.Width~Species,data=iris)
TukeyHSD(iris.mc)
#t-test
#Testing Normality
attach(iris.2)
hist(Sepal.Width)
shapiro.test(Sepal.Width)
#Test variances
leveneTest(Sepal.Width,Species,data=iris.2)
#The t-test
t.test(Sepal.Width~Species,data=iris.2)