Exercise 7.1: ANOVA on the lettuce dataset - solution

1 Background

In agriculture, it is important to have a high yield of crops. For lettuce plants, plants with many leaves are known to be preferred by the consumers.

One way to increase the number of leaves (or better, total leaf weight) is by using a fertilizer. Recently, there has been 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.

2 The lettuce dataset

The researchers hypothesize that biochar, compost and a combination of both biochar and compost can have a different 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 assess if the treatments have a different effect on the growth of lettuce plants. If so, they will use a post-hoc test (Tuckey test) to discover which specific treatments have an effect.

Load the required libraries

library(tidyverse)
theme_set(theme_bw()) # Change ggplot2 default theme

library(multcomp)

3 Data import

lettuce <- read_csv("https://raw.githubusercontent.com/statOmics/PSLSData/main/freshweight_lettuce.txt")
# Take a glimpse at the data
glimpse(lettuce)

## Rows: 28
## Columns: 3
## $ id          <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 1…
## $ treatment   <chr> "control", "control", "control", "control", "con…
## $ freshweight <dbl> 38, 34, 41, 43, 43, 29, 38, 59, 64, 57, 56, 50, …

First, we will set the treatment variable to a factor

lettuce <- lettuce %>%
  mutate(
    ## Convert treatment to factor
    treatment = as.factor(treatment),
    ## Relevels such that "control" is the reference
    treatment = relevel(treatment, ref = "control")
  )

4 Data exploration

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

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

boxplot <- lettuce %>%
  ggplot(aes(x = treatment, y = freshweight, fill = treatment)) +
  scale_fill_brewer(palette = "RdGy") +
  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 = mean, geom = "point",
    shape = 5, size = 3, color = "black"
  )
boxplot

Note that there are no clear outliers in the data. We can see that the mean freshweight 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?

5 ANOVA

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

5.1 Formulate null and alternative hypotheses

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.

5.2 Checking the assumptions for ANOVA

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

5.2.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.

5.2.2 Assumption of normality

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

ggplot(lettuce, aes(sample = freshweight)) +
  geom_qq() +
  geom_qq_line(aes(col = treatment), show.legend = FALSE) +
  scale_color_brewer(palette = "Dark2") +
  facet_wrap(vars(treatment))

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)

## sample 7 observations from a uniform distribution, repeat 6 times
df <- data.frame(x = runif(7 * 6), rep = rep(1:6, each = 7))
ggplot(df, aes(sample = x)) +
  geom_qq() +
  geom_qq_line(col = "darkorchid") +
  facet_wrap(vars(rep)) +
  ggtitle("QQ plots for uniformly distributed data")

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

Sample from normal distribution:

set.seed(5)
## sample 7 observations from a normal distribution, repeat 6 times
df <- data.frame(x = rnorm(7 * 6), rep = rep(1:6, each = 7))
ggplot(df, aes(sample = x)) +
  geom_qq() +
  geom_qq_line(col = "dodgerblue") +
  facet_wrap(vars(rep)) +
  ggtitle("QQ plots for normally distributed data")

Even though the data originates from a normal distribution, some deviations can be observed due to random chance.

Alternatively, 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 Q-Q plot of the residuals shows a slight skew to the left. Care should be taken when interpreting the results for this dataset. Due to the sample size, a better alternative might be a non-parametric approach.

5.2.3 Assumption of equal variances

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

boxplot

As a measure of variability, we may take the height of each 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 these assumptions.

5.3 Modeling

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

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.

5.4 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.

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|)    
## cobc - control == 0      16.143      2.766   5.837   <0.001 ***
## compost - control == 0   20.857      2.766   7.541   <0.001 ***
## refoak - control == 0    -1.857      2.766  -0.671    0.907    
## compost - cobc == 0       4.714      2.766   1.704    0.343    
## refoak - cobc == 0      -18.000      2.766  -6.508   <0.001 ***
## refoak - compost == 0   -22.714      2.766  -8.213   <0.001 ***
## ---
## 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     
## cobc - control == 0     16.1429   8.5101  23.7757
## compost - control == 0  20.8571  13.2243  28.4899
## refoak - control == 0   -1.8571  -9.4899   5.7757
## compost - cobc == 0      4.7143  -2.9185  12.3471
## refoak - cobc == 0     -18.0000 -25.6328 -10.3672
## refoak - compost == 0  -22.7143 -30.3471 -15.0815

5.5 Conclusion

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

The mean freshweight 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 freshweight 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.

---

This work is licensed under the CC BY-NC-SA 4.0 licence.