The focus of this tutorial is analysis of variance (ANOVA). ANOVA is a statistical approach to compare means of an outcome variable of interest across different groups. If there are only two groups, the testing problem could simply use a \(t\)-test. But ANOVA makes it possible to test a multiple null hypothesis that compares several groups simultaneously, which would otherwise be tested by several \(t\)-tests with larger type I error.
The assumptions for parametric (fixed-effects) ANOVA are:
A first diagnostic test is to plot the residuals against the predicted responses and against treatment levels. Normality can be checked with a qq-plot as usual. Fortunately, ANOVA works well when the assumptions are nearly satisfied.
The parametric one-way ANOVA has a simple form for \(y_{ij}\), the \(j\)th observed response of group \(i\): \[ y_{ij} = u + \mu_i+\epsilon_{ij} \] where \(u\) is the grand mean, each group \(i\) has unique group mean \(u+\mu_i\), and \(\epsilon_{ij}\) is i.i.d. \(N(0,\sigma^2)\). If there is more than one grouping variable, say 2 grouping variables as an example, with the same assumptions, the two-way ANOVA has a slightly more complex form for \(y_{ijk}\), the \(k\)th observed reponse of group \(ij\): \[ y_{ijk} = u + \mu_i + \eta_j + \gamma_{ij} + \epsilon_{ijk} \] where \(u\) is again the grand mean, the mean of group \(ij\) is \(\mu_i + \eta_j + \gamma_{ij}\), and \(\epsilon_{ijk}\) is i.i.d. \(N(0,\sigma^2)\). The terms in the mean of group \(ij\) can be interpreted as follows: \(\mu_i\) is the mean effect of the first grouping variable, \(\eta_j\) is the mean effect of the second grouping variable, and the interaction term \(\gamma_{ij}\) is “the failure of the response to one factor to be the same at different levels of another factor” (see PSU site). ANOVA analysis is simply attributing the grand variations in the outcome variables into the variation caused by the specific group means and the unexplainable random noises. And then ANOVA compares the variations attributable to specific group means and those attributed to the noises and see if the group means are responsible for a relatively large size of the variation.
To do inference for one-way ANOVA, one uses the \(F\)-statistic for one-way ANOVA: \[F=\frac{\text{variation between sample means}}{\text{variation within the samples}}.\] For more on the \(F\)-statistic, hypothesis testing, and \(p\)-values, see this Minitab Blog entry.
Sometimes the assumptions for the parametric ANOVA above are not satisfied, and we could instead turn to a nonparametric counterpart of ANOVA, called Kruskal-Wallis test. The Kruskal-Wallis test simply transforms the original outcome variable data into the ranks of the data and then tests whether group mean ranks are different. Based on normality, the parametric ANOVA uses F-test while the Kruskal-Wallis test uses permutation test instead, which typically has more power in non-normal cases. In this tutorial, we would briefly go over one-way ANOVA, two-way ANOVA, and the Kruskal-Wallis test in R, STATA, and MATLAB.
Since the ANOVA could only tell us whether the group means of all groups are different, we still need to identify which groups are actually different by doing multiple comparisons across different group pairs. For ANOVA results, a specific multiple comparison approach could be used, which is called the Tukey honest significance test. The test simply compares pairwise group means but just like many other multiple comparison (\(p\)-value adjustment) methods, it controls for family-wise error rate (FWER). We will also demonstrate how to do this post-estimation step in all three languages. For a deeper and detailed introduction to one-way ANOVA and Tukey’s HSD, please refer to this. For two-way ANOVA, the tutorial at Penn State here is helpful. And for Kruskal-Wallis test, refer to Section 6.1 in book Nonparametric Statistical Methods by Hollander and Wolfe for a more theoretical treatment.
The data set we use for the illustration is the diet data, which is a small survey of several male and female adults who take on three different diet plans. And the focus of the survey is people’s body weight changes during a six-week time period. A link to the data is here. The self-descriptive variables of interest are the response variable weightloss, and the two grouping factors diet and gender. For ANOVA, we are considering the mean of weightloss because of gender, diet, and lastly because of the effect of diet according to gender (or vice-a-versa).
Variables of Interest | Description | Data type |
---|---|---|
weightloss | weight loss in kg | numeric |
diet | diet type | categorical |
gender | gender | categorical |
In the R, STATA, and MATLAB tutorials below, we first do a one-way ANOVA for weight loss for only females and for only males across the three diets and make a graph that shows the confidence intervals (adjusted for multicomparisons) for the three weight loss group means. Next we do two-way ANOVA for weight loss for the whole dataset across the three diets and gender. Finally, we again do the initial task of one-way ANOVA for weight loss for only females and for only males across the three diets, but use the Kruskal-Wallis test instead of parametric ANOVA.
The most important commands are aov() in R, anova in Stata, and multcompare in MATLAB.
First we want to load the dataset into R and omit any missing values using read.csv()
and na.omit()
.
diet.R = read.csv("./Diet.csv", header = TRUE)
diet.R = na.omit(diet.R)
The na.omit()
omitted two data points, Person 25 and 26. These points were missing the value for their gender.
Then we need to tell R that Diet and gender are factors by using as.factor()
and factor()
.
diet.R$Diet = as.factor(diet.R$Diet)
diet.R$gender = factor(diet.R$gender,c(0,1))
Then we need to create the variable of interest ‘weightlost’.
diet.R$weightlost = diet.R$pre.weight - diet.R$weight6weeks
summary(diet.R)
## Person gender Age Height pre.weight
## Min. : 1.00 0:43 Min. :16.00 Min. :141.0 Min. :58.00
## 1st Qu.:19.75 1:33 1st Qu.:32.50 1st Qu.:163.8 1st Qu.:66.00
## Median :40.50 Median :39.00 Median :169.0 Median :72.00
## Mean :39.87 Mean :39.22 Mean :170.8 Mean :72.29
## 3rd Qu.:59.25 3rd Qu.:47.25 3rd Qu.:175.2 3rd Qu.:78.00
## Max. :78.00 Max. :60.00 Max. :201.0 Max. :88.00
## Diet weight6weeks weightlost
## 1:24 Min. :53.00 Min. :-2.100
## 2:25 1st Qu.:61.95 1st Qu.: 2.300
## 3:27 Median :68.95 Median : 3.700
## Mean :68.34 Mean : 3.946
## 3rd Qu.:73.67 3rd Qu.: 5.650
## Max. :84.50 Max. : 9.200
Here is a quick summary of our dataset so you can see what variables are in the dataset and have a brief summary of those variables.
We then split the dataset into subsets by ‘gender’ using subset()
so we can have separate datasets for the female and male data.
Diet.female = subset(diet.R, gender==0)
Diet.male = subset(diet.R, gender==1)
Now that we have all the setup work done we can move on to our ANOVA analysis.
To create our ANOVA model we use ‘Diet’ as the grouping variable and use aov()
to create our parametric one-way ANOVA model for the female dataset. We can then use summary()
to analyze this model.
anovaFemale = aov(weightlost~Diet, data = Diet.female)
summary(anovaFemale)
## Df Sum Sq Mean Sq F value Pr(>F)
## Diet 2 92.32 46.16 10.64 0.000197 ***
## Residuals 40 173.53 4.34
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
From the above output of the parametric one-way ANOVA for the female dataset, we can see that the p-value is much smaller than 0.05 so we conclude that there is at least one group that is statistically different from the other groups in the female dataset.
In order to check which group is statistically different, we run the Tukey post hoc test for pairwise comparison following a one-way ANOVA using TukeyHSD()
.
TukeyHSD(anovaFemale)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = weightlost ~ Diet, data = Diet.female)
##
## $Diet
## diff lwr upr p adj
## 2-1 -0.4428571 -2.3589312 1.473217 0.8406368
## 3-1 2.8300000 0.9461311 4.713869 0.0020846
## 3-2 3.2728571 1.3889883 5.156726 0.0003833
From the above result we can see that there is a statistically significant difference in weight loss between the ‘Diet 3’ group and the ‘Diet 1’ group and between the ‘Diet 3’ group and the ‘Diet 2’ group. Because of this we can determine that the ‘Diet 3’ group is statistically different from the other two groups.
Now let’s do the same thing for the male dataset.
anovaMale = aov(weightlost~Diet, data = Diet.male)
summary(anovaMale)
## Df Sum Sq Mean Sq F value Pr(>F)
## Diet 2 2.0 1.001 0.148 0.863
## Residuals 30 202.8 6.760
From the above output of the parametric one-way ANOVA for the male dataset, we can see that the p-value is larger than 0.05 so we fail to reject the null hypothesis.
Let’s confirm this by looking at the Tukey post hoc test.
TukeyHSD(anovaMale)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = weightlost ~ Diet, data = Diet.male)
##
## $Diet
## diff lwr upr p adj
## 2-1 0.4590909 -2.341516 3.259698 0.9141710
## 3-1 0.5833333 -2.161144 3.327810 0.8602563
## 3-2 0.1242424 -2.551325 2.799809 0.9928028
As you can see there is no statistically significant difference between the three groups in the male dataset.
The following plots help visualize these two Tukey post hoc tests. Here is the plot for the female data.
library(gplots)
##
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
##
## lowess
plotmeans(weightlost~Diet, data = Diet.female)
And here is the plot for the male data.
plotmeans(weightlost~Diet, data = Diet.male)
For both of these plots n
is the number of people that participated in each Diet.
A similar way we could analysis this data is to do parametric two-way ANOVA for the whole dataset. For this we still use aov()
but now we regard ‘gender’ and ’Diet as two grouping variables. And to analysis the model we will use Anova()
and we want it to be a Type III Sums of Squares ANOVA.
library(car)
## Loading required package: carData
anova2 = aov(weightlost~gender*Diet,data=diet.R)
Anova(anova2,type = "III")
## Anova Table (Type III tests)
##
## Response: weightlost
## Sum Sq Df F value Pr(>F)
## (Intercept) 130.23 1 24.2247 5.489e-06 ***
## gender 2.10 1 0.3906 0.5340081
## Diet 92.32 2 8.5861 0.0004626 ***
## gender:Diet 33.90 2 3.1532 0.0488423 *
## Residuals 376.33 70
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The gender:Diet interaction is statistically significant at the p = 0.04884 level. There was no statistically significant difference in weight loss between gender (p = 0.77854), but there were statistically significant differences between Diet groups (p = 0.01304).
Now let’s look at Nonparametric ANOVA using the Kruskal-Wallis Test. We regard ‘Diet’ as the grouping variable and use kruskal.test()
to do nonparametric one-way ANOVA, i.e. Kruskal-Wallis test for the female data.
kruskal.test(weightlost~Diet, data = Diet.female)
##
## Kruskal-Wallis rank sum test
##
## data: weightlost by Diet
## Kruskal-Wallis chi-squared = 14.545, df = 2, p-value = 0.0006945
We get a p-value that is much smaller than 0.05 so we can reject the null hypothesis and conclude that there is at least one group statistically different from the other groups in the female dataset. This is the same conclusion we got in the parametric one-way ANOVA for the female data.
Now we do the same thing for the male dataset.
kruskal.test(weightlost~Diet, data = Diet.male)
##
## Kruskal-Wallis rank sum test
##
## data: weightlost by Diet
## Kruskal-Wallis chi-squared = 0.62218, df = 2, p-value = 0.7326
We get a p-value that is larger than 0.05 so there is no statistically significant difference in weight loss between the three groups in the male dataset. This is the same result we got from the parametric one-way ANOVA for the male data.
Finally, after all of this analysis we can conclude that for this dataset that you can either analyze the data using either parametric or nonparametric ANOVA and receive the same results. Our conclusion is that one group, the Diet 3 group, had a significant difference in weight loss compared to the other groups in the female dataset. However, for the male dataset there was no significant difference in weight loss between the different Diet groups.
Firstly, we load the dataset into Stata by using import delimited
.
import delimited Diet.csv
To observe the data preliminarily, we check the summary of data by using summarize
.
summarize
All variables seem reasonable except ‘gender’. Its abnormal summary probably is due to missing data. So we take a deeper look at the variable ‘gender’ by viewing the first ten rows of the data using list
.
list gender in 1/10
The variable ‘gender’ does have missing values. We delete those observations with missing values in ‘gender’ by drop if
and change ‘gender’ from string to numeric by using destring
(pass the replace
option in order to overwrite an existing file).
drop if gender == " "
destring gender, replace
Now the observations with missing values are removed. Then we generate the variable of interest ‘weightloss’ by generate
. We only keep those variables that we need by keep
and save the processed dataset by save
as we will use it later. Note that in Stata, only a single dataset can be opened at a time. As we want to operate on multiple datasets, we need to switch between them from time to time. So it is important that we save the data with save
and load the data with use
(pass the clear
option to remove the existing data as Stata will refuse if you try and open a new dataset with unsaved changes in the existing data).
* Generate the variable 'weightloss'
generate weightloss = preweight - weight6weeks
* Keep only the needed variables and save the dataset
keep gender diet weightloss
save diet.dta
We divide the dataset by ‘gender’ using keep if
and save the male data and female data as two separate datasets.
keep if gender == 0
save diet_female.dta
use diet.dta, clear
keep if gender == 1
save diet_male.dta
We regard ‘diet’ as the grouping variable and use the anova
command to do parametric one-way ANOVA for the female data.
use diet_female.dta, clear
anova weightloss diet
To help visualize the test for the female data, use the anovaplot
command following the one-way ANOVA. Note that to use the anovaplot
command, we should run ssc install modeldiag
first to load the corresponding package.
* ssc install modeldiag
anovaplot
From the above outputs of parametric one-way ANOVA for the female data, we can see that the p-value is much smaller than 0.05 . So we should reject the null hypothesis and conclude that there is at least one group statistically different from other groups in females.
In order to check that which group is statistically different, we run the Tukey post hoc test for pairwise comparison following a one-way ANOVA using pwmean DependentVar, over(IndependentVar), mcompare(tukey) effects
.
pwmean weightloss, over(diet) mcompare(tukey) effects
From the result above, we can see that there is a statistically significant difference in weight loss between the ‘diet 3’ group & ‘diet 1’ group and ‘diet 3’ group & ‘diet 2’ group in females. The ‘diet 3’ group is statistically different from other two groups.
Do the same things for the male data.
* Parametric one-way ANOVA for male data
* Regard 'diet' as the grouping variable
use diet_male.dta, clear
anova weightloss diet
* To visualize the test
anovaplot
From the above output of the parametric one-way ANOVA for the male dataset, we can see that the p-value is much larger than 0.05 , so we fail to reject the null hypothesis.
Let’s confirm this by running the Tukey post hoc test for pairwise comparison following a one-way ANOVA.
* Run the Tukey post hoc test for pairwise comparison
pwmean weightloss, over(diet) mcompare(tukey) effects
The results indicate that there is no statistically significant difference in weight loss among the three diet groups in males.
We could also use the anova
command to do parametric two-way ANOVA for the whole dataset. Now we regard ‘gender’ and ‘diet’ as two grouping variables.
use diet.dta, clear
anova weightloss i.gender##i.diet
To help visualize the test for the whole dataset, use the anovaplot
command again. Here type out the variable names following anovaplot
to confirm the order of predictors for better effect of the plot.
anovaplot diet gender
The gender#diet interaction is statistically significant at the p = 0.0488 level. There was no statistically significant difference in weight loss between gender (p = 0.7785), but there were statistically significant differences between diet groups (p = 0.0130).
We regard ‘diet’ as the grouping variable and use the kwallis
command to do nonparametric one-way ANOVA, i.e. Kruskal-Wallis test for the female data.
use diet_female.dta, clear
kwallis weightloss, by(diet)
We get a p-value much smaller than 0.05 . So we reject the null hypothesis and get the conclusion that there is at least one group statistically different from other groups in females, which is the same as what we get in parametric one-way anova for the female data.
Now we do the same thing for the male data.
use diet_male.dta, clear
kwallis weightloss, by(diet)
We get a p-value much larger than 0.05 . The results indicate that there is no statistically significant difference in weight loss among the three diet groups in males, which is also the same as what we get in parametric one-way anova for the male data.
In summary, we conclude that for one-way parametric and nonparametric ANOVA, there is some difference in weight loss of females and no statistically significant difference in weight loss of males among the three diet groups for the “diet” data. And for two-way parametric ANOVA, the interaction between ‘gender’ and ‘diet’ is statistically significant. There was no statistically significant difference in weight loss between gender, but there were statistically significant differences between diet groups.
Firstly, we load the data into the working space and convert the data into the type dataset
, remove the observations with missing value, create the dependent variable of interest and take a brief look at the data:
% import dataset
filename = 'Diet.csv';
Diet = readtable(filename);
Diet = table2dataset(Diet);
% remove missing values
diet_miss = ismissing(Diet);
Diet = Diet(~any(diet_miss,2),:);
% Dependent variable of interest:
Diet.weightloss = Diet.pre_weight - Diet.weight6weeks;
% first 5 rows of the data:
Diet(1:5,:)
Please note that in MALTAB, the basic data form is array, which cannot be operated upon as the way dataset is operated, so we need to transform the data table into dataset object. And using ismissing
to generate a 0-1 matrix of missing value indicator for each data cell in the dataset. And then using negation ~
to drop those missing observations.
In MATLAB, after we import a dataset, we simply call for the column of it by .
operator. And to subset the data, we could either use the column index as we do above, or we could use Diet(,'gender')
to refer to columns.
Because we want to conduct ANOVA analysis and firstly, we want to analyze the data under the categorical variable diet, and in the second step we are going to analyze the data under grouping variables gender and diet, so we need to subset the full sample into female and male subsamples for the first-step analysis:
% Subset of data under each gender
Diet_female = Diet(Diet.gender==0,:);
Diet_male = Diet(Diet.gender==1,:);
Firstly, for both female and male subsamples, we carry out one-way ANOVA analysis, which aims at identifying whether different diet plans lead to dispersed weightloss performances, and in MATLAB, we simply use anova1
command to indicate one-way anova:
[p_fe_aov,table_fe_aov,stats_fe_aov] = anova1(Diet_female.weightloss,Diet_female.Diet);
In anova1
command, just as many other MATLAB commands, the output(things on the left hand side) could be multidimensional. and each of the output could either be a value, a vector, a structure object, a table and many others. Here the p_fe_aov is the p-value of the ANOVA, table_fe_aov is the ANOVA result table, and the stats_fe_aov is a structure object storing relevant statistics for the dataset being analyzed, which is used for post-estimation command.
The result of the previous one-way anova is:
To see which group is actually different and lead to such a small p-value, we use the Tukey method to do pairwise-comparison and p-value adjustment by using the multcompare
command upon the stats_fe_aov object, which is used specifically for the postestimation command here:
[c_fe_aov,m_fe_aov,h_fe_aov,gnames] = multcompare(stats_fe_aov)
Here the output on the left is just some miscellaneous test statistic tables and p-values vectors which we would not go into details here. A confidence interval graph will automatically pop out if we run the command:
The y-axis is the indicator variable diet
, and the x-axis represents the group-average weightloss. Here we could see that the females taking on diet plan 3 is statistically different from the other two groups. To see it more clearly, we could take a look at the c_fe_aov object, which contains the pair-wise comparison statistics:
Here we could see that the pairwise comparison between group 1 and 3 and the comparison between group 2 and 3 are both significant, so group 3 is statistically different from the other groups.
We did exactly the same on the male subsample and get:
[p_ma_aov,table_ma_aov,stats_ma_aov] = anova1(Diet_male.weightloss,Diet_male.Diet);
[c_ma_aov,m_ma_aov,h_ma_aov,gnames] = multcompare(stats_ma_aov)
To analyze the full sample, we could also do a two-way anova, which simply takes on two grouping variables, and in MATLAB, the anovan
command could be used for general multi-anova(MANOVA) problem:
[p_full,table_full,stats_full] = anovan(Diet.weightloss,{Diet.gender Diet.Diet},'model','interaction','varnames',{'gender' 'diet type'});
[c_full,m_full,h_full,gnames] = multcompare(stats_full,'Dimension',[1 2])
To briefly the two commands, the anovan
function takes on the dependent variable, the grouping variables vector, and with a few options. A lot of MATLAB functions include options via so-called key-value pairs. Here the two options are:
*(‘model’,‘interaction’): set the ANOVA model to include interaction effect of the two grouping variables
*(‘varnames’,{‘gender’,‘diet type’}): reset the label of the variables in the ANOVA table
And in the multcompare
command, there is another option which just indicates that the grouping variables are 2-dimensional. We could also add a key-value pair (‘ctype’,‘scheffe’/‘bonferroni’) to change the p-value adjustment method to scheffe method and bonferroni method. The default is Tukey honestly significant difference criterion, which is exactly what we want.
The results are intuitive:
It shows that only the female group taking on diet plan 3 is causing the overall difference, which is sort of equivalent to the two one-way anova’s results combined.
To get a more robust hypothesis test of the female and male subsamples, we could simply do the one-way anova test in a non-parametric way using the kruskal-wallis test where we simply implement it by kruskalwallis()
function:
[p_fe_kw,table_fe_kw,stats_fe_kw] = kruskalwallis(Diet_female.weightloss,Diet_female.Diet)
[p_ma_kw,table_ma_kw,stats_ma_kw] = kruskalwallis(Diet_male.weightloss,Diet_male.Diet)
The kruskalwallis()
function works in the same manner as anova1
which only asks the outcome variable and the grouping variable, the result is resonating that in parametric ANOVA. For female subsample:
And for male subsample:
And the multicomparison works in the same way:
[c_full,m_full,h_full,gnames] = multcompare(stats_fe_kw)
[c_full,m_full,h_full,gnames] = multcompare(stats_ma_kw)
Result for female subsample:
Result for male subsample:
The conclusion from non-parametric test is exactly the same as that in parametric test.
In general, the analysis in MATLAB on the diet dataset yields exactly the same result as the other two examples in STATA and R. We have found significant differences in females’ different diet groups but do not observe a similar pattern for males. And furthermore, when we consider the full dataset, it seems that the group that stands out is still the female group that makes the difference in the female subsample. The result holds for both parametric ANOVA and Kruskal-Wallis test.
The parametric ANOVA is powerful and has several post estimation commands for further analysis. But before doing the parametric ANOVA, especially when working with a dataset with small sample size, we need to carefully check the assumptions of ANOVA. Specifically, we should check whether data in each group follows normal distribution using Kolmogorov-Smirnov test and also test whether the variances of data are equal across groups by doing Levene test.
The other thing to consider is how to do pairwise comparisons for Kruskal-Wallis test. We could do pairwise Wilcoxon rank sum test for each group pair, but it seems unclear to us how to adjust the p-value to control for the overall FWER. The Tukey HSD applied to the parametric ANOVA object seems not applicable to Kruskal-Wallis since it requires the same parametric assumptions as those in parametric ANOVA.
Pennsylvania State University Online Course Statistics 500: Applied Statistics [see Section 11.2]
Myles Hollander and Douglas A. Wolfe. Nonparametric Statistical Methods, Wiley Series in Probability and Statistics, 1999. [see Section 6.1]
Paul Teetor. R Cookbook, O’Reilly Media, 2011. [See Sections 11.20, 11.23]