Maximum margin classifiers

  1. Run the following code to create a simple dataset.
X1 <- c(3,1,3,1,2,4,4)
X2 <- c(4,2,3,4,1,3,1)
Y <- as.factor(c(rep("Red",4),rep("Blue",3)))
dat <- data.frame(X1,X2,Y)
print(dat)
##   X1 X2    Y
## 1  3  4  Red
## 2  1  2  Red
## 3  3  3  Red
## 4  1  4  Red
## 5  2  1 Blue
## 6  4  3 Blue
## 7  4  1 Blue
  1. Plot the data. Color-code the points by their outcome value Y.
library(ggplot2)
qplot(x = X1, y = X2, data = dat, colour = Y)

  1. (3 points) Fit a maximal margin classifier with a linear decision boundary to this dataset using the svm function from e1071. Use the option scale=FALSE so that the model does not do post-processing of the training data. Plot the fitted model.
library(e1071)
svm_simple <- svm(Y ~ ., data=dat, cost=1000, kernel="linear", scale = FALSE)
plot(svm_simple, data=dat, formula = X2 ~ X1)

  1. (1 point) Use the coef function to extract model coefficients from the maximum margin classifier.
coef(svm_simple)
## (Intercept)          X1          X2 
##   0.9994192  -1.9996832   1.9998416
  1. (2 points) Suppose we are going to add another observation to the dataset and refit the maximum margin classifier. Give a new datapoint that, when added to this dataset, will not change the decision boundary. Refit the model to double check. Call this new datapoint point1.
point1 <- data.frame(X1=5, X2 = 2, Y="Blue")
dat_nochange <- rbind(dat, point1)
svm_nochange <- svm(Y ~ ., data=dat_nochange, cost=1000, kernel="linear", scale = FALSE)
plot(svm_nochange, data=dat, formula = X2 ~ X1)

print(paste("Are coefficients the same?", all(coef(svm_nochange) == coef(svm_simple))))
## [1] "Are coefficients the same? TRUE"
  1. (2 points) Now give a new datapoint that, when added to dataset dat, WILL change the decision boundary. Call this new datapoint point2. Pick an example point where the new dataset remains linearly separable.
point2 <- data.frame(X1=1, X2 = 1, Y="Blue")
dat_change <- rbind(dat, point2)
svm_change <- svm(Y ~ ., data=dat_change, cost=1000, kernel="linear", scale = FALSE)
plot(svm_change, data=dat, formula = X2 ~ X1)

print(paste("Are coefficients the same?", all(coef(svm_change) == coef(svm_simple))))
## [1] "Are coefficients the same? FALSE"

Support vector machines

  1. Read in the dementia data “dementia2.csv” into a data frame called dementia_dat. This data is from the UCSF Memory and Aging Center. The goal was to predict the type of dementia based on patterns of brain loss as measured through structural MRI.
dementia_dat <- read.csv("dementia2.csv")
  1. (1 point) We will now predict the dementia diagnosis based on the available predictors. The diagnosis is given in the multiclass outcome of MacCohort_kr. To start, generate a table of the elements in the MacCohort_kr variable.
table(dementia_dat$MacCohort_kr)
## 
##         AD      bvFTD    CONTROL nfPPAunspc        PSP      svPPA 
##        151         68        315         38         44         44
  1. (1 point) Set the random seed to 7 and then split the data into 2 sets (440 train, and 220 test)
set.seed(7)
groups <- c(rep("train",440),rep("test",220))
groups <- sample(groups,length(groups))
datTrain <- dementia_dat[groups=="train", ]
datTest <- dementia_dat[groups=="test", ]
  1. (3 points) Fit a support vector machine to predict MacCohort_kr using predictors Left_MTG_middle_temporal_gyrus, Left_Amygdala and Left_AIns_anterior_insula, with a radial basis kernel. Use 3-fold CV to tune the value of the C and sigma hyperparameters. Use the caret package with method=svmRadial and tune over the values C = 1e-2,1e-1,1,10,100,1000,10000 and sigma=1e-4,1e-3,1e-2,1e-1, 1,10. (You may need the kernlab package to run this.)
library(caret)
## Loading required package: lattice
train_control <- trainControl(method="cv", number=3)
svm_grid <-  expand.grid(C = c(1e-2,1e-1,1,10,100,1000,10000), 
                        sigma = c(1e-4,1e-3,1e-2,1e-1, 1,10))
cv_svm <- train(
  MacCohort_kr ~ Left_MTG_middle_temporal_gyrus + Left_Amygdala + Left_AIns_anterior_insula,
  data = datTrain,
  trControl=train_control,
  tuneGrid=svm_grid,
  method="svmRadial",
)
  1. (1 points) Determine the predictions of the SVM fit on the test dataset. What is the test accuracy?
predvals <- predict(cv_svm, newdata=datTest)
table(predvals)
## predvals
##         AD      bvFTD    CONTROL nfPPAunspc        PSP      svPPA 
##         51         23        131          0          0         15
print(paste("Accuracy", mean(predvals==datTest$MacCohort_kr)))
## [1] "Accuracy 0.640909090909091"
  1. (3 points) What is the accuracy of a model that randomly guesses the class label simply based on the fraction of observations in each class? Is the SVM doing better than random guessing?
true_fractions <- table(datTest$MacCohort_kr)/length(datTest$MacCohort_kr)
random_acc <- sum(true_fractions^2)
random_acc
## [1] 0.2943388
# The SVM's accuracy is higher than random guessing.

Classification and Regression Trees

  1. Read in the breast cancer imaging data “ispy1doctored.csv” into a data frame called dat. The iSPY1 trial was a prospective study to test if we could use MRIs to predict response to treatment and risk-of-recurrence in patients with stage 2 or 3 breast cancer receiving neoadjuvant chemotherapy (NACT).
ispy_dat <- read.csv("ispy1doctored.csv")
  1. Load the rpart package. Set the random seed to 10.
library(rpart)
set.seed(10)
  1. (1 point) Generate a histogram of the MRI_LD_Tfinal variable that will be our outcome to predict.
hist(ispy_dat$MRI_LD_Tfinal, xlab="MRI_LD_Tfinal", main="Histogram of MRI_LD_Tfinal", col="blue")

  1. (1 point) Split the dataset into a training set of size 70 and a testset consisting of the remaining data.
trainSet <- sample(1:nrow(ispy_dat),70)
ispy_datTrain <- ispy_dat[trainSet, ]
ispy_datTest <- ispy_dat[-trainSet, ]
  1. (3 points) Fit a regression tree to the training data with MRI_LD_Tfinal as outcome and all other variables as candidate predictors. Make sure that you specify the correct method for regression. Plot the fitted tree with text labels.
tree1 <- rpart(MRI_LD_Tfinal ~ ., method="anova", data=ispy_datTrain)

library(rpart.plot)
rpart.plot(tree1)

  1. (1 point) How many partitions does the fitted tree have?
# There are 4 partitions formed by the tree.
  1. (2 points) Determine the mean squared error (MSE) of the fitted tree based on the test data.
mean((ispy_datTest[ ,"MRI_LD_Tfinal"] - predict(tree1,newdata=ispy_datTest))^2)
## [1] 371.4719

Bagging

  1. Load the randomForest package and set the random seed to 4
library(randomForest)
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
## 
##     margin
set.seed(4)
  1. (2 points) Perform bagging with MRI_LD_Tfinal as outcome and all other variables as candidate predictors.
bag1 <- randomForest(MRI_LD_Tfinal ~ ., data=ispy_datTrain, mtry=16)
  1. (1 point) Determine the mean squared error (MSE) of the bagging based on the test data.
mean((ispy_datTest[ ,"MRI_LD_Tfinal"] - predict(bag1,newdata=ispy_datTest))^2)
## [1] 346.9801

Random Forests

  1. (2 points) Set the seed to 4 again and perform random forests with 6 candidate predictors per tree, MRI_LD_Tfinal as outcome and all other variables as candidate predictors.
set.seed(4)
rf <- randomForest(MRI_LD_Tfinal ~ ., data=ispy_datTrain, mtry=6, importance=T)
  1. (1 point) Determine the mean squared error (MSE) of the random forest based on the test data.
mean((ispy_datTest[ ,"MRI_LD_Tfinal"] - predict(rf, newdata=ispy_datTest))^2)
## [1] 349.0734
  1. (2 points) Plot the variable importance from the random forest. Use the importance measured in terms of the mean increase in mean squared error. Based on this importance measure, which variable is most important? Which variable is least important?
varImpPlot(rf, type=1)

# According to this plot, MRI_LD_T2 is most important and HR_HER2status is the least important.