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

```{r, message=FALSE, warning=FALSE}
library(tidyverse)
library(corrplot)
```

```{r}
myData <- read.table(file="../data/bats.csv", sep=";", skip=3, header=TRUE)
names(myData)
```

## Study of the relationship between brain weight and body mass

```{r message=FALSE, out.width="80%", fig.align = "center"}
phyto <- myData[(myData$Diet == 1),]
ggplot(phyto, aes(x=BOW, y=BRW)) + geom_point() +
  xlab("Body mass (g)") + ylab("Brain weight (µg)") +
  ggtitle("Brain weight with respect to body mass in phytophagous") +
  theme(plot.title = element_text(size=12))
```

It looks like body mass (in grams) has a clear influence on brain weight (in micrograms). This relation can be analysed fitting a linear model of brain weight as a function of body mass. The fitted line is also shown in next plot.

```{r message=FALSE, out.width="80%", fig.align = "center"}
reg1 <- lm(BRW ~ BOW, data=phyto)
ggplot(phyto, aes(x=BOW, y=BRW)) + geom_point() + geom_smooth(method="lm") +
  xlab("Body mass (g)") + ylab("Brain weight (µg)") +
  ggtitle("Brain weight with respect to body mass in phytophagous (with regression line)") +
  theme(plot.title = element_text(size=12))
```

With this expression, R estimates the $\beta$ coefficients of the formula $Y = \beta_0 + \beta_1X + \epsilon$, where $X$ is the vector containing the values of the measured body masses and $Y$ contains the measured brain weights. $\epsilon$ is random noise.

```{r}
summary(reg1)
```

In this case, a few results can be observed:

- The intercept $\beta_1$ is estimated as 623.45.
- The p-value for the employed test statistics is extremely small, lower than 2.2e-16.
- The null hypothesis of this test is $H_0 : \beta_1 = \beta_2 = 0$. This hypothesis is therefore rejected.
- Body mass has a significative relation with brain weight. The intercept is also important.
- The coefficient of determination $R^2$ of the fitted model is 0.95, so a big portion of the global variation is explained by the model.

Next, analysis of variance can be performed.

```{r}
anova(reg1)
```

This confirms that body weight is highly significative. In addition, the previous table contains information about the model residual, with the sum of residual squares being 4253838. The residuals have been computed during model fitting and can be plotted against the predicted values.

```{r out.width="80%", fig.align = "center"}
tmp <- data.frame(reg1$fitted.values, reg1$residuals)
ggplot(tmp, aes(x=reg1.fitted.values, y=reg1.residuals)) + geom_point() +
  xlab("Predicted") + ylab("Residuals") + ggtitle("Diagnostic graph") +
  theme(plot.title = element_text(size=12))
```

Clearly, something is not going well: most prediction values are below 3750 µg, except one with predicted brain weight around 10000 µg. As this last point is so far from the others, it influences the model results worsening the prediction for values in the middle. For now, we can consider it as an outlier.

Cook's distance can be used for outlier detection in the fitted model:

```{r out.width="80%", fig.align = "center"}
plot(reg1, 4)
```

Observation number 7 has distance much higher than the higher and probably corresponds to the previously identified point.

```{r}
myData[7,]
```

Indeed, `r myData[7,]$Species` has body mass `r myData[7,]$BOW` g and brain weight `r myData[7,]$BRW` µg.

We can therefore redo the previous analysis without taking it into account.
```{r message=FALSE, out.width="80%", fig.align = "center"}
phytobis <- phyto[which(phyto$BRW<8000),]
reg2 <- lm(BRW ~ BOW, data=phytobis)
ggplot(phytobis, aes(x=BOW, y=BRW)) + geom_point() + geom_smooth(method="lm") +
  xlab("Body mass (g)") + ylab("Brain weight (µg)") +
  ggtitle("Brain weight with respect to body mass in phytophagous (with regression line)") +
  theme(plot.title = element_text(size=12))
summary(reg2)
```

Clearly, now the regression line better fits intermediate point and is not affected by the removed extreme observation. The predicted coefficient for the body weight predictor is higher than before, as well as the coefficient of determination.

A few graphs can therefore plotted for model diagnosis, including the residuals agains predicted shown before. For the first model, including the "outlier":

```{r fig.height=7}
par(mfcol=c(2,2))
plot(reg1)
```
If the "outlier" is removed:
```{r fig.height=7}
par(mfcol=c(2,2))
plot(reg2)
```
Important differences can be noticed:

- The plots related to residuals vs fitted and scale-location are more "flat" for the second model, so it provides more accurate predictions than the first one.
- According to the Q-Q plot, it is more probable that the data used to fit the second model are drawn from a normal distribution, with respect to the first model.

Therefore, the second model better satisfies common assumptions of linear regression.


## Study of the contribution to the total weight of each part of the brain

Preivous analysis can be expanded, using other variables (related to brain parts) to explain brain weight. The possible explanatory variables are AUD (auditory nuclei volume), MOB (main olfactory bulb volume) and HIP (hippocampus volume).

```{r out.width="80%", fig.align = "center"}
phytoNum <- phyto[, c(4:8)]
mat.cor <- cor(phytoNum)
corrplot(mat.cor, type="upper")
```

We see high positive correlation in the couples, for example, (BOW, BRW) and (MOB, HIP), while the AUD variable is not particularly correlated with any other. Let's statistically analyse the correlation between some variables, using the Pearson test.

```{r}
cor.test(phyto$BRW, phyto$HIP)
cor.test(phyto$BRW, phyto$MOB)
cor.test(phyto$BRW, phyto$AUD)
```

Resulting p-values are very low if BRW is tested against HIP and MOB, while it is bigger (0.003) against AUD. What is shown in the previous figure is confirmed: correlation is important with HIP and MOB, less with AUD.

```{r}
regm <- lm(BRW ~ AUD + MOB + HIP, data=phytobis)
summary(regm)
anova(regm)
```

The estimated model predicts the $\beta$ coefficient of the regression model $Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_3$, where $Y$ is the BRW variable and $X_1$, $X_2$ and $X_3$ are AUD, MOB and HIP respectively. The objective is therefore to predict the brain mass of a bat from the volumes of its auditory nuclei, main olfactory bulb and hyppocampus. Theoretically, they may all have an impact on the brain mass, even though the main olfactory bulb can be present in the bat's nose and not in its brain.

The coefficient associated to the three variables are, respectively, 47.989, -2.444 and 15.981. The coefficient related to MOB is not significative, showing a high p-value, and could therefore be removed from the analysis.

Additional analysis can be conducted on this model. Here, the "batman" is taken into account.

```{r}
reg0 <- lm(BRW ~ 1, data=phyto)
step(reg0, scope=BRW ~ AUD + MOB + HIP, direction="forward")
```

According to the `help`, `step` performs a stepwise search for a best model according to the AIC criterion. In this case, intercept only is considered at the first step and all explanatory variables are included at the end. In sequence, the models `BRW ~ 1`, `BRW ~ HIP`, `BRW ~ HIP + MOB` and `BRW ~ HIP + MOB + AUD` are evaluated, with the latter having the lowest information loss.

## Link between volume of the auditory part and diet

```{r fig.show="hold", out.width="50%"}
myData$Diet_F <- as.factor(myData$Diet)
with(myData, plot(AUD~Diet))
with(myData, plot(AUD~Diet_F))
```

The two graphs basically show the same information, but in different ways. The first displays the single data points and it is clear that few examples are available for diets 2 and 4. The second displays the boxplots instead and shows data distributions for each different diet. Boxplots for diets 2 and 4 are condensed around the mean, while the ones with more data points are more spread.

```{r}
lm <- lm(AUD~Diet_F, data=myData)
anova(lm)
```

The p-value of the ANOVA test is quite high, 0.432, so it is not possible to conlude that there is a significant difference between the different diets. The bats' diets don't have therefore an impact on the volume of their auditory part of the brain.
