# Returning to the sleep data analysis


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

```

```{r, message=FALSE}
library(tidyverse)
library(aqm)
library(plotly)
```

## Introduction

The sleep data set is used for the assignment. Your task is to come up with some meaningful scientific questions that can be addressed through applying graphical and statistical analysis. To help with this we can look at the relationship between brain weight and body weight and derive an additional variable from doing so.

## Read in the data

The data are included in the AQM package. This makes it simple to begin the analysis.

```{r}
data(sleep)
d<-sleep

```



For most analyses it is preferable to make sure that any categorical variables are considered to be factors. The code below does this. 

```{r}
d$Diet<- as.factor(d$Diet)
d$Species<- as.factor(d$Species)
```

Now we already know that these variables need log transforming. Make sure that you do this, as analysing the raw brain and body weights will not make much sense from a statistical perspective

```{r}
d$LogBrainWt <- log10(d$BrainWt)
d$LogBodyWt <- log10(d$BodyWt)
```

## Using ggplot for visualising regression

```{r}
g0 <- ggplot(d, aes(x=LogBodyWt, y=LogBrainWt, label=Species)) # set up the aesthetics
g1 <- g0 + geom_point() + geom_smooth(method="lm") # Add the geometries
g1<-g1 +labs(x="Log10 body weight measured in grams", y="Log10 brain weight measured in grams")
g1
```

We can label the points, but this is going to look very messy.

```{r}
g1 + geom_label()
```

## Using plotly for interactive data exploration

Although this cannot be used in a publication, interactive data exploration using plotly is the easiest way to solve this problem when you are exploring the data yourself.


```{r}
ggplotly(g1)
```

## Fitting the regression

We can fit the regression in the usual manner.

```{r}
mod <- lm(data=d, LogBrainWt~LogBodyWt)

```

Note the warning about not being interested in elements such as p-values. Although we can ask R for a summary of the model we have no interest in the significance. The $R^2$ may be of some limited interest. However the real reason for fitting the model is to obtain the residuals in this case.

```{r}
summary(mod)
```

## Obtaining the residuals

The residuals can be used in additional analyses as an indicator of brain weight **after controlling for body weight**. Note that they remain on a logarithmic scale.

```{r}
d$brain_index <- residuals(mod)
```

## Some tricks to plot the residuals

The fct_reorder function that is part of the tidyverse can be handy for reordering a factor.

```{r}
d$Species <- fct_reorder(d$Species, d$brain_index, min)

```

 Now we can get a rather "busy" but useful visualisation of the residuals. Notice the use of geom_col to produce columns. I.e. we get a barplot.  The bars can be coloured by diet for added effect. The default colours are not very nice and can be changed quite easily, but I will leave it there for now.
 
```{r}
g1<-ggplot(d, aes(x=Species,y=brain_index, label=Species, fill=Diet)) + geom_col() +  coord_flip()
g1
```

## Further work

This analysis shows how regression can be used to "control for" the effect of one variable. The residuals from the model can be thought of informally as a braininess index. The higher the value the larger the brain size is after taking into account body weight. So you now have a derived variable that can be related to diet or sleeping patterns if you want to do so using the methods you have already seen.

