One-way ANOVA PlantGrowth
Note: This article shows an R example on how to conduct a one-way analysis of variance (ANOVA) on the dataset PlantGrowth. For general information about ANOVA please refer to the following article: ANOVA.
In short: In this article an exemplary one-way ANOVA is conducted on the dataset PlantGrowth which examines if the yields of the plants dependent on different treatments differ. The results are further investigated with post-hoc tests and visualised in boxplots. The code as well as the outputs are provided with explanations.
Contents
The Dataset[edit]
the dataset PlantGrowth is included in the R base package and contains data on the yield of crops under two different treatments and one control group.
The variables we will use in our analysis are:
- weight: a numeric variable containing the dried weight of the plants
- group: a factor containing the treatments “ctrl”,”trt1”, “trt2”
Examining the Dataset[edit]
Before conducting the analysis, we can examine the dataset and check for the different assumptions of an ANOVA to make sure that the ANOVA produces valid results. The assumptions of the ANOVA are:
- Normality of the dependent variable (weight)
- Equal variance for each factor
- Independent measurements
- Equal size in each treatment group
Before running the code, you might need to install the packages using the command install.packages() eg. install.packages("multcompView").
The command head() displays the first entries of the data frame, the command str() display the internal structure, while summary() displays information about the different variables.
To check for the equal size of groups, the summary() command can be used. The independent measurements are a property of the data, which is given in this case. The normality of the dependent variable can be analysed with a shapiro-wilk test while the equal variances can either be tested with a levene test or visualised with boxplots (Fig. 2).
#Install or load necessary libraries based on your requirements library(ggplot2) library(dplyr) library(multcompView) library(dplyr) #Inspect the dataset head(PlantGrowth) str(PlantGrowth) summary(PlantGrowth) #check for assumptions #normality of dependent variable shapiro.test(PlantGrowth$weight) #historgram (Fig. 1) hist(plantGrowth$weight) #equal variance between groups boxplot(PlantGrowth$weight~PlantGrowth$group)
Output[edit]
> head(PlantGrowth)
weight group
1 4.17 ctrl
2 5.58 ctrl
3 5.18 ctrl
4 6.11 ctrl
5 4.50 ctrl
6 4.61 ctrl
> str(PlantGrowth)
'data.frame': 30 obs. of 2 variables:
$ weight: num 4.17 5.58 5.18 6.11 4.5 4.61 5.17 4.53 5.33 5.14 ...
$ group : Factor w/ 3 levels "ctrl","trt1",..: 1 1 1 1 1 1 1 1 1 1 ...
> summary(PlantGrowth)
weight group
Min. :3.590 ctrl:10
1st Qu.:4.550 trt1:10
Median :5.155 trt2:10
Mean :5.073
3rd Qu.:5.530
Max. :6.310
> shapiro.test(PlantGrowth$weight)
Shapiro-Wilk normality test
data: PlantGrowth$weight
W = 0.98268, p-value = 0.8915
The summary shows that we have an equal size of measurements for each group (10 per group). The Shapiro-Wilk test delivers a non-significant p-value which means we cannot reject the 0-hypothesis that the data is normally distributed. We thus assume a normal distribution of the yield. The boxplots shows the variance within the yield for each factor. The variances are not quite the same which means other variables potentially influence the yield. The results of the following analysis therefore need to be interpreted with care.
ANOVA[edit]
The ANOVA compares the means of the weight for the different treatment groups. This way we can find out if the treatment has an effect on the weight and thus the growth of the plants.
A one-way ANOVA is performed to test if the means of the different treatment groups, the yields dependent on the treatments, differ.
#Perform one-way ANOVA anova_model <- aov(weight ~ group, data = PlantGrowth) #summary of ANOVA summary(anova_model)
Output[edit]
> summary(anova_model)
Df Sum Sq Mean Sq F value Pr(>F)
group 2 3.766 1.8832 4.846 0.0159 *
Residuals 27 10.492 0.3886
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The summary of the ANOVA model delivers the following results:
- Df: The Degrees of freedom are a measure of the Independent values that can vary in our design. For each factor they are calculated by n-1.
- Sum Sq: The sum of squares is a measure of how much variance is explained by the independent variable. The group Sum Sq plus the residuals Sum Sq equals to the total variation. If you divide the Sum Sq of one variable by the total variation, you get the explained variation by variable group. The treatment explains 3.766/14.258 of the variance within the yield.
- Pr: The P - value indicates the probability that the 0 hypothesis is true. Since the p-value is below 0.05, we can reject the 0-hypothesis that there is no significant difference between the different treatment groups and thus assume a significant effect of the treatment on the yield.
Post-hoc Test[edit]
An ANOVA simply analyses if there is a difference between all the means. A post-hoc test is testing for the differences of the means in the individual groups, correcting for errors in the significance levels when we perform multiple tests. If multiple tests are performed within the same treatment group (for example different t-tests testing for the difference between each of the treatment groups) the probability to falsely reject a 0-hypothesis and deem a non-significant result as significant is increased.
for the post-hoc test the tukeyHSD test is used. It only needs to be applied if the ANOVA itself delivers a significant result.
With the use of the library multicompView, letters are assigned to each treatment group which indicate if they differ significantly. If two groups share no letter, there is a significant difference between them.
# post-hoc test tukey_res <- TukeyHSD(anova_model) tukey_res #Create compact letter display for group differences tukey_letters <- multcompView::multcompLetters4(anova_model, tukey_res) tukey_letters
To illustrate the need for a tukey test, t-tests for each of the group differences are performed and compared with the results of the tukeyHSD test.
#t-tests as comparisons
t.test(weight ~ group,
data = subset(PlantGrowth, group %in% c("ctrl", "trt1")))
t.test(weight ~ group,
data = subset(PlantGrowth, group %in% c("ctrl", "trt2")))
t.test(weight ~ group,
data = subset(PlantGrowth, group %in% c("trt2", "trt1")))
Output[edit]
> tukey_res
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = weight ~ group, data = PlantGrowth)
$group
diff lwr upr p adj
trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
trt2-ctrl 0.494 -0.1972161 1.1852161 0.1979960
trt2-trt1 0.865 0.1737839 1.5562161 0.0120064
> tukey_letters
$group
trt2 ctrl trt1
"a" "ab" "b"
The adjusted p-values indicate that there only is a significant difference between the two treatment groups. This can also be seen in the letter display as trt1 and trt2 share no letter. For our dataset this means that the treatment applied to the plants in both cases has no effect in comparison to the control group. Only the two treatments are significantly different from each other which implies that one treatment has a negative effect on the yield while the effect of the other is positive. Because of the minor effects of the treatments, these would not be used to increase the crop yield.
The t-tests deliver the following results:
> t.test(weight ~ group,
+ data = subset(PlantGrowth, group %in% c("ctrl", "trt1")))
Welch Two Sample t-test
data: weight by group
t = 1.1913, df = 16.524, p-value = 0.2504
alternative hypothesis: true difference in means between group ctrl and group trt1 is not equal to 0
95 percent confidence interval:
-0.2875162 1.0295162
sample estimates:
mean in group ctrl mean in group trt1
5.032 4.661
> t.test(weight ~ group,
+ data = subset(PlantGrowth, group %in% c("ctrl", "trt2")))
Welch Two Sample t-test
data: weight by group
t = -2.134, df = 16.786, p-value = 0.0479
alternative hypothesis: true difference in means between group ctrl and group trt2 is not equal to 0
95 percent confidence interval:
-0.98287213 -0.00512787
sample estimates:
mean in group ctrl mean in group trt2
5.032 5.526
> t.test(weight ~ group,
+ data = subset(PlantGrowth, group %in% c("trt2", "trt1")))
Welch Two Sample t-test
data: weight by group
t = -3.0101, df = 14.104, p-value = 0.009298
alternative hypothesis: true difference in means between group trt1 and group trt2 is not equal to 0
95 percent confidence interval:
-1.4809144 -0.2490856
sample estimates:
mean in group trt1 mean in group trt2
4.661 5.526
Generally, the p-values of the t-tests are much lower than the ones that the tukey test delivers. In the case of “crtl” and “trt2” this leads to the difference being regarded as significant in the t-test while the tukey tests deems it as non-significant due to the correction for multiple testing.
Plotting the Results[edit]
The results of the ANOVA and the post hoc test are visualised in a boxlpot including error bars which display the letters that were created with the TukeyHSD test above the bars to visualise the ANOVA results and display the information which treatment groups differ from each other.
#boxplots with error bars and tukeyHSD letters (Fig. 3)
#Prepare letter positions
letters_tukey <- data.frame(
group = c("ctrl", "trt1", "trt2"),
letter = c("ab", "b", "a")
)
pos_tukey <- PlantGrowth %>%
group_by(group) %>%
summarise(
y_pos = max(weight) + 0.2 # shift above boxplots
)
letters_tukey <- left_join(letters_tukey, pos_tukey, by="group")
# Plot
ggplot(PlantGrowth, aes(group, weight)) +
stat_boxplot(geom = 'errorbar', width=0.5) +
geom_boxplot(outlier.shape = 1) +
stat_summary(fun.y = mean, geom="point", size=2) +
stat_summary(fun.data = mean_se, geom="errorbar") +
geom_text(data = letters_tukey,
aes(x = group, y = y_pos, label = letter),
size = 5) +
labs(
title = "Effect of Treatment on Plant Growth",
x = "Treatment Group",
y = "Mean Plant Weight (g)"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
axis.title.x = element_text(face = "bold"),
axis.title.y = element_text(face = "bold")
)
Diagnostic plots[edit]
To check for the validity of the results, the diagnostic plots can be examined. The 4 diagnostic plots are printed in a 2*2 grid
#Optional: Diagnostic plots par(mfrow = c(2, 2)) plot(anova_model) par(mfrow = c(1, 1))
The residuals vs. fitted should ideally be equally distributed for each factor. In this model there seems to be roughly the same amount of variance for each factor although the distribution could be more equal. In the Q-Q plot, the residuals should roughly follow the line. We can still see the factors here which is unavoidable in an ANOVA. The diagnostic plots used here are hence more appropriate for a regression analysis. Alternatively, one can look at the distribution of the residuals on all factor levels, ideally with boxplots or histograms.