---
title: 'Homework #4: Boosting, Dimension reduction, Clustering'
output: html_document
---

```{r setup, echo=FALSE}
set.seed(7)
```

1. Load the iSPY1 dataset by running the following. Notice that there is a `fakeSite` column with values ranging from 1-8. We've added this column to simulate combining data from 8 different sites.
```{r data}
ispy_dat <- read.csv("ispy1doctored_site.csv")
```

2. (1 point) Hold out sites 7 and 8 for testing. Store the training data in `datTrain` and the test data in `datTest`.
```{r}
datTrain <- ispy_dat[ispy_dat$fakeSite < 7, ]
datTest <- ispy_dat[ispy_dat$fakeSite >= 7, ]
```

3. (2 points) In this homework, we will manually implement site-wise 3-fold cross-validation (i.e. two sites per fold) rather than using the `caret` package. Create binary masks to split the rows in `datTrain` into 3 folds, where sites 1 and 2 are in the first fold, sites 3 and 4 are in the second fold, and sites 5 and 6 are in the third fold. Create a list with these three binary masks. It may be helpful to reference this list later in this homework.
```{r}
fold1 <- datTrain$fakeSite %in% c(1,2)
fold2 <- datTrain$fakeSite %in% c(3,4)
fold3 <- datTrain$fakeSite %in% c(5,6)
foldList <- list(fold1, fold2, fold3)
```


## Boosting

4. Load the `gbm` package.
```{r}
library(gbm)
```

5. (3 points) We will select the hyperparameters for a gradient boosted model using site-wise three-fold CV. Create a function named `fit_fold` that takes as input the fold number `fold_idx`, number of trees, and interaction depth. The function will fit a gradient boosted model using `gbm` that predicts `MRI_LD_Tfinal` using all the predictors except for `fakeSite`. Train on all the folds except for the `fold_idx`-th one. The function should output the mean squared error of the fitted model on the held out fold.  Use the folds you made in question 3. Fix the shrinkage hyperparameter in `gbm` as 0.01.
```{r}
# Using either Gaussian or poisson is accepted for homework. Using a poisson distribution can be preferable because it only outputs positive predictions. Still, the MSEs from using either method are quite similar.
fit_fold <- function(fold_idx, n.trees=50, interaction.depth=1) {
  subTrainDat <- datTrain[!foldList[[fold_idx]],]
  subValDat <- datTrain[foldList[[fold_idx]],]
  gbmTemp <- gbm(MRI_LD_Tfinal ~ . - fakeSite, data=subTrainDat, distribution="poisson", n.trees=n.trees, interaction.depth=interaction.depth, shrinkage = 0.01)
  gbmPred <- predict(gbmTemp, newdata=subValDat, type="response", n.trees=n.trees)
  mean((gbmPred - subValDat$MRI_LD_Tfinal)^2)
}
```

6. (4 points) Using the function you made in question 5, tune the number of trees and interaction depth using site-wise three-fold CV. Search over the values `n.trees=100, 200, 400, 800, 1600` and `interaction.depth=1, 2`. Which hyperparameter values minimize the cross-validated mean squared error?
```{r}
# Function takes in hyperparameters and outputs the CV error
do_cv <- function(hyperparams) {
  cv_errs <- sapply(seq(3), fit_fold, n.trees=hyperparams["n.trees"], interaction.depth=hyperparams["interaction.depth"])  
  mean(cv_errs)
}

# Create grid for tuning
tuneGrid <- expand.grid(
  n.trees=c(100, 200, 400, 800, 1600),
  interaction.depth=c(1,2))

# Do CV
cvRes <- apply(tuneGrid, 1, do_cv)
tuneGrid["cv_err"] = cvRes

bestHyperparam <- tuneGrid[which.min(cvRes),]
bestTrees <- bestHyperparam$n.trees
bestInteraction <- bestHyperparam$interaction.depth
print(paste("Best n.trees:", bestTrees))
print(paste("Best interaction.depth:", bestInteraction))

#note with set.seed(7) the best interaction.depth is 1. However, with various seeds, sometimes it is 2. Best n.trees is always 400.
```

7. (2 points) Plot the cross-validated error with respect to `n.trees`. Keep `interaction.depth` fixed at 1.

```{r}
tuneGrid$cv_err <- cvRes
subTuneGrid <- tuneGrid[tuneGrid$interaction.depth == 1,]
plot(subTuneGrid$n.trees, subTuneGrid$cv_err, type = "b")
```

8. (1 point) Refit the gradient boosted model on all the training data (`datTrain`) using the hyperparameters that minimized the CV error.
```{r}
gbm_final <- gbm(MRI_LD_Tfinal ~ . - fakeSite, data=datTrain, distribution="gaussian", n.trees=bestTrees, interaction.depth=bestInteraction, shrinkage = 0.01)
```

9. (1 point) Evaluate the MSE of the fitted model on the test data.
```{r}
predTestGbm <- predict(gbm_final, newdata = datTest, n.trees=bestTrees)
mean((predTestGbm - datTest$MRI_LD_Tfinal)^2)
```


## Kmeans

10. (1 point) Let's perform kmeans on the iSPY data. Remove the columns "race", "HR_HER2status", and "fakeSite". Call this new dataset `ispy_subdat`.
```{r}
library(dplyr)
ispy_subdat <- ispy_dat %>% select(-c("race", "HR_HER2status", "fakeSite"))
```

11. (3 points) Tune the number of clusters used in K-means. To do this, use the function `fviz_nbclust` from the `factoextra` library. The function `fviz_nbclust` determines and visualizes the optimal number of clusters using different methods (within cluster sums of squares, average silhouette and gap statistics). Plot the average silhouette with respect to the number of clusters by passing in the argument `method="silhouette"`. What is the optimal number of clusters according to the silhouette statistic?
```{r}
library(factoextra)
fviz_nbclust(ispy_subdat,
             FUNcluster = kmeans,
             method = "silhouette") 
# Optimal number of clusters is 2 according to the silhouette statistic.
```

12. (3 points) Refit k-means using the optimal number of clusters with 15 random initializations. Use the function `fviz_cluster()` to plot the clusters from K-means. Observations are represented by points in the plot, using principal components if $p > 2$. An ellipse is drawn around each cluster.

```{r}
tuned_km <- kmeans(ispy_subdat, centers=2, nstart=15)
fviz_cluster(tuned_km, data=ispy_subdat)
```