Exercise 7.2: ANOVA on the cuckoo dataset - solution

1 Cuckoo dataset

The common cuckoo does not build its own nest: it prefers to lay its eggs in another birds’ nest. It is known, since 1892, that the type of cuckoo bird eggs are different between different locations. In a study from 1940, it was shown that cuckoos return to the same nesting area each year, and that they always pick the same bird species to be a “foster parent” for their eggs.

Over the years, this has lead to the development of geographically determined subspecies of cuckoos. These subspecies have evolved in such a way that their eggs look as similar as possible as those of their foster parents.

The cuckoo dataset contains information on 120 Cuckoo eggs, obtained from randomly selected “foster” nests. For these eggs, researchers have measured the length (in mm) and established the type (species) of foster parent. The type column is coded as follows:

• type=1: Meadow pipit
• type=2: Tree pipit
• type=3: Dunnock
• type=4: European robin
• type=5: White wagtail
• type=6: Eurasian wren

2 Goal

The researchers want totest if the type of foster parent has an effect on the average length of the cuckoo eggs.

In theory, they want to study this for all six species. Previously, we looked at a single pairwise comparison between the European robin and the Eurasian wren with a t-test. Here, we will analyse all types simultaneously with ANOVA.

Load the required libraries

library(tidyverse)

2.1 Import the data

Cuckoo <- read_tsv("https://raw.githubusercontent.com/statOmics/PSLS21/data/Cuckoo.txt")
head(Cuckoo)

3 Data tidying

It seems that the tpye column is a double rather than a factor. Let’s fix this:

Cuckoo <- Cuckoo %>%
  mutate(type = as.factor(type)) %>%
  mutate(type = recode(type, 
           "1" = "meadow", 
           "2" = "tree", 
           "3" = "dunnock",  
           "4" = "robin", 
           "5" = "wagtail", 
           "6" = "wren"))

4 Data exploration

How many birds do we have for each type?

Cuckoo %>%
  count(type)

Visualize the data

Cuckoo %>%
  ggplot(aes(x=type,y=length,fill=type)) +
    theme_bw() +
    scale_fill_brewer(palette="RdGy") +
    geom_boxplot() +
    geom_jitter(width = 0.2) +
    ggtitle("Boxplot of the length of eggs per type") +
    ylab("length (mm)") + 
    stat_summary(fun=mean, geom="point", shape=5, size=3, color="black", fill="black")

5 ANOVA

To study if the observed differences in average egg length between the different foster bird types are significant, we may perform an ANOVA.

5.1 Formulate null and alternative hypothesis

The null hypothesis of ANOVA states that: \(H0\): The mean egg length is equal between the different bird types.

The alternative hypothesis of ANOVA states that: \(HA\): The mean egg length for at least one bird type is different from the mean egg length in at least one other bird type.

5.2 Parameter estimation

mod1 <- lm(length~type, Cuckoo)

5.3 Check assumptions

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 (length) must be normally distributed (in all groups)
3. The variability within all groups is similar

The first assumption is met, as we may assume that there are no specific patterns of correlation between the randomly selected nests.

To check the normality assumption, we will use QQ plots.

plot(mod1, which = 2)

Cuckoo %>% 
  ggplot(aes(sample=length)) +
  geom_qq() +
  geom_qq_line() + 
  facet_grid(~type)

There seem to be no clear deviations from normality.

The third assumption of equal variances seems to be met based on the visualization with the boxplots (see above).

As such, we may proceed with the ANOVA analysis.

5.3.1 Simulate to train your skills in assessing assumptions

set.seed(1031)
sigma <- mod1 %>% sigma 
dataHlp <- Cuckoo

simModels <- list()
plotList <- list()
plotListQQ1 <- list() 

par(mfrow=c(3,3))
for (i in 1:9)
{
  nobs <- Cuckoo %>% nrow
  dataHlp$ySim <- mod1$fit + rnorm(nobs, sd = sigma)
  simModels[[i]] <- lm(ySim~type, dataHlp)
  plot(simModels[[i]], which = 2)
  plotList[[i]] <- dataHlp %>%  
    ggplot(aes(type,ySim)) + 
    geom_boxplot(outlier.shape = NA) +
    geom_jitter()
  plotListQQ1[[i]] <- dataHlp %>%  
    filter(type == "meadow")  %>% 
    ggplot(aes(sample = ySim)) + 
    geom_qq() +
    geom_qq_line() + 
    xlab("QQ for meadow")
}

do.call(gridExtra::grid.arrange, plotList)

do.call(gridExtra::grid.arrange, plotListQQ1)

The most deviating QQ-plot of the residuals in the nine simulation runs:

par(mfrow=c(1,1))
plot(simModels[[8]], which = 2)

5.4 ANOVA model

mod1_anova <- anova(mod1)
mod1_anova

The p-value of the ANOVA analysis is extremely significant (p-value < 0.001), so we reject the null hypothesis that the mean egg length is equal between the different bird types. We can say that the mean egg length is significantly different between at least two bird types on the 5% significance level.

Based on this analysis, we do not yet know between which particular bird types there is a significant difference. To study this, we will perform the Tuckey post-hoc analysis.

5.5 Post-hoc analysis

We will perform a post-hoc analysis, to look at the difference in egg length between each pairwise comparison of bird types. Importantly, with this strategy, the p-values will be correctly adjusted for multiple testing.

The null hypothesis for each pairwise test is that there is no difference in the mean egg length between both bird types.

The alternative hypothesis for each pairwise test states that there is indeed a difference in the mean egg length between both bird types.

We will also calculate the confidence interval on the mean differences.

library(multcomp, quietly = TRUE)
mcp <- glht(mod1, linfct = mcp(type = "Tukey"))
summary(mcp)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = length ~ type, data = Cuckoo)
## 
## Linear Hypotheses:
##                        Estimate Std. Error t value Pr(>|t|)   
## tree - meadow == 0      0.59000    0.39263   1.503  0.65686   
## dunnock - meadow == 0   0.62143    0.40301   1.542  0.63147   
## robin - meadow == 0     0.07500    0.38332   0.196  0.99996   
## wagtail - meadow == 0   0.40333    0.39263   1.027  0.90518   
## wren - meadow == 0     -1.37000    0.39263  -3.489  0.00847 **
## dunnock - tree == 0     0.03143    0.48939   0.064  1.00000   
## robin - tree == 0      -0.51500    0.47330  -1.088  0.88202   
## wagtail - tree == 0    -0.18667    0.48087  -0.388  0.99878   
## wren - tree == 0       -1.96000    0.48087  -4.076  0.00114 **
## robin - dunnock == 0   -0.54643    0.48195  -1.134  0.86268   
## wagtail - dunnock == 0 -0.21810    0.48939  -0.446  0.99764   
## wren - dunnock == 0    -1.99143    0.48939  -4.069  0.00115 **
## wagtail - robin == 0    0.32833    0.47330   0.694  0.98170   
## wren - robin == 0      -1.44500    0.47330  -3.053  0.03198 * 
## wren - wagtail == 0    -1.77333    0.48087  -3.688  0.00433 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

confint(mcp)

## 
##   Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = length ~ type, data = Cuckoo)
## 
## Quantile = 2.8889
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##                        Estimate lwr      upr     
## tree - meadow == 0      0.59000 -0.54429  1.72429
## dunnock - meadow == 0   0.62143 -0.54285  1.78571
## robin - meadow == 0     0.07500 -1.03239  1.18239
## wagtail - meadow == 0   0.40333 -0.73096  1.53763
## wren - meadow == 0     -1.37000 -2.50429 -0.23571
## dunnock - tree == 0     0.03143 -1.38238  1.44524
## robin - tree == 0      -0.51500 -1.88234  0.85234
## wagtail - tree == 0    -0.18667 -1.57589  1.20255
## wren - tree == 0       -1.96000 -3.34922 -0.57078
## robin - dunnock == 0   -0.54643 -1.93875  0.84589
## wagtail - dunnock == 0 -0.21810 -1.63190  1.19571
## wren - dunnock == 0    -1.99143 -3.40524 -0.57762
## wagtail - robin == 0    0.32833 -1.03901  1.69567
## wren - robin == 0      -1.44500 -2.81234 -0.07766
## wren - wagtail == 0    -1.77333 -3.16255 -0.38411

5.6 Conclusion

There association between the mean length of a cuckoo’s egg and the foster bird species is an extremely significant (p << 0.001).

The mean length of cuckoo’s eggs in nests of the Eurasian wern are smaller as compared to those from all other bird species in the study:

• the meadow pipit (adjusted p-value = 0.009, mean difference = -1.4mm, 95% CI [-2.5, -0.2])
• the tree pipit (adjusted p-value = 0.001, mean difference = -2mm, 95% CI [-3.3, -0.6])
• the dunnock (adjusted p-value = 0.001, mean difference = -2mm, 95% CI [-3.4, -0.6])
• the European robin (adjusted p-value = 0.03, mean difference = -1.4mm, 95% CI [-2.8, -0.1])
• the white wagtail (adjusted p-value = 0.004, mean difference = -1.8mm, 95% CI [-3.2, -0.4])

The mean length of the cuckoo bird’s eggs for all other pairwise comparisons between foster bird species are not statistically significantly different at the 5% level.

---

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