```{r setup, include=FALSE}
library("knitr")
# library("caret")       # for general model fitting
library("rpart")       # for fitting decision trees
library("ipred")  
library("foreach")     # for parallel processing with for loops
library("ggplot2")
library("rpart.plot")
library("doParallel")  # for parallel backend to foreach
library("data.table")
library("tidymodels")
library("randomForest")

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

## Reading the data

For this practical session on logistic regression we are using a dataset on the relationship between cleft lip in dogs (Nova Scotia Duck Tolling Retriever, NSDTR) and SNP genotypes ([data](https://datadryad.org/stash/dataset/doi:10.5061/dryad.j8r8q); [paper](https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1005059)).

The public dataset downloaded from Dryad is a *Plink* `.tped/.tfam` file. The data have been preprocessed:

-   filtered (SNP quality, call-rate, MAF)
-   imputed (imputation of missing genotypes using LHCI: localised haplotype-clustering imputation)
-   selected (only SNPs on chromosomes 25, 26, 27, 28, 29)

```{r}
dogs <- fread("../data/dogs_imputed.raw")
dogs <- dogs %>%
  select(-c(IID,FID,PAT,MAT,SEX))

dogs$PHENOTYPE = dogs$PHENOTYPE-1
```

## Trees and bagging

#### Classification tree

Let's start with a classification tree. We split the data in training and testing sets:

```{r subset, echo=FALSE}
seed = 137
set.seed(seed)

##Training data
n = nrow(dogs) ## sample size
n_training = round(0.75*n,0)
n_test = n - n_training
training_records <- sample(n,n_training)
x_train <- dogs[training_records,]
x_test <- dogs[-training_records,-1]
y_test <- as.factor(dogs[-training_records,]$PHENOTYPE)
print(y_test)
```

```{r}
table(x_train$PHENOTYPE)
```

Then fit the classification tree model to the training set:

-   function `rpart` (from the *rpart* package)
-   method = "class" (classification)
-   `minsplit`: minimum number of observations (records) in a split (splits with fewer than 3 records will not be considered/attempted)

```{r classification_tree}
##Model without bagging
no_bag_model <- rpart(PHENOTYPE ~ ., data = x_train, method = "class", control = rpart.control(minsplit = 3))

rpart.plot(no_bag_model)
```

In the *tree-plot* we see:

-   the **classification**
-   the modelled class probability $\rightarrow$ $P(y=1|x)$ (i.e. the probability of the class conditioned on the node)
-   the percentage of observations used at that node

With the trained model we can make predictions in the test set:

```{r}
predictions <- predict(no_bag_model, x_test, type = "class")
table(predictions, y_test)
```

#### Bagging

We now take 10 **bootstrapped samples** of the data, and fit a classification tree to each sample. We use OOB observations to measure the prediction error. At each bagging replicate (bootstrapped sample) typically $1/3$ of the observations are not used for training (out-of-bag -OOB- observations), and can naturally be used to measure the prediction error:

1.  Fitting (multiple) tree models:

```{r bagging}
## fitting tree models to bootstrapped data (bagging)
n_model = 10
bagged_models <- list()
oob <- list()
for (i in 1:n_model) {
  
  records = seq(1,nrow(dogs))
  new_sample <- sample(records, size=nrow(dogs), replace=TRUE) ## bootstrapping here!!
  oob <- c(oob, list(records[!(records %in% new_sample)]))
  bagged_models <- c(
    bagged_models,
    list(rpart(PHENOTYPE ~ ., dogs[new_sample,], control = rpart.control(minsplit=3))))
}
```

2.  Making predictions from fitted models:

```{r}
##Getting estimate from the bagged model
oob_results = data.frame("rep"=NULL, "error"=NULL)
i = 1
for (bag_model in bagged_models) {

  preds = predict(bag_model, dogs[oob[[i]],-1]) ## predict by applying the trained model on the OBB data
  obs = dogs$PHENOTYPE[oob[[i]]]
  oob_results = rbind.data.frame(oob_results, data.frame("rep"=i, "error"=1-sum(preds == obs)/length(preds))) ## check obs == preds
  i = i+1
}
```

```{r}
print(oob_results)
```

```{r}
print(paste("average OOB error:" , mean(oob_results$error)))
```

The average OOB error is `r mean(oob_results$error)`.

Alternatively, we can proceed with splitting the data in **training and testing sets**, fit the classification trees to bootstrapped replicates of the training data, and then estimate the error rate on the test data.

Also here we will have OOB observations, but the error rate in this case would tend to be overestimated since the training set is smaller than the original dataset.

```{r bagging2}
n_model = 10
bagged_models <- list()
oob <- list()
for (i in 1:n_model) {
  
  new_sample <- sample(training_records, size=length(training_records), replace=TRUE)
  oob <- c(oob, list(training_records[!(training_records %in% new_sample)]))
  bagged_models <- c(
    bagged_models,
    list(rpart(PHENOTYPE ~ ., dogs[new_sample,], control = rpart.control(minsplit=3))))
}
```

```{r}
##Getting estimate from the bagged model
bagged_result = NULL
i = 0
for (bag_model in bagged_models) {
 
  if (is.null(bagged_result)) {
    bagged_result = predict(bag_model, x_test)
  } else {
   bagged_result = (i*bagged_result+predict(bag_model, x_test))/(i+1) ## % of times each observation was classified as 1
  }
i = i+1
}

diff = abs(bagged_result - (as.integer(y_test)-1)) ## absolute differences between predictions and observations (averaged over the n replicates of the model)
error = sum(n_model*diff)/(length(y_test)*n_model) ## average error over number of predictions (size of test set time n. of replicates)
print(paste("The bagging error is", error))
```

```{r}
df <- data.frame("pred"=bagged_result,"obs"=as.integer(y_test)-1)
head(df)
```

Total average test error rate is `r error`. The breakdown in average FPR and FNR over the bagging replicates is shown below:

```{r}
df %>%
  group_by(obs) %>%
  summarise(N=n(),err=sum(10*abs(pred-obs))/(N*10))
```

Final bagging classification.

```{r}
preds = ifelse(bagged_result > 0.5,1,0)
table(preds,y_test)
```

## Random Forest

We can now make the next step from bagging (averaging over bootstrapped trees) to Random Forest (adding randomly sampled subsets of features to bootstrapped data):

-   first, split the data in training and test sets (we keep y and x separate here)
-   then we fit a RF model to the training set, without tuning any hyperparameters (for now)
-   the function `randomForest` internally evaluates the model on OOB data

```{r rf, echo=FALSE}
##Training data
set.seed(123)
n = nrow(dogs) ## sample size
n_training = round(0.8*n,0)
n_test = n - n_training
training_records <- sample(n,n_training)
x_train <- as.matrix(dogs[training_records,-1])
y_train <- as.factor(dogs[training_records,]$PHENOTYPE)
x_test <- as.matrix(dogs[-training_records,-1])
y_test <- as.factor(dogs[-training_records,]$PHENOTYPE)
print(y_test)
```

```{r}
rf.dogs <- randomForest(x = x_train, y = y_train, ntree = 10)
print(rf.dogs)
```

#### Let's tune it up a little!

Now we use a convenient built-in function to tune hyperparameters in random forest: the `tuneRF` function, from the *randomForest* package:

-   we tune just the number of variables to sample in each tree, `mtry`
-   mtryStart: starting value of `mtry` (usually: $\sqrt{p}$ for classification; $p/3$ for regression)
-   stepFactor: at each iteration, `mtry` is inflated (or deflated) by this value
-   improve: the (relative) improvement in OOB error must be by this much for the search to continue

```{r}
tuned_rf <- tuneRF(x_train, y_train,             ## DATA
                   mtryStart = 25, ntreeTry=50,  ## HYPERPARAMETERS
                   stepFactor=2, improve=0.01,
                   trace=TRUE, plot=TRUE) ## we take a smaller starting mtry to better see the plot

head(tuned_rf)
```

```{r}
num_row <- which(tuned_rf[,2] == min(tuned_rf[,2])) ## row corresponding to min OOB error
num_vars <- tuned_rf[num_row,1]  ## mtry value corresponding to min OOB error
print(num_vars)
```

Now we fit a final RF model with the tuned value for `mtry` (number of variables subsampled in each tree), and build a small forest with 100 trees:

```{r}
rf.dogs <- randomForest(x = x_train, y = y_train, ntree = 100, mtry = num_vars, importance = TRUE)
rf.dogs
```

#### Variable importance

From RF models, we can extract variable importance (here we take the top 20 variables ordered for importance):

```{r}
imp <- as.data.frame(importance(rf.dogs))
summary(imp)
top_vars <- imp %>%
  arrange(desc(MeanDecreaseGini)) %>%
  top_n(20) %>%
  select(MeanDecreaseGini) %>%
  mutate(name_var = gsub("\\_.*$","",row.names(.)))

top_vars$name_var <- factor(top_vars$name_var, levels = rev(top_vars$name_var))
```

```{r varimp}
p <- ggplot(top_vars, aes(x=name_var,y=MeanDecreaseGini)) + geom_bar(stat="identity")
p <- p + theme(axis.text.x = element_text(angle = 90, hjust = 1))
p <- p + coord_flip()
p
```

#### Predictions

We are now ready to evaluate our tuned RF model on the **test data**:

```{r}
probs <- predict(rf.dogs, newdata = x_test, type = "prob")
predictions <- ifelse(probs[,"1"] > 0.5, 1, 0 )
table(predictions, y_test)
```

```{r label="roc"}
library("ROCit")

load("ROClogit.RData")
load("ROClasso.RData")

ROCit_rf <- rocit(score=probs[,"1"],class=y_test)
plot(ROCit_rf, col = c(1,"gray50"), 
     legend = FALSE, YIndex = FALSE)
lines(ROCit_lasso$TPR~ROCit_lasso$FPR,col = 2, lwd = 2)
lines(ROCit_logit$TPR~ROCit_logit$FPR,col = 3, lwd = 3)
legend("bottomright", col = c(1,2,3),
       c("RF","Lasso", "Logistic regression"), lwd = 2)
```
