---
title: "Correlation and Regression"
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/06C.png)


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


# Introduction

## Aims

You will learn how to chose between correlation and regression as well as applying, interpreting and reporting 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 principles of correlation and of regression (MLO 1)
* Apply (appropriately), interpret and evaluate the legitimacy of, both 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

Correlation and Regression: [pdf](../lects/web09Correlation and Regression.pdf) (recommended) / [pptx](../lects/web09Correlation and Regression.pptx)


# Exercises

## Getting started

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

![R](../../pics/R.png) Make a new project with File | New Project and chose New directory and then New project. Be purposeful about where you create it by using the Browse button. I suggest using your 17C folder. Give the Project (directory) a name, perhaps "regress_correl"

![R](../../pics/R.png) Make a new folder 'raw_data' where you will later save data files.

![R](../../pics/R.png) Make a new folder 'figures' where you will later save your figures.

![R](../../pics/R.png) Make a [new script](../../pics/newscript.png) file called analysis.R or similar 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)
```



## Pearson's Correlation 

The data given in [height.txt](../data/height.txt) are the heights of eleven brother and sister pairs. 

![R](../../pics/R.png) Save a copy of [height.txt](../data/height.txt) to your `raw_data` folder and import it.
    
```{r echo=FALSE}
#---CODING ANSWER---
# note: your directory structure differes from mine so you'll  have "raw_data/height.txt"
height <- read.table("../data/height.txt", header = TRUE)

```

### Exploring

![Q](../../pics/Q.png) What type of variables are 'brother' and  'sister'? What are the implications for the test? 
<!--
#---THINKING ANSWER---
They have decimal places and height is a variable we would expect to be normally distributed. It's probable that a parametric correlation is fine. However they are quite small samples, making it difficult to show that so doing a non-parametric would not be especially over cautious. 
-->


![Q](../../pics/Q.png) Do a quick plot of the data. We don't have a causal relationship here so either varaible can go on the *x*-axis.
   
```{r,fig.width=4,fig.height=4}
ggplot(height, aes(x = sister, y = brother) ) +
  geom_point()
  
```

![Q](../../pics/Q.png) Remembering that one of the assumptions for parametric correlation is that any correlation should be linear, what do you conclude from the plot? 
<!--
#---THINKING ANSWER---
I'm not totally convinced these look linear - but then I wouldn't be convinced of any other relationship either. There are not always 'right' answers to these questions - sometimes we have to make a judgement .
--> 

### Applying, interpreting and reporting

We will do a parametric correlation in any case. 

![R](../../pics/R.png) We can carry out a Pearson's product moment correlation with:
```{r}
cor.test(height$sister, height$brother, method = "pearson")
```
  
![Q](../../pics/Q.png) What do you conclude from the test?  
<!--
#---THINKING ANSWER---
There is no significant correlation between the heights of brothers and sisters - the p value is greater than 0.05
 -->

### Illustrating

![R](../../pics/R.png) Create a better figure for our data using:
```{r,fig.width=4,fig.height=4}
fig1 <- ggplot(height, aes(x = sister, y = brother)) +
  geom_point() +
  xlim(120, 180) +
  ylim(120, 190) +
  xlab("Heights of sister (cm)") +
  ylab("Heights of brother (cm)") +
  theme_classic()

fig1
```

Figure 1. Correlation in height between brothers and sisters.

![R](../../pics/R.png) Use `ggsave()` to save your figure to file. You can use what ever format you prefer (png, jpg, tiff eps). You may want to look up [a previous week's workshop](06One-wayAnovaAndKW.html).

```{r echo = FALSE}
# you will need figures/figure1.png to save in your figures directory
# i choose to do a png but you may have a different preference
ggsave("periwinkle.png",
       plot = fig1,
       width = 3.5,
       height = 3.5,
       units = "in",
       dpi = 300)

```




## Effect of sample size on correlation

Now we will explore the effect of sample size on the value of the correlation coefficient and its significance. 

![R](../../pics/R.png) Create a dataset with twice the number of observation like this: 

```{r}
height2 <- rbind(height, height)
str(height2)
```

Each pair of values will apear twice.

![R](../../pics/R.png) Now repeat the correlation with height2
```{r echo=FALSE}
#---CODING ANSWER---
  cor.test(height2$sister, height2$brother, method = "pearson")
```
   
![Q](../../pics/Q.png) What do you conclude? What does this tell you about the sensitivity of correlation to sample size? 
<!--
#---THINKING ANSWER---
The correlation coefficient is the same but now it's significant!
When you are carrying out correlations you should take note of the sample size. Often you see people carrying out correlations on very large datasets and, for example finding a correlation coefficient of r = 0.03 which is highly significant and drawing the conclusion that the correlation is very strong/meaningful. It isn't. Yes it is highly statistically significant, but that doesn't necessarily mean it is biologically significant. You need to both consider both the r value and the p value - especially with large datasets.
--> 

## Spearman's rank Correlation

Since our brother-sister dataset is so small we might very reasonably have chosen to do a non-parametric correlation. This is very easy to do.

![R](../../pics/R.png) We just need to change the method: 

```{r}
cor.test(height$sister, height$brother, method = "spearman")
```
  ![Q](../../pics/Q.png) What do you conclude?  
<!--
#---THINKING ANSWER---
the correlation coefficient is different (r = 0.501), as is the p value (0.116). Still not significant - as we would expect as non-parametric tests are (usually) more conservative than parametric tests.
-->

## Linear Regression

The data in [plant.xlsx](../data/plant.xlsx) is a set of observations of plant growth over two months. The researchers planted the seeds and harvested, dried and weighed a plant each day from day 10 so all the data points are independent of each other.
There's a good package for reading in Excel files, `readxl`. It is one of the [tidyverse](https://www.tidyverse.org/) packages and will be installed when you install `tidyverse` but it is not one of the core packages that get loaded with `library(tidyverse)` therefore we need to library it. I recommend putting all your library statements together at the top of the file.

```{r}
library(readxl)
```

![R](../../pics/R.png) Save a copy of [plant.xlsx](../data/plant.xlsx) to your `raw_data` folder and import it.
    
```{r echo=FALSE}
#---CODING ANSWER---
# note: your directory structure differes from mine so you'll  have "raw_data/height.txt"
height <- read.table("../data/plant.xlsx", header = TRUE)

```


Excel workbooks can have multiple sheets. You can list the sheets in the workbook using `excel_sheets()`

![R](../../pics/R.png) What sheets are there in [plant.xlsx](../data/plant.xlsx)?
```{r eval=FALSE}
excel_sheets("raw_data/plant.xlsx")
```

```{r echo=FALSE}
## note: your directory structure differs from mine so you'll  have "raw_data/plant.txt"
excel_sheets("../data/plant.xlsx")
```


![R](../../pics/R.png) To read read the data in we use the `read_excel()` function
```{r eval=FALSE}

plant <- read_excel("raw_data/plant.xlsx", sheet = "plant")

```

```{r echo = FALSE}

## note: your directory structure differes from mine so you'll  have "raw_data/plant.txt"
plant <- read_excel("../data/plant.xlsx", sheet = "plant")

```
By default, the `read_excel()` function will read the first sheet in the workbook but the sheet option allows us to specify a particular sheet.

  
![Q](../../pics/Q.png) What type of variables do you have? Which is the response and which is the explanatory? What is the null hypothesis?
<!--
#---THINKING ANSWER---
Day is discrete (but ordered) and is the explanatory; mass looks continuous and is the response. The null hypothesis can be expressed as one of the following: day does not explain the variation in mass; there is no linear relationship between day and mass; slope of the best fitting straight line is zero
-->



### Exploring

![R](../../pics/R.png) Do a quick plot of the data (as shown)
```{r echo=FALSE}
#---CODING ANSWER---
ggplot(plant, aes(x = day, y = mass)) +
  geom_point()
```

![Q](../../pics/Q.png) What are the assumptions of linear regression? Do these seem to be met?
<!--
#---THINKING ANSWER---
linear relationship - looks ok on the scatter
each y is drawn from a normal distribution for each x and these normal distributions have the same variance - looks ok too. We'd expect y to be continuous and the scatter around a best fitting line we imagine is about even.
The values of y are independent of each other - yes it was a different plant each day (in fact plant was destroyed). Measuring the height of the SAME plants each day would violate this assumption
 x has been chosen by the experimenter - yes.
-->

### Applying, interpreting and reporting

![R](../../pics/R.png) We now carry out a regression assigning the result of the `lm()` procedure to a variable and examining it with `summary()`.
```{r}
mod <- lm(data = plant, mass ~ day)
summary(mod)
```

The Estimates in the Coefficients table give the intercept (first line) and the slope (second line) of the best fitting straight line. The _p_-values in the same table are tests of whether that coefficient is different from zero. 

The _F_ value and _p_-value in the last line are a test of whether the model as a whole explains a significant amount of variation in the dependent variable. For a single linear regression this is exactly equivalent to the test of the slope against zero.

![Q](../../pics/Q.png) What is the equation of the line? What do you conclude from the analysis? 
<!--
#---THINKING ANSWER---
 mass = 1.606*day - 8.683
 the slope of the line is significantly different from zero / the day explains a significant amount of the variation in mass (ANOVA: F = 88.4; d.f. = 1,49; p < 0.0001). 
 -->
 
![Q](../../pics/Q.png) Does the line go through (0,0)? 
<!--
#---THINKING ANSWER---
yes it could - whilst the intercept is -8.683 it is not significantly different from zero thus it could be zero 
-->

![Q](../../pics/Q.png) What percentage of variation is explained by the line? 
<!--
#---THINKING ANSWER---
64% 
-->


### Checking selection

![R](../../pics/R.png) Check the assumptions of the test by looking at the distribution of the 'residuals'.
```{r assumptionperi, include=FALSE}
#---CODING AND THINKING ANSWER---
plot(mod, which = 1)
hist(mod$residuals)
shapiro.test(mod$residuals)
#These look ok
```

### Illustrating
We want a figure with the points and a best fitting straight line.

![R](../../pics/R.png) Create a figure with using both `geom_point()` and `geom_smooth()`

 
```{r,fig.width=4,fig.height=4}
ggplot(plant, aes(x = day, y = mass)) +
  geom_point() +
  geom_smooth(method = lm, se = FALSE, colour = "black") +
  ylab("Mass (g)") +
  ylim(0, 120) +
  xlim(0, 65) +
  xlab("Day") +
  theme_classic()
  
```

Figure 2. The relationship between day since planting and mass. 

![R](../../pics/R.png) Can you workout how to add the equation of the line to the figure?

```{r, echo = FALSE, fig.width=4,fig.height=4}
fig2 <- ggplot(plant, aes(x = day, y = mass)) +
  geom_point() +
  geom_smooth(method = lm, se = FALSE, colour = "black") +
  ylab("Mass (g)") +
  ylim(0, 120) +
  xlim(0, 65) +
  xlab("Day") +
  annotate("text", x = 15, y = 100, label = expression(italic("y = 1.6 x - 8.69"))) +
  theme_classic()
  
```


![R](../../pics/R.png) Save your figure to your 'figures' folder.
```{r echo = FALSE}
# you will need figures/figure1.png to save in your figures directory
# i choose to do a png but you may have a different preference
ggsave("fig2.png",
       plot = fig2,
       width = 3.5,
       height = 3.5,
       units = "in",
       dpi = 300)

```

![W](../../pics/W.png) Close your project, locate the folder in Windows Explorer and zip it by doing: Right click, Send to | Compressed (zipped) folder. Then email it to your neighbour and check they can open your unzipped project and have all the code work.


# Independent study


## Analyses

Decide how to analyse the following data sets. In each, use an RStudio project with logical directory structure. Organise and comment your code well and include your reasoning and decisions. Write your conclusions as comments in a form suitable for including in a report. Create and save an appropriate figure.

### Effect of anxiety status and sporting performance
The data in [sprint.txt](../data/sprint.txt) are from an investigation of the effect of anxiety status and sporting performance. A group of 40 100m sprinters undertook a psychometric test to measure their anxiety shortly before competing. The data are their anxiety scores and the 100m times achieved. What you do conclude from these data?


```{r include=FALSE}
#---CODING AND THINKING ANSWER---
# this example is designed to emphasise the importance of plotting your data first
sprint <- read.table("../data/sprint.txt", header = TRUE)
# Anxiety is discrete but ranges from 16 to 402 meaning the gap between possible measures is small and 
# the variable could be treated as continuous if needed. Time is a contrinuous measure that has decimal places and which we would expect to follow a normal distribution 

# explore with a plot
ggplot(sprint, aes(x = anxiety, y = time) ) +
  geom_point()

# A scatterplot of the data clearly reveals that these data are not linear. There is a good relationship between the two variables but since it is not linear, a correlation (either parametric or non-parametric) is inappropriate. A Pearson's correlation comes out at close to zero and NS (since the scatter around a linear correlation would be great).

```  


### Juvenile hormone in stag beetles
The concentration of juvenile hormone in stag beetles is known to influence mandible growth. Groups of stag beetles were injected with different concentrations of juvenile hormone (arbitrary units) and their average mandible size (mm) determined. The experimenters planned to analyse their data with regression. The data are in [stag.txt](../data/stag.txt) 


```{r include=FALSE}
#---CODING AND THINKING ANSWER---
#read the data in and check the structure
stag <- read.table("../data/stag.txt",header=T)
str(stag)

# jh is discrete but order and has been chose by the experimeter - it is the explanatory variable.  
# the response is mandible size which has decimal placesand is something we would expect to be 
# normally distrubted. So far, common sense suggests the assumptions of regression are met.

# exploratory plot
ggplot(stag, aes(x = jh, y = mand)) +
  geom_point()
# looks linear-ish on the scatter
# regression still seems appropriate
# we will check the other assumptions after we have run the lm

# build the statistical model
mod <- lm(data = stag, mand ~ jh)

# examine it
summary(mod)
# mand = 0.032*jh + 0.419
# the slope of the line is significantly different from zero / the jh explains a significant amount of the variation in mand (ANOVA: F = 16.63; d.f. = 1,14; p = 0.00113).
# the intercept is 0.419 and differs significantly from zero 

# checking the assumption
plot(mod, which = 1) 
# we're looking for the variance in the residuals to be equal along the x axis.
# with a small data set there is some apparent heterogeniety but it doesn't look too.
# 
hist(mod$residuals)
# We have some skew which gain might be partly a rewsult od a small sample size.
shapiro.test(mod$residuals) # the also test not sig diff from normal

# On balance the use of regression is probably justifiable
# but ideally the experiment would be better if multiple individauls were measure at
# each of the chosen juvenile hormone levels.

# a better plot
ggplot(stag, aes(x = jh, y = mand) ) +
  geom_point() +
  geom_smooth(method = lm, se = FALSE, colour = "black")+
  ylab("Mandible size (mm)") +
  ylim(0, 2) +
  xlim(0, 32) +
  xlab("Juvenile hormone (arbitrary units)") +
  theme_classic()
#Figure x. The effect of juvenile hormone injections on the mandible size of stag beetles.   
```

# The Code files

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

[Rmd file](08CorrelationAndRegression.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/08CorrelationAndRegression.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/exampleprac8.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)

[One-way ANOVA and Kruskal-Wallis](06One-wayAnovaAndKW.html)

* 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)

[Two-way ANOVA](07Two-wayAnova.html)

* Explain the rationale behind ANOVA and complete a partially filled ANOVA table (MLO 1 and 4)
* Read in data formatted for other statistical packages (MLO 3)
* Apply (appropriately), interpret and evaluate the legitimacy of, two-way ANOVA in R (MLO 2, 3 and 4)
* Explain the meaning of a significant interaction (MLO 4)
* Summarise and illustrate with appropriate figures test results scientifically (MLO 3 and 4)
* Use RStudio projects (MLO 4)


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