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.
ispy_dat <- read.csv("ispy1doctored_site.csv")
  1. (1 point) Hold out sites 7 and 8 for testing. Store the training data in datTrain and the test data in datTest.
datTrain <- ispy_dat[ispy_dat$fakeSite < 7, ]
datTest <- ispy_dat[ispy_dat$fakeSite >= 7, ]
  1. (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.
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

  1. Load the gbm package.
library(gbm)
## Loaded gbm 2.1.5
  1. (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.
# 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)
}
  1. (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?
# 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))
## [1] "Best n.trees: 400"
print(paste("Best interaction.depth:", bestInteraction))
## [1] "Best interaction.depth: 1"
#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.
  1. (2 points) Plot the cross-validated error with respect to n.trees. Keep interaction.depth fixed at 1.
tuneGrid$cv_err <- cvRes
subTuneGrid <- tuneGrid[tuneGrid$interaction.depth == 1,]
plot(subTuneGrid$n.trees, subTuneGrid$cv_err, type = "b")

  1. (1 point) Refit the gradient boosted model on all the training data (datTrain) using the hyperparameters that minimized the CV error.
gbm_final <- gbm(MRI_LD_Tfinal ~ . - fakeSite, data=datTrain, distribution="gaussian", n.trees=bestTrees, interaction.depth=bestInteraction, shrinkage = 0.01)
  1. (1 point) Evaluate the MSE of the fitted model on the test data.
predTestGbm <- predict(gbm_final, newdata = datTest, n.trees=bestTrees)
mean((predTestGbm - datTest$MRI_LD_Tfinal)^2)
## [1] 419.1062

Kmeans

  1. (1 point) Let’s perform kmeans on the iSPY data. Remove the columns “race”, “HR_HER2status”, and “fakeSite”. Call this new dataset ispy_subdat.
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
ispy_subdat <- ispy_dat %>% select(-c("race", "HR_HER2status", "fakeSite"))
  1. (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?
library(factoextra)
## Loading required package: ggplot2
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
fviz_nbclust(ispy_subdat,
             FUNcluster = kmeans,
             method = "silhouette") 

# Optimal number of clusters is 2 according to the silhouette statistic.
  1. (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.
tuned_km <- kmeans(ispy_subdat, centers=2, nstart=15)
fviz_cluster(tuned_km, data=ispy_subdat)