3.training_testing.Rmd

---
title: "training and test sets"
author: "Filippo Biscarini"
date: "7/6/2021"
output: html_document
---

```{r setup, include=FALSE}
library("broom")
library("knitr")
library("ggplot2")
library("corrplot")
library("reshape2")
library("tidyverse")
library("tidymodels")
library("data.table")

knitr::opts_chunk$set(echo = TRUE)
```

## Reading the data

For this practical session on linear regression we are using a dataset on the relationship between the age of wild bats and DNA methylation at specific CpG sites ([data](https://datadryad.org/stash/dataset/doi:10.5061/dryad.rn0198d); [paper](https://onlinelibrary.wiley.com/doi/abs/10.1111/1755-0998.12925)).

The public dataset downloaded from Dryad is an `.xlsm` file, and can be read into *R* using the `read.xlsx()` function from the **xlsx** package.

```{r label="reading_data"}
ch4 <- readxl::read_excel("../data/DNA methylation data.xlsm")
ch4 <- na.omit(ch4[,-c(1,3)])
head(ch4)
```

## Training and testing sets

We now need to subset the data into a **training set** and a **testing set** We can randomly assign records from the original dataset to the two subsets, for instance $80\%$ to the training set and $20\%$ to the test set:

```{r pressure, echo=FALSE}
# set.seed(295)
n = nrow(ch4) ## sample size
n_training = round(0.8*n,0)
n_test = n - n_training

training_records <- sample(n,n_training)
training_set <- ch4[training_records,]
test_set <- ch4[-training_records,]
```

```{r}
## sanity checks
nrow(training_set)
nrow(test_set)
```

## Fitting the multiple linear regression model

Note: we use only the **training data** to fit the model!

```{r}
fit <- lm(Age ~ ., data = training_set)
tidy(fit)
r_squared <- summary(fit)$r.squared
print(paste("
cbind.data.framR^2", round(r_squared,4)))
```

The coefficient of determination $R^2$ is `r r_squared`

```{r}
data.frame("coefficient"=coef(fit))
```

[Question: How do we interpret the model coefficients?]{style="color:red"}

### The test error

```{r}
predictions <- predict(fit, newdata = test_set[,-1], interval="none", type = "response", na.action=na.pass)

cbind.data.frame("test_age"=test_set[,1], "predictions"=round(predictions,2))
```

```{r}
r_pearson = cor(test_set$Age,predictions)
r_spearman = cor(test_set$Age,predictions, method = "spearman")
mse = mean((test_set$Age-predictions)^2)
rmse = sqrt(mse)
nrmse = sqrt(mse)/mean(test_set$Age)

res.test = data.frame("metric"=c("MSE","RMSE","NRMSE","r pearson","r spearman"),
           "value"=c(mse,rmse,nrmse,r_pearson,r_spearman))
res.test
```

```{r}
ggplot(data.frame("test_age"=test_set$Age,"predictions"=predictions), aes(x=predictions,y=test_age)) + geom_point()
```

### The training error

We can compare the test error measured above with the training error, i.e. the model performance measured on the training data:

```{r}
predictions <- predict(fit, newdata = training_set[,-1], interval="none", type = "response", na.action=na.pass)
```

```{r}
r_pearson = cor(training_set$Age,predictions)
r_spearman = cor(training_set$Age,predictions, method = "spearman")
mse = mean((training_set$Age-predictions)^2)
rmse = sqrt(mse)
nrmse = sqrt(mse)/mean(training_set$Age)

res.train = data.frame("metric"=c("MSE","RMSE","NRMSE","r pearson","r spearman"),
           "value"=c(mse,rmse,nrmse,r_pearson,r_spearman))

res.train
```

```{r}
ggplot(data.frame("training_y"=training_set$Age,"predictions"=predictions), aes(x=predictions,y=training_y)) + geom_point()
```

```{r}
res.train |> inner_join(res.test, by = "metric") |> rename(training = value.x, test = value.y) |>
  mutate(tst_trn_pct = round(100*(test-training)/training, 3))
```

## Exercise 3.1

-   try a different training/testing split
-   train the model
-   measure the performance of the model

#### 1. split the data

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

#### 2. train the model

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

#### 3. evaluate the model on the test data

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

## Validation set approach

We now implement the validation set approach. First, we choose the relative proportions for the training and test sets, e.g. $70\%$ and $30\%$. We feed the original dataset and the chose proportions in a function to implement the validation set approach for our regression problem, and then we repeat this 20 times, to get an estimate of the variability of the predictor $\hat{f}(x)$.

```{r validation_set, echo=FALSE}
validation_approach <- function(dataset,proportion) {
  
  ## prepare the subsets
  n = nrow(dataset) ## sample size
  n_training = round(proportion*n,0)
  n_test = n - n_training

  training_records <- sample(n,n_training)
  training_set <- ch4[training_records,]
  test_set <- ch4[-training_records,]
  
  ##fit the model
  fit <- lm(Age ~ ., data = training_set)
  
  ## evaluate the model
  predictions <- predict(fit, newdata = test_set[,-1], interval="none", type = "response", na.action=na.pass)
  r_pearson = cor(test_set$Age,predictions)
  r_spearman = cor(test_set$Age,predictions, method = "spearman")
  mse = mean((test_set$Age-predictions)^2)
  rmse = sqrt(mse)
  nrmse = sqrt(mse)/mean(test_set$Age)

  # res <- data.frame("metric"=c("MSE","RMSE","NRMSE","r pearson","r spearman"),
  #          "value"=c(mse,rmse,nrmse,r_pearson,r_spearman))
  res <- data.frame(
           "value"=c(mse,rmse,nrmse,r_pearson,r_spearman))
  
  return(res)
}

n_replicates = 20
results <- replicate(n=n_replicates, validation_approach(ch4,0.7),simplify = FALSE)

df <- do.call(cbind.data.frame, results)
df <- transpose(df)
names(df) <- c("MSE","RMSE","NRMSE","r pearson","r spearman")
head(df,10)
```

We can now summarise and visualise the results:

```{r}
mDF <- reshape2::melt(df, variable.name = "metric")
mDF %>%
  group_by(metric) %>%
  summarise(N=n(), avg=mean(value), std=sd(value))
```

```{r}
p <- ggplot(mDF, aes(x = factor(1), y = value))
p <- p + geom_jitter(aes(colour=metric)) 
p <- p + facet_wrap(~metric, scales = "free")
p <- p + xlab("")
p
```

## k-fold cross-validation

Now, we implement 5-fold cross-validation. We set k=5 and then assign each observation from the original dataset to a fold:

```{r}
n = nrow(ch4)
k = 5
folds <- sample(x = seq(1,k), size = n, replace = TRUE)
table(folds)
```

We now fit and evaluate the model *k* times, and save results:

```{r}
cv.error = data.frame("fold"=NULL,"MSE"=NULL,"RMSE"=NULL,
                   "NRMSE"=NULL,"r_pearson"=NULL,"r_spearman"=NULL)
for(i in 1:k) {
  
  print(paste("using validation fold",i))
  training_set <- ch4[folds!=i,]
  fit <- lm(Age ~ ., data = training_set)
  
  test_set <- ch4[folds==i,]
  predictions <- predict(fit, newdata = test_set[,-1], interval="none", type = "response", na.action=na.pass)
  
  r_pearson = cor(test_set$Age,predictions) 
  r_spearman = cor(test_set$Age,predictions, method = "spearman")
  mse = mean((test_set$Age-predictions)^2)
  rmse = sqrt(mse)
  nrmse = sqrt(mse)/mean(test_set$Age)
  
  res = data.frame("fold"=i,"MSE"=mse,"RMSE"=rmse,
                   "NRMSE"=nrmse,"r_pearson"=r_pearson,"r_spearman"=r_spearman)
  cv.error <- rbind.data.frame(cv.error,res)
}
```

```{r}
print(cv.error)
```

```{r}
mCV <- reshape2::melt(cv.error, id.vars = "fold", variable.name = "metric")
mCV %>%
  group_by(metric) %>%
  summarise(N=n(),avg = mean(value), std = sd(value))
```

```{r}
p <- ggplot(mCV, aes(x = factor(1), y = value))
p <- p + geom_jitter(aes(colour=metric)) 
p <- p + facet_wrap(~metric, scales = "free")
p <- p + xlab("")
p
```

## Exercise 3.2

-   implement cross-validation with a different value for k (e.g. k=4, k=10, etc.): what happens to the estimated accuracy?

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