---
title: "One-way ANOVA and Kruskal-Wallis"
author: "Emma Rand"
output:
  html_document:
    depth: 3
    highlight: pygments
    theme: lumen
    toc: yes
    toc_float:
      collapsed: no
      smooth_scroll: no
  pdf_document:
    toc: yes
  word_document: default
---
![](../../pics/17C.png)

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,
                      message = FALSE,
                      warning = FALSE)
```

# Introduction

## Aims

In this practical you will learn how to run one-way ANOVA and Kruskal-Wallis in R, interpret their output and report the results. You will also practice choosing between them.


## Learning Outcomes 

By actively following the lecture and practical and carrying out the independent study the successful student will be able to:

* Explain the rationale behind ANOVA and complete a partially filled ANOVA table (MLO 1 and 2)
* Apply (appropriately), interpret and evaluate  the legitimacy of, one-way ANOVA and Kruskal-Wallis including post-hoc tests in R (MLO 2, 3 and 4)
* Summarise and illustrate with appropriate R figures test results scientifically (MLO 3 and 4)


## Philosophy

Workshops are not a test. It is expected that you often don't know how to start, make a lot of mistakes and need help.
Do not be put off and don't let what you can not do interfere with what you can do. You will benefit from collaborating with others and/or discussing your results.

The lectures and the workshops are [closely integrated](../BIO00017C.html) and it is expected that you are familar with the lecture content before the workshop. You need not understand every detail as the workshop should build and consolidate your understanding. You may wish to refer to the slides as you work through the workshop schedule.

## Slides

One-way ANOVA and Kruskal-Wallis: [pdf](../lects/web07Oneway ANOVA and Kruskal-Wallis.pdf) (recommended) / [pptx](../lects/web07Oneway ANOVA and Kruskal-Wallis.pptx)


# Exercises

## Getting started

![W](../../pics/W.png) Start RStudio from the [Start menu](../../pics/startmenu.png).

![R](../../pics/R.png) In RStudio, [set your working directory](../../pics/menusetwd.png) to the folder you created previously for your 17C Data Analysis work.

![R](../../pics/R.png) Make a [new script](../../pics/newscript.png) file called workshop6.R to carry out the rest of the work.

![R](../../pics/R.png) You probably want to load the `tidyverse` with `library(tidyverse)`.
```{r echo=FALSE, results='hide'}
library(tidyverse)
```

## Myoglobin in seal muscle 

The myoglobin concentration of skeletal muscle of three species of seal in grams per kilogram of muscle was determined and the data are given in [seal.txt](../data/seal.txt). We want to know if there is a difference between species.  Each row represents an individual seal. The first column gives the myoglobin concentration and the second column indicates species. 

### Exploring

![R](../../pics/R.png) Save a copy to your working directory, read it in and check the structure. Do you need to check back to [Week 3: Testing, Data types and Reading in data](02TestingDataTypesReadingInData.html#import_data_2)
```{r import, echo=FALSE, results="hide"}
#---CODING ANSWER---
seal <- read.table("../data/seal.txt", header = T)
str(seal)
```


![R](../../pics/R.png) Do a quick plot of the data. Do you need to check back to [Week 7: One- and two-sample tests](05OneAndTwoSampleTests.html#egg_laying_in_parasitic_wasps)
```{r echo=FALSE, results="hide"}
ggplot(data = seal, aes(x = species, y = myoglobin)) +
  geom_violin()
```


![R](../../pics/R.png) Create a data frame called `sealsummary` that contains the means, standard deviations, sample sizes and standard errors for each species. You may want to look this up in the [lecture notes](../lects/web07Oneway%20ANOVA%20and%20Kruskal-Wallis.pdf) or last week's [workshop](05OneAndTwoSampleTests.html#egg_laying_in_parasitic_wasps#summarising_the_data)
```{r waspsummary, echo = FALSE, results='hide'}
#---CODING ANSWER---
sealsummary <- seal %>%
  group_by(species) %>%
  summarise(mean = mean(myoglobin),
            std = sd(myoglobin),
            n = length(myoglobin),
            se = std/sqrt(n))

```


### Applying, interpreting and reporting
We can now carry out a one-way ANOVA. When doing an ANOVA it is common to assign the result of the `aov()` procedure to a variable, then examine it using `summary()`

![R](../../pics/R.png) Carry out an ANOVA and examine the results with:
```{r anovatest}
mod <- aov(data = seal, myoglobin ~ species)
summary(mod)
```

Remember: the tilde (~) means test the values in myoglobin when grouped by the values in species. Or explain myoglobin with species.

![Q](../../pics/Q.png) What do you conclude so far from the test? Write your conclusion in a form suitable for a report. 
<!--
#---THINKING ANSWER---
There is a significant difference in myoglobin concentration between Seal species (ANOVA: F = 5.35; d.f. = 2, 87; p = 0.006).
-->

The ANOVA is significant but which means differ? We need a post-hoc test. A post-hoc ("after this") test is done after (and only after) a significant ANOVA test. The ANOVA tells you at least two of means differ, the post-hoc test tells you where the differences are. There are several possible post-hoc tests. A popular option is the Tukey Honest Significant Difference test.

![R](../../pics/R.png) Carry out a Tukey HSD with:

```{r posthoc}
TukeyHSD(mod)
```

Each row is a comparison between means. The 'diff' column is the difference between the means the 'p adj' column indicates whether that difference is significant.

A plot can be used to visualise the result of the post hoc and make it easier to understand

![R](../../pics/R.png) Use the generic `plot()` function to plot the post-hoc like this:
```{r poshocplot, fig.height=5}
plot(TukeyHSD(mod), cex.axis = 0.7)   #cex.axis just changes the size of the axis labels
# I could also use TukeyHSD(mod) %>% plot(cex.axis = 0.7)
```

You may want to ![zoom](../../pics/zoom.png). This shows the confidence intervals on the differences between the means and the dashed vertical line is at zero so if the confidence interval crosses that line (_i.e._, the C.I. includes zero) then the means do not differ significantly.

![Q](../../pics/Q.png) What do you conclude from the test?  
<!--
#---THINKING ANSWER---
The significant difference lies between Harbour and Bladdernose.
-->

### Checking selection

We need to examine the residuals. Unlike with `t.test` we don't have to calculate them - the object which is created by `aov()` contains a variable called `$residuals` 

![R](../../pics/R.png) Use the `shapiro.test()` on the model residuals
```{r normalityafter}
shapiro.test(mod$residuals)
```

![R](../../pics/R.png) Plot the model residuals against the fitted values like this:
```{r normalityafter2}
plot(mod, which = 1)
```
The group means are the fitted (or predicted) values; each residual is the difference between the mean and the actual value

![Q](../../pics/Q.png) What to you conclude?
<!--
#---THINKING ANSWER---
the residuals seem normally distributed and the variances look similar
-->

### Illustrating

Producing a figure to go with a significant ANOVA in a report using ggplot2

Look at the we did in the [last practical](05OneAndTwoSampleTests.html#illustrating) for egg laying behaviour of wasps - notice how the fundamental structure of the plot is the same as we require here, there is just a different number of groups. We will again use both our `seal` and `sealsummary` dataframes.


![R](../../pics/R.png) Create the plot:   
```{r sealfig, fig.width = 5, fig.height = 5}
ggplot() +
  geom_point(data = seal, aes(x = species, y = myoglobin),
             position = position_jitter(width = 0.1, height = 0),
             colour = "gray50") +
  geom_errorbar(data = sealsummary, 
                aes(x = species, ymin = mean - se, ymax = mean + se),
                width = 0.3) +
  geom_errorbar(data = sealsummary, 
                aes(x = species, ymin = mean, ymax = mean),
                width = 0.2) +
  ylab(expression("Myoglobin concentration g "*Kg^{-1})) +
  ylim(0, 80) +
  scale_x_discrete(labels = c("Bladdernose", "Harbour", "Weddell"), 
                   name = "Seal Species") +
  theme_classic()
```
   
Notice the use of `expression()` in `ylab()` to allow you to specify special characters. This figure is good but it would be nice to show the result of the post-hoc host by annotating the figure. Again, ggplot makes this fairly easy.

![R](../../pics/R.png) We can add annotations to illustrate the post-hoc results:
```{r sealfigannotated, fig.width = 5, fig.height = 5}
ggplot() +
  geom_point(data = seal, aes(x = species, y = myoglobin),
             position = position_jitter(width = 0.1, height = 0),
             colour = "gray50") +
  geom_errorbar(data = sealsummary, 
                aes(x = species, ymin = mean - se, ymax = mean + se),
                width = 0.3) +
  geom_errorbar(data = sealsummary, 
                aes(x = species, ymin = mean, ymax = mean),
                width = 0.2) +
  ylab(expression("Myoglobin concentration g "*Kg^{-1})) +
  ylim(0, 80) +
  scale_x_discrete(labels = c("Bladdernose", "Harbour", "Weddell"), 
                   name = "Seal Species") +
  annotate("segment", x = 1, xend = 2, 
           y = 72, yend = 72,
           colour = "black") +
  annotate("segment", x = 2, xend = 2, 
           y = 72, yend = 70,
           colour = "black") +
  annotate("segment", x = 1, xend = 1, 
           y = 72, yend = 70,
           colour = "black") +
  annotate("text", x = 1.5,  y = 74, 
           label = "*", size = 8) +
  theme_classic()
  
```

Figure 1. The mean myoglobin concentration of three seal species. Error bars are $\pm$ 1 S.E. 

![Q](../../pics/Q.png) Write up the result is a form suitable for a report.
<!-- There is a significant difference in myoglobin concentration between Seal species (ANOVA: F = 5.35; d.f. = 2, 87; p = 0.006). Post-hoc testing revealed that difference to be between the Harbour Seal with the highest myoglobin concentrations (49.0 +/- 1.51)` ) and the Bladdernose Seal with the lowest (42.3 +/- 1.46). -->

## Leafminers on Birch

Larvae of the Ambermarked birch leafminer, _Profenusa thomsoni_, feed on the interior leaf tissues of Birch (Betula) species. They do not normally kill the tree but can weaken it making it susceptible to attack from other species. Researchers are interested in whether there is a difference in the rates at which white, grey and yellow birch are attacked. They introduce adult female _P.thomsoni_ to a green house containing 30 young trees (ten of each type) and later count the egg laying events on each tree. The data are in [leaf.txt](../data/leaf.txt).

### Exploring

![R](../../pics/R.png) Save a copy to your working directory, read it in and check the structure
```{r echo = FALSE, results = "hide"}
#---CODING ANSWER---
leaf <- read.table("../data/leaf.txt", header = T)
str(leaf)
```

![R](../../pics/R.png) Calculate the medians, means and sample sizes. 
```{r echo = FALSE,results='hide'}
#---CODING ANSWER---
leaf %>% 
  group_by(birch) %>% 
  summarise(mean = mean(eggs),
            median = median(eggs),
            n = length(eggs))
```

![Q](../../pics/Q.png) Why is a Kruskal-Wallis appropriate in this case?
<!--
#---THINKING ANSWER---
the data are counts hence they are not continuous. We can sometimes treat count data as continuous if there is a large number of possible values and the distance between each value relative to the range is small but this is not the case here. In addition, with only ten values in each group it is impossible to judge whether the data are normally distributed
-->

### Applying, interpreting and reporting

To carry out a Kruskal-Wallis we use:
```{r}
kruskal.test(data = leaf, eggs ~ birch)
```

![Q](../../pics/Q.png) What do you conclude from the test? Write your conclusion in a form suitable for a report. 
<!--
#---THINKING ANSWER---
There is a significant difference between birch types in the number of laying events by _P.thomsoni_ (Kruskal-Wallis: chi-squared = 6.34; d.f. = 2; p = 0.042).
-->

The significant Kruskal-Wallis tells us at least two of the groups differ but where do the differences lie? A post-hoc multiple comparison test for a significant Kruskal-Wallis exists in the pgirmess package.

![R](../../pics/R.png) Load the package using `library()`
```{r echo = FALSE, results = "hide", warning = FALSE, message = FALSE}
library(pgirmess)
```
Now run the test with:
```{r}
kruskalmc(data = leaf, eggs ~ birch)
```
![Q](../../pics/Q.png) What do you conclude from the test? 
   <!--
#---THINKING ANSWER---
The signifcant difference is between yellow birch and white birch (the highest and lowest number of eggs)
-->

![Q](../../pics/Q.png) Write up the result is a form suitable for a report.
<!-- There is a significant difference in the rate of leafminer attack between birch species (Kruskal-Wallis: chi-squared = 6.34, df = 2, p-value = 0.042). Post-hoc testing revealed that difference to be between the White Birch with the lowest number of eggs (median = 20) and Yellow Birch with the highest (median = 54). -->

### Illustrating

![R](../../pics/R.png) A box plot is an appropriate choice for illustrating a Kruskal-Wallis. Can you produce a figure like this?

Tip - Look at the grouse-nematode figure we did in the [last practical](05OneAndTwoSampleTests.html#illustrating17) for two-sample Wilcoxon.

```{r echo = FALSE, fig.width = 5, fig.height = 5}
#---CODING ANSWER---
ggplot(leaf, aes(x = birch, y = eggs) ) +
  geom_boxplot() +
  xlab("Birch") +
  ylab("Number of eggs") +
  ylim(0, 105) +
  annotate("segment", x = 2, xend = 3, 
           y = 100, yend = 100,
           colour = "black") +
  annotate("segment", x = 2, xend = 2, 
           y = 100, yend = 97,
           colour = "black") +
  annotate("segment", x = 3, xend = 3, 
           y = 100, yend = 97,
           colour = "black") +
  annotate("text", x = 2.5,  y = 102, 
           label = "*", size = 8) +
  theme_classic()
```

## Saving figures

A useful function for save figures is `ggsave()`. It has arguments for the size of the image and resolution. See the [`ggsave()` reference page](https://ggplot2.tidyverse.org/reference/ggsave.html). Since I often make more than one figure, I might set these arguments first:

```{r}
# figure saving settings
units = "in"
fig_w <- 3.2
fig_h <- fig_w
dpi <- 300
device <- "tiff" # this is format often required by journals; you may want png or jpg
```

Then I would normally assign the ggplot to a variable to pass to the `ggsave()` function.

```{r}
fig1 <- ggplot(leaf, aes(x = birch, y = eggs) ) +
  geom_boxplot() +
  xlab("Birch") +
  ylab("Number of eggs") +
  ylim(0, 105) +
  annotate("segment", x = 2, xend = 3, 
           y = 100, yend = 100,
           colour = "black") +
  annotate("segment", x = 2, xend = 2, 
           y = 100, yend = 97,
           colour = "black") +
  annotate("segment", x = 3, xend = 3, 
           y = 100, yend = 97,
           colour = "black") +
  annotate("text", x = 2.5,  y = 102, 
           label = "*", size = 8) +
  theme_classic()


ggsave("fig1-birch.tif", 
       plot = fig1, 
       device = device,
       width = fig_w, 
       height = fig_w,
       units = units,
       dpi = dpi)
```


# Independent study

Decide which test you should use to analyse the following data sets. In each case give the reasons for your choice of test and state the null hypothesis. Write your conclusions in a form suitable for including in a report. Can you make figures?

#### Effects of fitness and heat acclimatisation

Sports scientists were investigating the effects of fitness and heat acclimatisation on the sodium content of sweat. They measured the sodium content of the sweat ($\mu mol\mspace{3mu}l^{-1}$) of three groups of individuals: unfit and unacclimatised (UU); fit and unacclimatised(FU); and fit and acclimatised (FA). The are in [sweat.txt](../data/sweat.txt). Is there a difference between the groups in the sodium content of their sweat?
   
```{r include=FALSE}
#---CODING AND THINKING ANSWER---
#read in the data and look at structure
sweat <- read.table("../data/sweat.txt", header = T)
str(sweat)

# quick plot of the data
ggplot(data = sweat, aes(x = gp, y = na)) +
  geom_boxplot()
ggplot(data = sweat, aes(x = gp, y = na)) +
  geom_point()

# Since the sample sizes are small and not the same in each group and the 
# variance in the FA gp looks a bit lower, I'm leaning to a non-parametric test K-W.
# However, don't panic if you decided to do an anova

# calculate some summary stats 
sweatsum <- sweat %>% 
  group_by(gp) %>% 
  summarise(mean = mean(na),
            std = sd(na),
            n = length(na),
            median = median(na))


# Kruskal-Wallis
kruskal.test(data = sweat, na ~ gp)
# We can say there is a difference between the groups in the sodium 
# content of their sweat (chi-squared = 11.9802, df = 2, p-value = 0.002503).
# Unfit and unacclimatised people have most salty sweat, 
# Fit and acclimatised people the least salty.

# a post-hoc test to see where the sig differences lie:

kruskalmc(data = sweat, na ~ gp)
# Fit and acclimatised people (FA) have significantly less sodium in their
#  sweat than the unfit and unacclimatised people (UU). 
# Fit and unacclimatised (FU) people have sodium concentrations 
# more similar to the FA group but don't reach significance 
# for being different to UU. See figure 1.

ggplot(sweat, aes(x = gp, y = na) ) +
  geom_boxplot() +
  xlab("Group") +
  ylab(expression("Sodium"*mu*"mol"*l^{-1})) +
  scale_x_discrete(labels = c("Fit Acclimatised", 
                              "Fit Unacclimatised", 
                              "Unfit Unacclimatised"), 
                   name = "") +
  ylim(0, 100) +
  annotate("segment", x = 1, xend = 3, 
           y = 90, yend = 90,
           colour = "black") +
  annotate("segment", x = 1, xend = 1, 
           y = 90, yend = 87,
           colour = "black") +
  annotate("segment", x = 3, xend = 3, 
           y = 90, yend = 87,
           colour = "black") +
  annotate("text", x = 2,  y = 93, 
           label = "**", size = 8) +
  theme_classic()

#Figure 1. Sodium content of sweat for three groups: Fit and acclimatised
#(FA), Fit and unacclimatised (FU) and Unfit and unacclimatised (UU). Heavy lines
#indicate the median, boxes the interquartile range and whiskers the range. 

```

#### Insecticide effectiveness
The data are given in [biomass.txt](../data/biomass.txt) are taken from an experiment in which the insect pest biomass (g) was measured on plots sprayed with water (control) or one of five different insecticides. Do the insecticides vary in their effectiveness? What advice would you give to a person:
 - currently using insecticide E?
 - trying to choose between A and D?
 - trying to choose between C and B?
   
```{r include=FALSE}
#---CODING AND THINKING ANSWER---
######################################################################
#                                                                    #
#   A comparison of the effects of five insect pesticides on the     #
#   insect biomass in treated plots.                                 #
#                                                                    #     
######################################################################

######################################################################
#                             Introduction                           #
######################################################################

# The data are given in biomass.txt are taken from an experiment 
# in which the insect pest biomass (g) was measured on plots sprayed 
# with water (control) or one of five different insecticides. 
# The goal of the analysis was to determine if the insecticides 
# vary in their effectiveness and specifically advise on:
#   - use of insecticide E
#   - the choice between A and D
#   - the choice between C and B

# The data are organised with an insecticide treatment group in
# each column:
# 'data.frame':	10 obs. of  6 variables:
# $ WaterControl: num  350 324 359 255 208 ...
# $ A           : num  159 146 116 135 137 ...
# $ B           : num  150.1 154.4 69.5 150.7 212.6 ...
# $ C           : num  80 266.4 161.2 161.4 51.2 ...
# $ D           : num  267 110 221 160 198 ...
# $ E           : num  350 320 359 255 208 ...

######################################################################
#                       Import and tidy data                         #
######################################################################

# data are in ../data
biom <- read.table("../data/biomass.txt", header = T)

# check structure
str(biom)

# 'data.frame':	10 obs. of  6 variables:
# $ WaterControl: num  350 324 359 255 208 ...
# $ A           : num  159 146 116 135 137 ...
# $ B           : num  150.1 154.4 69.5 150.7 212.6 ...
# $ C           : num  80 266.4 161.2 161.4 51.2 ...
# $ D           : num  267 110 221 160 198 ...
# $ E           : num  350 320 359 255 208 ...

# The data are organised with an insecticide treatment group in
# each column. Put the data into tidy format.

biom <- gather(biom, key = spray, value = biomass)

######################################################################
#                     Exploratory Analysis                           #
######################################################################

# quick plot of the data
ggplot(data = biom, aes(x = spray, y = biomass)) +
  geom_boxplot()


# summary statistics
biomsum <- biom %>% 
  group_by(spray) %>% 
  summarise(mean = mean(biomass),
            median = median(biomass),
            sd = sd(biomass),
            n = length(biomass),
            se = sd / sqrt(n))

# conclusion: the sample sizes are equal, 10 is a smallish but
# reasonable sample size
# the means and medians are similar to each other (expected for
# normally distributed data), A has a smaller variance 

# We have one explanatory variable, "spray" comprising 6 levels
# Biomass has decimal places and we would expect such data to be 
# normally distributed therefore one-way ANOVA is the desired test
# - we will check the assumptions after building the model

######################################################################
#                     Statistical Analysis                           #
######################################################################

# Carrying out an ANOVA
model <- aov(data = biom, biomass ~ spray)
summary(model)
# There is a very highly signifcant effect of spray identity on pest 
# biomass (F = 26.5; d.f., 5, 54; p < 0.001).

# Carrying out a Tukey Honest Signifcant differences test
# to see where differences bewteen spray treatments lie
# ordering can make it easier to understand. Plot included 
TukeyHSD(model, ordered = T)
plot(TukeyHSD(model, ordered = T), cex.axis = 0.5)

# signifcant comparisions
#                   diff        lwr       upr     p adj
# D-A             76.505  11.729841 141.28016 0.0118570
# E-A            175.515 110.739841 240.29016 0.0000000
# WaterControl-A 175.915 111.139841 240.69016 0.0000000
# E-C            155.710  90.934841 220.48516 0.0000000
# WaterControl-C 156.110  91.334841 220.88516 0.0000000
# E-B            154.323  89.547841 219.09816 0.0000001
# WaterControl-B 154.723  89.947841 219.49816 0.0000000
# E-D             99.010  34.234841 163.78516 0.0004759
# WaterControl-D  99.410  34.634841 164.18516 0.0004477

# Note smaller values (insect biomass) indicates more 
# effective control

# All sprays are better than the water control except E. 
# This is probably the most important result.
# What advice would you give to a person currently using insecticide E?
# Don't bother!! It's no better than water. Switch to any of 
# the other sprays
# 

# Sorting the summary table by the mean can make it easier to 
# interpret the results
arrange(biomsum, mean)
# What advice would you give to a person currently
#   + trying to choose between A and D? Choose A because A has sig lower
#   insect biomass than D 
#   + trying to choose between C and B? It doesn't matter because there is 
#   no differnce in insect bbiomass. Use other criteria to chose (e.g., price)

# We might report this like:
# There is a very highly signifcant effect of spray identity on pest 
# biomass (F = 26.5; d.f., 5, 54; p < 0.001). Post-hoc testing 
# showed E was no more effective than the control; A, C and B were 
# all better than the control but could be equally as good as each
# other; D would be a better choice than the control or E but 
# worse than A. See figure 1


######################################################################
#                                 Figure                             #
######################################################################

# uses the summary data
# I reordered the bars to make is easier for me to annotate with

ggplot() +
  geom_point(data = biom, aes(x = reorder(spray, biomass), y = biomass),
             position = position_jitter(width = 0.1, height = 0),
             colour = "gray50") +
  geom_errorbar(data = biomsum, 
                aes(x = spray, ymin = mean - se, ymax = mean + se),
                width = 0.3) +
  geom_errorbar(data = biomsum, 
                aes(x = spray, ymin = mean, ymax = mean),
                width = 0.2) +
  ylim(0, 520) +
  ylab("Pest Biomass (units)") +
  xlab("Spray treatment") +
  # E and control are one group
  annotate("segment", x = 4.5, xend = 6.5, 
           y = 397, yend = 397,
           colour = "black", size = 1) +
  annotate("text", x = 5.5,  y = 385, 
           label = "N.S", size = 4) +
  # WaterControl-D and E-D    ***
  annotate("segment", x = 4, xend = 5.5, 
           y = 410, yend = 410,
           colour = "black") +
  annotate("segment", x = 4, xend = 4, 
           y = 410, yend = 400,
           colour = "black") +
  annotate("segment", x = 5.5, xend = 5.5, 
           y = 410, yend = 400,
           colour = "black") +
  annotate("text", x = 4.5,  y = 420, 
           label = "***", size = 5) +
  # WaterControl-B ***
  annotate("segment", x = 3, xend = 5.5, 
         y = 440, yend = 440,
         colour = "black") +
  annotate("segment", x = 3, xend = 3, 
           y = 440, yend = 430,
           colour = "black") +
  annotate("segment", x = 5.5, xend = 5.5, 
           y = 440, yend = 430,
           colour = "black") +
  annotate("text", x = 4,  y = 450,
           label = "***", size = 5) +
  # WaterControl-C ***
  annotate("segment", x = 2, xend = 5.5, 
           y = 475, yend = 475,
           colour = "black") +
  annotate("segment", x = 2, xend = 2, 
           y = 475, yend = 465,
           colour = "black") +
  annotate("segment", x = 5.5, xend = 5.5, 
           y = 475, yend = 465,
           colour = "black") +
  annotate("text", x = 3.5,  y = 485, 
           label = "***", size = 5) +
  # WaterControl-A ***
  annotate("segment", x = 1, xend = 5.5, 
         y = 510, yend = 510,
         colour = "black") +
  annotate("segment", x = 1, xend = 1, 
           y = 510, yend = 500,
           colour = "black") +
  annotate("segment", x = 5.5, xend = 5.5, 
           y = 510, yend = 500,
           colour = "black") +
  annotate("text", x = 3.5,  y = 520, 
           label = "***", size = 5) +  
# A-D ***
  annotate("segment", x = 1, xend = 4, 
         y = 330, yend = 330,
         colour = "black") +
  annotate("segment", x = 1, xend = 1, 
           y = 330, yend = 320,
           colour = "black") +
  annotate("segment", x = 4, xend = 4, 
           y = 330, yend = 320,
           colour = "black") +
  annotate("text", x = 2.5,  y = 335, 
           label = "*", size = 5) +
  theme_classic()

# Figure 1. The mean pest biomass following various insecticide treatments.
# Error bars are +/- 1 S.E. Significant comparisons are indicated.
```


# The Code files

These contain answers and code even though they do not appear on the webpage itself.

[Rmd file](06One-wayAnovaAndKW.Rmd)
The Rmd file is the file I use to compile the practical. Rmd stands for R markdown allow R code and ordinary text to be inter weaved to produce well-formatted reports including webpages.

[Plain script file](../scripts/06One-wayAnovaAndKW.R)
This is plain script (.R) version of the practical


# Script example
This is an example of a well formatted analysis script for one of the independent study problems.

[Script example](../scripts/exampleprac6.R)

# Objectives from previous sessions

[Introduction to module and RStudio](01IntroductionToModuleAndRStudio.html)

* to explain why we need statistical tests and the logic of hypothesis testing (MLO 1)
* use the R command line as a calculator and to assign variables (MLO 3)
* create and use the basic data types in R (MLO 3)
* find their way around the RStudio windows (MLO 3)
* create, use and save a script file to run r commands (MLO 3)
* search and understand manual pages (MLO 3)

[Testing, Data types and reading in data](02TestingDataTypesReadingInData.html)

* to able to explain what response and explanatory variables are, distinguish between data types and describe how these impact choice of test (MLO 1 and 2)
* demonstrate the process of hypothesis testing with an example and evaluate potential inferences (MLO 1 and 2)
* read in data in to RStudio, create simple summaries and plots using manual pages where necessary (MLO 3)
* create neat reports in Word which include text and figures (MLO 4)

[Goodness of Fit and Contingency chi-squared tests](03ChiSquaredTests.html)

* recognise when to use chi-squared Goodness of Fit and Contingency tests (MLO 2)
* be able to carry out, interpret and report scientifically both types of test by hand and in R (MLO 3 and 4)

[Calculating summary statistics, probabilities and confidence intervals](04CalculatingSummaryStatisticsProbabilitiesAndCI.html)

* Explain the properties of 'normal distributions' and their use in statistics (MLO 1 and 2)
* Define, select and calculate with R probabilities, quantiles and confidence intervals (MLO 3 and 4)

[One- and two- sample tests](05OneAndTwoSampleTests.html)

* Explain dependent and independent samples (MLO 2)
* Select, appropriately, t-tests and their non-parametric equivalents (MLO 2)
* Apply, interpret and evaluate the legitimacy of  the tests in R (MLO 3 and 4)
* Summarise and illustrate with appropriate R figures test results scientifically (MLO 3 and 4)

![](../../pics/17Cend.png)
