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

------------------------------------------------------------------------



```{r}
library(MASS)
library(ggplot2)
library(dplyr)
library(gridExtra)
```

# Animals

In this exercise you will work with the `Animals` data from the package `MASS`. 
These data include two variables measuring the brain and body weight of various 
animals. The central question in this exercise is: "Can we predict the weight of 
an animal's brain from the weight of the animal's body? The exercise shows that
with some simple data transformations and adding a predictor we can go from 0\% 
variance explained to almost 100\%!

## Linear model

a. Check the help file for `Animals`, and display the row names to get an idea 
of what kinds of animals are in the data.

```{r}
rownames(Animals)
```

b. Let's start with with some data visualizations. Display the boxplots showing 
the distributions of `brain` and `body`. 

```{r fig.asp = 0.4}
grid.arrange(
  ggplot(Animals) +
    geom_boxplot(aes(y = brain)),
  ggplot(Animals) +
    geom_boxplot(aes(y = body)), 
  nrow = 1)
```

c. Display a scatter plot with `body` on the x-axis and `brain` on the y-axis, 
and add a linear regression line including the confidence band. Interpret the 
plot.

```{r fig.width = 5, fig.height = 3}
ggplot(Animals, aes(x = body, y = brain)) +
  geom_point() + 
  geom_smooth(method = "lm") +
  theme_minimal()
```

d. Let's see how the linear regression model performs. Fit the linear model 
$m_0$ predicting `brain` from `body`; display its summary, and interpret the 
parameter estimates and the R-squared. 

```{r}
# type your code here
```

e. Display the diagnostic residual plots for the previous model, and interpret 
the results.

```{r fig.asp = 1.2}
par(mfcol = c(2, 2)) # this prepares a plotting area for 4 plots in a 2x2 display
# type your code here to make the plot of the lm class object:

```

## Data transformations

Since the variables `brain` and `body` are substantially skewed, a data 
transformation may improve the fit of the model.

a. Apply transformations to normalize the distributions of the variables `body` 
and `brain`, and check the result with boxplots. 

```{r fig.asp = 0.4}
grid.arrange(
  ggplot(Animals) +
    geom_boxplot(aes(y = log(brain))),
  ggplot(Animals) +
     geom_boxplot(aes(y = log(body))),
    nrow = 1)
```

b. Display a scatter plot of the transformed versions of `brain` and `body`, and 
add a linear regression line. 

```{r fig.width = 5, fig.height = 3}
ggplot(Animals, aes(x = log(body), y = log(brain))) +
  geom_point() + 
  geom_smooth(method = "lm") +
  theme_minimal()
```

d. Fit the linear model `m1` with the transformed `body` predicting the 
transformed `brain`; display its summary, and interpret the parameter estimates 
and the R-squared. 

```{r}
# Type your code here
```

e. Display the diagnostic plots for `m1` and interpret. Are the model 
assumptions fulfilled?

```{r fig.asp = 1.2}
# Type your code here
```

## Adding a predictor

Additional predictors can help to satisfy the assumptions. For the data at hand, 
the residuals plots show three outliers that have one thing at common; they are 
prehistoric species. Let's see what happens if we add a predictor to the model 
that distinguishes between extinct and other species.   

a. Make the factor `species` with levels "extinct" for the prehistoric animals 
and "other" for the others, and add this factor to the data frame `Animals`.

```{r}
Animals <- mutate(Animals, 
                  species = recode_factor(rownames(Animals),
                                          Brachiosaurus = "extinct",
                                          Triceratops   = "extinct",
                                          Dipliodocus   = "extinct", 
                                          .default      = "other")
                  )
```

**OR**

```{r eval = FALSE}
Animals <- mutate(Animals, 
                  species = case_when(
                    rownames(Animals) %in% c("Brachiosaurus",
                                             "Triceratops",
                                             "Dipliodocus") ~ "extinct",
                    TRUE ~ "other")
                  )
```

b. Fit the model `m2` by adding the predictor `species` to model `m1`, display 
its summary, and interpret the parameter estimates and the R-squared.

```{r}
# Type your code here
```

c. Display and interpret the diagnostic plots for `m2`.

```{r fig.asp = 1.2}
# Type your code here
```

The diagnostic plots show outliers that have a higher (logarithm of) `brain` 
than predicted by the model. These outliers belong to the class of species know 
as primates. This result suggests that the predictor `species` could be extended 
with the category `primates`. 

d. Add the level "primate" to the variable `species` (check which of the animals 
are primates). 

```{r}
Animals  <- mutate(Animals, species = recode_factor(rownames(Animals),
                                  'Brachiosaurus' = "extinct",
                                  'Triceratops'   = "extinct",
                                  'Dipliodocus'   = "extinct", 
                                  'Human'         = "primate",
                                  'Chimpanzee'    = "primate",
                                  'Gorilla'       = "primate",
                                  'Rhesus monkey' = "primate",
                                  'Potar monkey'  = "primate",
                                  .default        = "other")
                   )
```


**OR**

```{r eval = FALSE}
Animals <- mutate(Animals, 
                  species = case_when(rownames(Animals) %in% c('Human', 
                                                               'Chimpanzee', 
                                                               'Gorilla',
                                                               'Rhesus monkey', 
                                                               'Potar monkey') ~ "primate", 
                                      TRUE ~ species)
                  )
```

e. Fit the model `m3` with the modified predictor `species`; display its 
summary, and interpret the parameter estimates and the R-squared.

```{r}
# Type your code here
```

f. Display and interpret the diagnostic plots for `m3`. What do the residual 
plots say about the brain weight of the Gorilla?

```{r fig.asp = 1.2}
# Type your code here
```

g. Test the fit of the models `m1` to `m3` with the `anova()` function.

```{r}
# Type your code here
```


------------------------------------------------------------------------

End of practical
