Tutorial 7.1: ANOVA in the lettuce dataset

To experiment with the concepts of ANOVA, we will work with a dataset on the freshweight of lettuce plants.

1 The lettuce dataset

In agriculture, it is important to have a high yield of crops. For lettuce plants, this means having as many leaves as possible per plant. This is known to preferred by the consumers.

One way to increase the number of leaves (or better, total leaf weigth) is by using a fertilizer. Recently, there is a tendency to rely more on natural fertilizers, such as compost. Near Ghent, the institute of research for agriculture and fishery is testing new, natural fertilization methods. One of these new fertilizers is called biochar. Biochar is a residual product from pyrolysis, a process in which biomass is burned under specific conditions (such as high pressure) in order to produce energy. Biochar is similar to charcoal, but has some very useful properties, such as retaining water in the soil. It also has a positive influence on the soil microbiome.

The researcher want to find out if biochar, compost and a combination of both biochar and compost have an influence on the growth of lettuce plants. To this end, they grew up lettuce plants in a greenhouse. The pots were filled with one of four soil types;

1. Soil only (control)
2. Soil supplemented with biochar (refoak)
3. Soil supplemented with compost (compost)
4. Soil supplemented with both biochar and compost (cobc)

The dataset freshweight_lettuce.txt contains the freshweight (in grams) for 28 lettuce plants (7 per condition). The researchers want to use an ANOVA test to find out whether or not there is an effect of one or more of the treatments on the growth of lettuce plants. If so, they will use a post-hoc test (Tuckey test) to find which specific treatments have an effect.

Load the required libraries

library(tidyverse)
library(car)

2 Data import

lettuce <- read_csv("https://raw.githubusercontent.com/GTPB/PSLS20/master/data/freshweight_lettuce.txt")
#Take a glimpse at the data
glimpse(lettuce)

## Observations: 28
## Variables: 3
## $ id          <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…
## $ treatment   <chr> "control", "control", "control", "control", "control", "c…
## $ freshweight <dbl> 38, 34, 41, 43, 43, 29, 38, 59, 64, 57, 56, 50, 64, 62, 3…

First, we will set the treatment variable to a factor

## treatment to factor
lettuce <- lettuce %>%
  mutate(treatment = as.factor(treatment))

3 Data exploration

## Count the number of observations per treatment
lettuce %>%
  count(treatment)

## # A tibble: 4 x 2
##   treatment     n
##   <fct>     <int>
## 1 cobc          7
## 2 compost       7
## 3 control       7
## 4 refoak        7

Now let’s make a boxplot displaying the freshweight of each treatment condition:

lettuce %>%
  ggplot(aes(x=treatment,y=freshweight,fill=treatment)) + 
  scale_fill_brewer(palette="RdGy") +
  theme_bw() +
  geom_boxplot(outlier.shape=NA) + 
  geom_jitter(width = 0.2) +
  ggtitle("Boxplot of the freshweigth for each treatment condition") +
  ylab("freshweight (gram)") + 
  stat_summary(fun.y=mean, geom="point", shape=5, size=3, color="black", fill="black")

Note that there are no clear outliers in the data. We can see that the mean fresweight is very comparable between the control and refoak treatments and between the compost and cobc treatments. We can also see that the mean freshweight is much higher in the cobc and control treatments than in the control and refoak treatments. But is this observed difference significant?

3.1 Checking the assumptions for ANOVA

To study whether or not the observed difference between the average freshweight values of the different treatment groups are significant, we may perform an ANOVA.

The null hypothesis of ANOVA states that: \(H0\): The mean freshweigth is equal between the different treatment groups.

The alternative hypothesis of ANOVA states that: \(HA\): The mean freshweigth for at least one treatment group is different from the mean freshweight in at least one other treatment group.

Before we may proceed with the analysis, we must make sure that all assumptions for ANOVA are met. ANOVA has three assumptions:

1. The observations are independent of each other (in all groups)
2. The data (freshweigth) must be normally distributed (in all groups)
3. The variability within all groups is similar

3.1.1 Assumption of independence

The first assumption is met; we started of with 28 lettuce plants and we randomly submitted them to one of four treatment conditions. There is no reason to believe that the plants display systematic differences between treatment groups, other than the actual treatment.

Question: What would happen if all of the control and refoak plants were grown in 1 greenhouse and the other two conditions were grown in another greenhouse?

3.1.2 Assumption of normality

For the second assumption, we must check normality in each group.

## get qqplots for each individual treatment group
par(mfrow = c(2,2))
for(i in levels(lettuce$treatment)){
  qqPlot(subset(lettuce,treatment == i)$freshweight, main = i, ylab = "")
}

According to these plots, the normality assumption is met for each group. However, we need to be careful with relying on these plots, as they are only based on 7 observations per group. Indeed, checking the assumptions of normality (and equal variances) on such small datasets can be tricky. For instance, when sampling 7 datapoints from data that has an underlying distribution that is not normally distributed, the small sample may be normally distributed by chance. Equivalently, when sampling 7 datapoints from a normal distribution, the small sample may appear not normally distributed. We can show this as follows.

Simulate from uniform distribution:

set.seed(5)
par(mfrow = c(3,2),mar=c(2, 4, 2, 0.2))
for (i in 1:6) {
  x <- runif(7) ## sample 7 observations from a uniform distribution
  qqPlot(x,ylab = "")
}

Even though the data originates from a uniform distribution, several samples meet the normality requirements by chance.

Sample from normal distribution:

set.seed(5)
par(mfrow = c(3,2),mar=c(2, 4, 2, 0.2))
for (i in 1:6) {
  x <- rnorm(7) ## sample 7 sobservations from a normal distribution
  qqPlot(x,ylab = "")
}

Even though the data originates from a normal distribution, several samples do not meet the normality requirements by chance.

As an alternative, we may generate a qq-plot of the residuals of a linear model. These residuals should be normally distributed if the data for each treatment condition is normally distributed.

fit <- lm(freshweight~treatment, lettuce)
plot(fit, which = 2,col = fit$model$treatment)
legend('bottomright', levels(fit$model$treatment), text.col = 1:4)

The data seems to be approximately normally distributed, with a slight left skew.

3.1.3 Assumption of equal variances

We can check the assumption of equal variance with a boxplot:

lettuce %>%
  ggplot(aes(x=treatment,y=freshweight,fill=treatment)) + 
  scale_fill_brewer(palette="RdGy") +
  theme_bw() +
  geom_boxplot(outlier.shape=NA) + 
  geom_jitter(width = 0.2) +
  ggtitle("Boxplot of the freshweigth for each treatment condition") +
  ylab("freshweight (gram)") + 
  stat_summary(fun.y=mean, geom="point", shape=5, size=3, color="black", fill="black")

As a measure of variability, we may take the height of each boxplot’s box. This is the interval between the 25% and 75% quantile. Here we can see that this interval, as well as the length of the whiskers, is approximately equal for most groups. However, the variability of cobc does seem to be quite a bit larger than the variability in the refoak group.

With this little observations per group, it is very difficult to reliably assess the assumptions of normality and equal variances. In this tutorial, we will assume that all assumptions are met. In a later tutorial, “Kruskal_lettuce_plants.Rmd”, we will see would happen if we decide to reject the assumptions.

4 Data analysis

4.1 ANOVA

fit <- lm(freshweight~treatment, data=lettuce)
fit_anova <- anova(fit)
fit_anova

## Analysis of Variance Table
## 
## Response: freshweight
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## treatment  3 2732.14  910.71  34.015 8.308e-09 ***
## Residuals 24  642.57   26.77                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The p-value of the ANOVA analysis is extremely significant (p-value = 8.308e-09), so we reject the null hypothesis that the mean freshweigth is equal for the different treatment groups. We can say that the mean freshweigth is significantly different between at least two treatment groups on the 5% significance level.

Based on this analysis, we do not yet know between which particular groups there is a significant difference. To assess this, we will perfrom the Tuckey post-hoc analysis.

4.2 Post-hoc analysis

We will perform a post-hoc analysis to look at the difference in fresweigth between each pairwise comparison of treatment groups. Importantly, with this strategy, the p-values will be automatically adjusted for multiple testing.

The null hypothesis for each pairwise test is that \(H0\) there is no difference in the average freshweight values between both groups.

The alternative hypothesis for each pairwise test states that \(HA\) there is indeed a difference in the average freshweight values between both groups.

library(multcomp, quietly = TRUE)
mcp <- glht(fit, linfct = mcp(treatment = "Tukey"))
mcp_summary <- summary(mcp)
mcp_summary

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = freshweight ~ treatment, data = lettuce)
## 
## Linear Hypotheses:
##                        Estimate Std. Error t value Pr(>|t|)    
## compost - cobc == 0       4.714      2.766   1.704    0.343    
## control - cobc == 0     -16.143      2.766  -5.837   <0.001 ***
## refoak - cobc == 0      -18.000      2.766  -6.508   <0.001 ***
## control - compost == 0  -20.857      2.766  -7.541   <0.001 ***
## refoak - compost == 0   -22.714      2.766  -8.213   <0.001 ***
## refoak - control == 0    -1.857      2.766  -0.671    0.907    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

#We will also calculate the confidence interval on the mean differences.
mcp_confint <- confint(mcp)
mcp_confint

## 
##   Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = freshweight ~ treatment, data = lettuce)
## 
## Quantile = 2.7597
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##                        Estimate lwr      upr     
## compost - cobc == 0      4.7143  -2.9185  12.3471
## control - cobc == 0    -16.1429 -23.7757  -8.5101
## refoak - cobc == 0     -18.0000 -25.6328 -10.3672
## control - compost == 0 -20.8571 -28.4899 -13.2243
## refoak - compost == 0  -22.7143 -30.3471 -15.0815
## refoak - control == 0   -1.8571  -9.4899   5.7757

5 Conclusion

We have found an extremely significant dependence (p-value = 8.30810^{-9}), between the mean freshweigth and the treatment condition on the global 5% significance level.

The mean fresweight of plants grown with compost is extremely significantly higher as compared to in control plants (Tuckey test, adjusted p-value < 0.0001, 20.9g higher, 95% CI: [13.2; 28.5]) and as compared to refoak plants (Tuckey test, adjusted p-value < 0.0001, 22.7g higher, 95% CI: [15.1; 30.3]).

The mean fresweight of plants grown with cobc is extremely significantly higher as compared to in control plants (Tuckey test, adjusted p-value < 0.0001, 16.1g higher, 95% CI: [8.5; 23.8]) and as compared to refoak plants (Tuckey test, adjusted p-value < 0.0001, 18.0g higher, 95% CI: [10.4; 25.6]).

We do not find enough evidence to claim a difference in mean freshweight between refoak and control plants or between cobc and compost plants.

We may conclude that supplementing soil with compost or with both compost and biochar will have a positive effect on the freshweigth of lettuce plants.

Question: can we extrapolate these results to another vegetable for instance, spinache?