Heart 9 - SVM (Classification)
Nao Mimoto - Dept. of Statistics : The University of Akron – Class Web Page
1. Load Data and Prepare
library(tidyverse)
Heart <- read_csv(file="https://nmimoto.github.io/datasets/heart.csv")
# Rename medv column as resp.
Heart2 <- Heart %>%
select(-"index") %>% # remove col named "index"
rename(resp=AHD) %>% # rename the column
relocate(resp) %>% # move "resp" to 1st column
mutate(resp=as.factor(resp), # Turn these columns to <factor> instead of <double>
ChestPain=as.factor(ChestPain),
Thal=as.factor(Thal),
Sex=as.factor(Sex),
Fbs=as.factor(Fbs),
RestECG=as.factor(RestECG),
ExAng=as.factor(ExAng))
print(Heart2, width=1000)## # A tibble: 303 × 14
## resp Age Sex ChestPain RestBP Chol Fbs RestECG MaxHR ExAng Oldpeak
## <fct> <dbl> <fct> <fct> <dbl> <dbl> <fct> <fct> <dbl> <fct> <dbl>
## 1 No 63 1 typical 145 233 1 2 150 0 2.3
## 2 Yes 67 1 asymptomatic 160 286 0 2 108 1 1.5
## 3 Yes 67 1 asymptomatic 120 229 0 2 129 1 2.6
## 4 No 37 1 nonanginal 130 250 0 0 187 0 3.5
## 5 No 41 0 nontypical 130 204 0 2 172 0 1.4
## 6 No 56 1 nontypical 120 236 0 0 178 0 0.8
## 7 Yes 62 0 asymptomatic 140 268 0 2 160 0 3.6
## 8 No 57 0 asymptomatic 120 354 0 0 163 1 0.6
## 9 Yes 63 1 asymptomatic 130 254 0 2 147 0 1.4
## 10 Yes 53 1 asymptomatic 140 203 1 2 155 1 3.1
## Slope Ca Thal
## <dbl> <dbl> <fct>
## 1 3 0 fixed
## 2 2 3 normal
## 3 2 2 reversable
## 4 3 0 normal
## 5 1 0 normal
## 6 1 0 normal
## 7 3 2 normal
## 8 1 0 normal
## 9 2 1 reversable
## 10 3 0 reversable
## # ℹ 293 more rows##
## No Yes
## 164 139######################################################
###
### When we are doing SVM or NN, we need to turn all variables into numbers.
### The data will be scaled inside svm() automatically.
###
######################################################
print(Heart, width=1000)## # A tibble: 303 × 15
## index Age Sex ChestPain RestBP Chol Fbs RestECG MaxHR ExAng Oldpeak
## <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 63 1 typical 145 233 1 2 150 0 2.3
## 2 2 67 1 asymptomatic 160 286 0 2 108 1 1.5
## 3 3 67 1 asymptomatic 120 229 0 2 129 1 2.6
## 4 4 37 1 nonanginal 130 250 0 0 187 0 3.5
## 5 5 41 0 nontypical 130 204 0 2 172 0 1.4
## 6 6 56 1 nontypical 120 236 0 0 178 0 0.8
## 7 7 62 0 asymptomatic 140 268 0 2 160 0 3.6
## 8 8 57 0 asymptomatic 120 354 0 0 163 1 0.6
## 9 9 63 1 asymptomatic 130 254 0 2 147 0 1.4
## 10 10 53 1 asymptomatic 140 203 1 2 155 1 3.1
## Slope Ca Thal AHD
## <dbl> <dbl> <chr> <chr>
## 1 3 0 fixed No
## 2 2 3 normal Yes
## 3 2 2 reversable Yes
## 4 3 0 normal No
## 5 1 0 normal No
## 6 1 0 normal No
## 7 3 2 normal Yes
## 8 1 0 normal No
## 9 2 1 reversable Yes
## 10 3 0 reversable Yes
## # ℹ 293 more rows# Rename medv column as resp.
Heart3 <- Heart %>%
select(-"index") %>% # remove col named "index"
rename(resp=AHD) %>% # rename the column
relocate(resp) %>% # move "resp" to 1st column
mutate(resp=as.factor(resp)) # Turn these columns to <factor> instead of <double>
print(Heart3, width=1000)## # A tibble: 303 × 14
## resp Age Sex ChestPain RestBP Chol Fbs RestECG MaxHR ExAng Oldpeak
## <fct> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 No 63 1 typical 145 233 1 2 150 0 2.3
## 2 Yes 67 1 asymptomatic 160 286 0 2 108 1 1.5
## 3 Yes 67 1 asymptomatic 120 229 0 2 129 1 2.6
## 4 No 37 1 nonanginal 130 250 0 0 187 0 3.5
## 5 No 41 0 nontypical 130 204 0 2 172 0 1.4
## 6 No 56 1 nontypical 120 236 0 0 178 0 0.8
## 7 Yes 62 0 asymptomatic 140 268 0 2 160 0 3.6
## 8 No 57 0 asymptomatic 120 354 0 0 163 1 0.6
## 9 Yes 63 1 asymptomatic 130 254 0 2 147 0 1.4
## 10 Yes 53 1 asymptomatic 140 203 1 2 155 1 3.1
## Slope Ca Thal
## <dbl> <dbl> <chr>
## 1 3 0 fixed
## 2 2 3 normal
## 3 2 2 reversable
## 4 3 0 normal
## 5 1 0 normal
## 6 1 0 normal
## 7 3 2 normal
## 8 1 0 normal
## 9 2 1 reversable
## 10 3 0 reversable
## # ℹ 293 more rows## # A tibble: 303 × 12
## resp Age Sex RestBP Chol Fbs RestECG MaxHR ExAng Oldpeak Slope Ca
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 No 63 1 145 233 1 2 150 0 2.3 3 0
## 2 Yes 67 1 160 286 0 2 108 1 1.5 2 3
## 3 Yes 67 1 120 229 0 2 129 1 2.6 2 2
## 4 No 37 1 130 250 0 0 187 0 3.5 3 0
## 5 No 41 0 130 204 0 2 172 0 1.4 1 0
## 6 No 56 1 120 236 0 0 178 0 0.8 1 0
## 7 Yes 62 0 140 268 0 2 160 0 3.6 3 2
## 8 No 57 0 120 354 0 0 163 1 0.6 1 0
## 9 Yes 63 1 130 254 0 2 147 0 1.4 2 1
## 10 Yes 53 1 140 203 1 2 155 1 3.1 3 0
## # ℹ 293 more rows##----------
## remove rows with NA.
Orig <- Heart3
# Check for N/A in data. Remove if there's any.
summary(Orig)## resp Age Sex RestBP Chol
## No :164 Min. :29.00 Min. :0.0000 Min. : 94.0 Min. :126.0
## Yes:139 1st Qu.:48.00 1st Qu.:0.0000 1st Qu.:120.0 1st Qu.:211.0
## Median :56.00 Median :1.0000 Median :130.0 Median :241.0
## Mean :54.44 Mean :0.6799 Mean :131.7 Mean :246.7
## 3rd Qu.:61.00 3rd Qu.:1.0000 3rd Qu.:140.0 3rd Qu.:275.0
## Max. :77.00 Max. :1.0000 Max. :200.0 Max. :564.0
##
## Fbs RestECG MaxHR ExAng
## Min. :0.0000 Min. :0.0000 Min. : 71.0 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:133.5 1st Qu.:0.0000
## Median :0.0000 Median :1.0000 Median :153.0 Median :0.0000
## Mean :0.1485 Mean :0.9901 Mean :149.6 Mean :0.3267
## 3rd Qu.:0.0000 3rd Qu.:2.0000 3rd Qu.:166.0 3rd Qu.:1.0000
## Max. :1.0000 Max. :2.0000 Max. :202.0 Max. :1.0000
##
## Oldpeak Slope Ca
## Min. :0.00 Min. :1.000 Min. :0.0000
## 1st Qu.:0.00 1st Qu.:1.000 1st Qu.:0.0000
## Median :0.80 Median :2.000 Median :0.0000
## Mean :1.04 Mean :1.601 Mean :0.6722
## 3rd Qu.:1.60 3rd Qu.:2.000 3rd Qu.:1.0000
## Max. :6.20 Max. :3.000 Max. :3.0000
## NA's :4## [1] 303 12## [1] 4## [1] 299 122. Data Separation
# Load the file from my website
source("https://nmimoto.github.io/R/ML-00.txt") # load CreateCV()
# Use this function to create CV dataset
my.seed = 131
CreateCV(Orig, numFolds=5, seed=my.seed)##
## number of folds: 5
## number of obs: 299
## training data: 49 x 5 = 245
## size of test data: 54
##
## Global variables created:
## Train.set Train.resp
## Test.set Test.resp
## CV.train[[k]] CV.train.resp[[k]]
## CV.valid[[k]] CV.valid.resp[[k]]# Output (all as Data Frame):
# Train.set Train.resp
# Test.set Test.resp
# CV.train[[[[k]]]] CV.train.resp[[[[k]]]]
# CV.valid[[[[k]]]] CV.valid.resp[[[[k]]]]3. SVM with Linear Kernel using In-House CV
library (e1071) # install.packages('e1071", repos='https://cran.case.edu/')
AUCs <- MSE.valid <- matrix(0, 5, 2)
colnames(AUCs) = c("Train AUC", "Valid AUC")
for (k in 1:5) {
Fit00 <- e1071::svm(resp~., data=CV.train[[k]], # <----- Change model here
kernel="linear", cost=5,
# kernel ="radial", gamma = .5, cost=1,
scale=TRUE)
summary(Fit00)
#- Training Confusion Matrix
Train.prob0 = predict(Fit00, CV.train[[k]], decision.values=TRUE)
Train.prob = attributes(Train.prob0)$decision.values
#- Predict in Validation Set
Valid.prob0 = predict(Fit00, CV.valid[[k]], decision.values=TRUE)
Valid.prob = attributes(Valid.prob0)$decision.values
#- Check the training set accuracy
library(caret)
library(pROC)
Train.prob1 = as.numeric(Train.prob)# (Train.prob-min(Train.prob))/(max(Train.prob)-min(Train.prob)))
Valid.prob1 = as.numeric(Train.prob) #(Valid.prob-min(Valid.prob))/(max(Valid.prob)-min(Valid.prob)))
AUCs[k,] <- round(c(auc(factor(as.matrix(CV.train.resp[[k]])), as.vector(Train.prob), levels=c("No", "Yes")),
auc(factor(as.matrix(CV.valid.resp[[k]])), as.vector(Valid.prob), levels=c("No", "Yes"))), 4)
}
AUCs## Train AUC Valid AUC
## [1,] 0.9328 0.7224
## [2,] 0.8834 0.9478
## [3,] 0.8953 0.9024
## [4,] 0.9110 0.7642
## [5,] 0.8957 0.9259## Av.Train AUC Av.Valid AUC
## 0.90364 0.85254###
### Make decision about best model based on Av Valid AUC.
### Better Cost parameter value must be found.
### After Best Cost value is found, final Train/Test fit must be performed.
###4. SVM with Linear Kernel using tune() on Train.set
library (e1071) # install.packages('e1071", repos='https://cran.case.edu/')
# Pick best Cost with automatic 5-fold CV
set.seed (my.seed)
Tuned01 = e1071::tune(svm, resp~., data=Train.set, kernel ="linear",
ranges=list(
cost=c(0.01, 0.1, 1, 5, 10, 100, 1000)), # <- try these cost parameters
scale=TRUE,
tunecontrol=tune.control(cross=5))
summary(Tuned01)##
## Parameter tuning of 'svm':
##
## - sampling method: 5-fold cross validation
##
## - best parameters:
## cost
## 0.1
##
## - best performance: 0.2040816
##
## - Detailed performance results:
## cost error dispersion
## 1 1e-02 0.2122449 0.07301446
## 2 1e-01 0.2040816 0.06920745
## 3 1e+00 0.2122449 0.06860301
## 4 5e+00 0.2122449 0.06388766
## 5 1e+01 0.2163265 0.06549717
## 6 1e+02 0.2163265 0.06549717
## 7 1e+03 0.2163265 0.06549717##
## Call:
## best.tune(METHOD = svm, train.x = resp ~ ., data = Train.set, ranges = list(cost = c(0.01,
## 0.1, 1, 5, 10, 100, 1000)), tunecontrol = tune.control(cross = 5),
## kernel = "linear", scale = TRUE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 0.1
##
## Number of Support Vectors: 119
##
## ( 60 59 )
##
##
## Number of Classes: 2
##
## Levels:
## No YesFinal Fit with Training Set vs Test Set using the best model
Final.fit <- Tuned01$best.model # result from the auto CV tuning
# Training Confusion Matrix
Train.prob0 = predict(Final.fit, Train.set, decision.values=TRUE)
Train.prob = attributes(Train.prob0)$decision.values
head(Train.prob) # Check col name for No/Yes or Yes/No## No/Yes
## 1 0.1415496
## 2 -3.8959628
## 3 -2.4653451
## 4 -0.2926215
## 5 1.5751539
## 6 1.5569651# Testing Confusion Matrix
Test.prob0 = predict(Final.fit, Test.set, decision.values=TRUE)
Test.prob = attributes(Test.prob0)$decision.values
head(Train.prob) # Check col name for No/Yes or Yes/No## No/Yes
## 1 0.1415496
## 2 -3.8959628
## 3 -2.4653451
## 4 -0.2926215
## 5 1.5751539
## 6 1.5569651# Pick a threshold
threshold <- .7
Train.pred = ifelse(Train.prob > threshold, "No", "Yes")
#Train.pred = ifelse(Train.prob > threshold, "Yes", "No") #-
Test.pred = ifelse(Test.prob > threshold, "No", "Yes")
#Test.pred = ifelse(Test.prob > threshold, "Yes", "No")
CM.train <- caret::confusionMatrix(factor(Train.pred), factor(as.matrix(Train.resp)), positive="Yes")
CM.train## Confusion Matrix and Statistics
##
## Reference
## Prediction No Yes
## No 95 11
## Yes 41 98
##
## Accuracy : 0.7878
## 95% CI : (0.7312, 0.8372)
## No Information Rate : 0.5551
## P-Value [Acc > NIR] : 2.201e-14
##
## Kappa : 0.5817
##
## Mcnemar's Test P-Value : 5.781e-05
##
## Sensitivity : 0.8991
## Specificity : 0.6985
## Pos Pred Value : 0.7050
## Neg Pred Value : 0.8962
## Prevalence : 0.4449
## Detection Rate : 0.4000
## Detection Prevalence : 0.5673
## Balanced Accuracy : 0.7988
##
## 'Positive' Class : Yes
## CM.test <- caret::confusionMatrix(factor(Test.pred), factor(as.matrix(Test.resp)), positive="Yes")
CM.test## Confusion Matrix and Statistics
##
## Reference
## Prediction No Yes
## No 17 4
## Yes 8 25
##
## Accuracy : 0.7778
## 95% CI : (0.644, 0.8796)
## No Information Rate : 0.537
## P-Value [Acc > NIR] : 0.0002203
##
## Kappa : 0.5481
##
## Mcnemar's Test P-Value : 0.3864762
##
## Sensitivity : 0.8621
## Specificity : 0.6800
## Pos Pred Value : 0.7576
## Neg Pred Value : 0.8095
## Prevalence : 0.5370
## Detection Rate : 0.4630
## Detection Prevalence : 0.6111
## Balanced Accuracy : 0.7710
##
## 'Positive' Class : Yes
## # ROC curve and AUC
layout(matrix(1:2,1,2))
pROC::plot.roc(factor(as.matrix(Train.resp)), as.vector(Train.prob), levels=c("No", "Yes"))
# point corresponding to CM.train
abline(h=CM.train[["byClass"]][["Sensitivity"]], v=1-CM.train[["byClass"]][["Specificity"]], col="red")
auc.train = pROC::auc(factor(as.matrix(Train.resp)), as.vector(Train.prob), levels=c("No", "Yes"))
text(.2, .2, paste("AUC=",round(auc.train, 3)))
# Testing
pROC::plot.roc(factor(as.matrix(Test.resp)),
as.vector(Test.prob), levels=c("No", "Yes"))
# point corresponding to CM.train
abline(h=CM.test[["byClass"]][["Sensitivity"]], v=1-CM.test[["byClass"]][["Specificity"]], col="red")
auc.test = pROC::auc(factor(as.matrix(Test.resp)), as.vector(Test.prob), levels=c("No", "Yes"))
text(.2, .2, paste("AUC=",round(auc.test, 3)))
## [1] 0.8973287 0.8979310###
### AUC.test is the estimate for this model's prediction power. (with threshold = .5)
### ROC curve represents the overall capability of this model.
### By changing threshold, we can be at anywhere on the curve.
###5. SVM with Radial Kernel using Auto-CV on Train.set
library (e1071) # install.packages('e1071", repos='https://cran.case.edu/')
# Pick best Cost with automatic 5-fold CV
set.seed (my.seed)
Tuned02 = e1071::tune(svm, resp~., data=Train.set, kernel ="radial",
ranges=list(
gamma = 2^(-1:4),
cost=c(0.001, 0.01, 0.1, 1, 5, 10, 100, 1000)),
scale=TRUE,
tunecontrol=tune.control(cross=5))
summary(Tuned02)##
## Parameter tuning of 'svm':
##
## - sampling method: 5-fold cross validation
##
## - best parameters:
## gamma cost
## 0.5 5
##
## - best performance: 0.244898
##
## - Detailed performance results:
## gamma cost error dispersion
## 1 0.5 1e-03 0.4448980 0.07272869
## 2 1.0 1e-03 0.4448980 0.07272869
## 3 2.0 1e-03 0.4448980 0.07272869
## 4 4.0 1e-03 0.4448980 0.07272869
## 5 8.0 1e-03 0.4448980 0.07272869
## 6 16.0 1e-03 0.4448980 0.07272869
## 7 0.5 1e-02 0.4448980 0.07272869
## 8 1.0 1e-02 0.4448980 0.07272869
## 9 2.0 1e-02 0.4448980 0.07272869
## 10 4.0 1e-02 0.4448980 0.07272869
## 11 8.0 1e-02 0.4448980 0.07272869
## 12 16.0 1e-02 0.4448980 0.07272869
## 13 0.5 1e-01 0.4448980 0.07272869
## 14 1.0 1e-01 0.4448980 0.07272869
## 15 2.0 1e-01 0.4448980 0.07272869
## 16 4.0 1e-01 0.4448980 0.07272869
## 17 8.0 1e-01 0.4448980 0.07272869
## 18 16.0 1e-01 0.4448980 0.07272869
## 19 0.5 1e+00 0.2693878 0.06356086
## 20 1.0 1e+00 0.4081633 0.09462876
## 21 2.0 1e+00 0.4489796 0.07215375
## 22 4.0 1e+00 0.4448980 0.07272869
## 23 8.0 1e+00 0.4448980 0.07272869
## 24 16.0 1e+00 0.4448980 0.07272869
## 25 0.5 5e+00 0.2448980 0.04786139
## 26 1.0 5e+00 0.3755102 0.09851035
## 27 2.0 5e+00 0.4448980 0.06980666
## 28 4.0 5e+00 0.4489796 0.07215375
## 29 8.0 5e+00 0.4448980 0.07272869
## 30 16.0 5e+00 0.4448980 0.07272869
## 31 0.5 1e+01 0.2448980 0.04786139
## 32 1.0 1e+01 0.3755102 0.09851035
## 33 2.0 1e+01 0.4448980 0.06980666
## 34 4.0 1e+01 0.4489796 0.07215375
## 35 8.0 1e+01 0.4448980 0.07272869
## 36 16.0 1e+01 0.4448980 0.07272869
## 37 0.5 1e+02 0.2448980 0.04786139
## 38 1.0 1e+02 0.3755102 0.09851035
## 39 2.0 1e+02 0.4448980 0.06980666
## 40 4.0 1e+02 0.4489796 0.07215375
## 41 8.0 1e+02 0.4448980 0.07272869
## 42 16.0 1e+02 0.4448980 0.07272869
## 43 0.5 1e+03 0.2448980 0.04786139
## 44 1.0 1e+03 0.3755102 0.09851035
## 45 2.0 1e+03 0.4448980 0.06980666
## 46 4.0 1e+03 0.4489796 0.07215375
## 47 8.0 1e+03 0.4448980 0.07272869
## 48 16.0 1e+03 0.4448980 0.07272869##
## Call:
## best.tune(METHOD = svm, train.x = resp ~ ., data = Train.set, ranges = list(gamma = 2^(-1:4),
## cost = c(0.001, 0.01, 0.1, 1, 5, 10, 100, 1000)), tunecontrol = tune.control(cross = 5),
## kernel = "radial", scale = TRUE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 5
##
## Number of Support Vectors: 226
##
## ( 120 106 )
##
##
## Number of Classes: 2
##
## Levels:
## No YesFinal Fit with Training Set vs Test Set using the best model
Final.fit <- Tuned02$best.model # results from CV tuning
# Training Confusion Matrix
Train.prob0 = predict(Final.fit, Train.set, decision.values=TRUE)
Train.prob = attributes(Train.prob0)$decision.values
head(Train.prob) # Check col name for No/Yes or Yes/No## No/Yes
## 1 0.9999361
## 2 -1.0001508
## 3 -0.9998412
## 4 1.0001208
## 5 1.0000147
## 6 0.9999346# Testing Confusion Matrix
Test.prob0 = predict(Final.fit, Test.set, decision.values=TRUE)
Test.prob = attributes(Test.prob0)$decision.values
head(Train.prob) # Check col name for No/Yes or Yes/No## No/Yes
## 1 0.9999361
## 2 -1.0001508
## 3 -0.9998412
## 4 1.0001208
## 5 1.0000147
## 6 0.9999346# Pick a threshold
threshold <- .7
Train.pred = ifelse(Train.prob > threshold, "No", "Yes")
#Train.pred = ifelse(Train.prob > threshold, "Yes", "No") #-
Test.pred = ifelse(Test.prob > threshold, "No", "Yes")
#Test.pred = ifelse(Test.prob > threshold, "Yes", "No")
CM.train <- caret::confusionMatrix(factor(Train.pred), factor(as.matrix(Train.resp)), positive="Yes")
CM.train## Confusion Matrix and Statistics
##
## Reference
## Prediction No Yes
## No 136 0
## Yes 0 109
##
## Accuracy : 1
## 95% CI : (0.9851, 1)
## No Information Rate : 0.5551
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 1
##
## Mcnemar's Test P-Value : NA
##
## Sensitivity : 1.0000
## Specificity : 1.0000
## Pos Pred Value : 1.0000
## Neg Pred Value : 1.0000
## Prevalence : 0.4449
## Detection Rate : 0.4449
## Detection Prevalence : 0.4449
## Balanced Accuracy : 1.0000
##
## 'Positive' Class : Yes
## CM.test <- caret::confusionMatrix(factor(Test.pred), factor(as.matrix(Test.resp)), positive="Yes")
CM.test## Confusion Matrix and Statistics
##
## Reference
## Prediction No Yes
## No 10 0
## Yes 15 29
##
## Accuracy : 0.7222
## 95% CI : (0.5836, 0.8354)
## No Information Rate : 0.537
## P-Value [Acc > NIR] : 0.0041711
##
## Kappa : 0.4173
##
## Mcnemar's Test P-Value : 0.0003006
##
## Sensitivity : 1.0000
## Specificity : 0.4000
## Pos Pred Value : 0.6591
## Neg Pred Value : 1.0000
## Prevalence : 0.5370
## Detection Rate : 0.5370
## Detection Prevalence : 0.8148
## Balanced Accuracy : 0.7000
##
## 'Positive' Class : Yes
## #- ROC curve and AUC
layout(matrix(1:2,1,2))
pROC::plot.roc(factor(as.matrix(Train.resp)), as.vector(Train.prob), levels=c("No", "Yes"))
# point corresponding to CM.train
abline(h=CM.train[["byClass"]][["Sensitivity"]], v=1-CM.train[["byClass"]][["Specificity"]], col="red")
auc.train = pROC::auc(factor(as.matrix(Train.resp)), as.vector(Train.prob), levels=c("No", "Yes"))
text(.2, .2, paste("AUC=",round(auc.train, 3)))
#- Testing
pROC::plot.roc(factor(as.matrix(Test.resp)),
as.vector(Test.prob), levels=c("No", "Yes"))
# point corresponding to CM.train
abline(h=CM.test[["byClass"]][["Sensitivity"]], v=1-CM.test[["byClass"]][["Specificity"]], col="red")
auc.test = pROC::auc(factor(as.matrix(Test.resp)), as.vector(Test.prob), levels=c("No", "Yes"))
text(.2, .2, paste("AUC=",round(auc.test, 3)))
## [1] 1.0000000 0.88551726. SVM with Polynomial Kernel using Auto-CV on Train.set
library (e1071) # install.packages('e1071", repos='https://cran.case.edu/')
# Pick best Cost with automatic 5-fold CV
set.seed (my.seed)
Tuned03 = e1071::tune(svm, resp~., data=Train.set, kernel ="polynomial",
ranges=list(
degree = 2^(-1:3),
gamma = 2^(2:5),
coef0=c(-2, 0,2^(-1:2)),
cost=c(0.0001, 0.001, 0.01, 0.1, 1, 5, 10, 100)),
scale=TRUE,
tunecontrol=tune.control(cross=5))
summary(Tuned03)
#--- Extract the model with best cost
# plot(resp~., Tuned03$best.model, data=Train.set)
summary(Tuned03$best.model)Final Fit with Training Set vs Test Set using the best model
Final.fit <- Tuned03$best.model # results from CV tuning
# Training Confusion Matrix
Train.prob0 = predict(Final.fit, Train.set, decision.values=TRUE)
Train.prob = attributes(Train.prob0)$decision.values
head(Train.prob) # Check col name for No/Yes or Yes/No
# Testing Confusion Matrix
Test.prob0 = predict(Final.fit, Test.set, decision.values=TRUE)
Test.prob = attributes(Test.prob0)$decision.values
head(Train.prob) # Check col name for No/Yes or Yes/No
# Pick a threshold
threshold <- .7
Train.pred = ifelse(Train.prob > threshold, "No", "Yes")
#Train.pred = ifelse(Train.prob > threshold, "Yes", "No") #-
Test.pred = ifelse(Test.prob > threshold, "No", "Yes")
#Test.pred = ifelse(Test.prob > threshold, "Yes", "No")
CM.train <- caret::confusionMatrix(factor(Train.pred), factor(as.matrix(Train.resp)), positive="Yes")
CM.train
CM.test <- caret::confusionMatrix(factor(Test.pred), factor(as.matrix(Test.resp)), positive="Yes")
CM.test
#- ROC curve and AUC
layout(matrix(1:2,1,2))
pROC::plot.roc(factor(as.matrix(Train.resp)), as.vector(Train.prob), levels=c("No", "Yes"))
# point corresponding to CM.train
abline(h=CM.train[["byClass"]][["Sensitivity"]], v=1-CM.train[["byClass"]][["Specificity"]], col="red")
auc.train = pROC::auc(factor(as.matrix(Train.resp)), as.vector(Train.prob), levels=c("No", "Yes"))
text(.2, .2, paste("AUC=",round(auc.train, 3)))
#- Testing
pROC::plot.roc(factor(as.matrix(Test.resp)),
as.vector(Test.prob), levels=c("No", "Yes"))
# point corresponding to CM.train
abline(h=CM.test[["byClass"]][["Sensitivity"]], v=1-CM.test[["byClass"]][["Specificity"]], col="red")
auc.test = pROC::auc(factor(as.matrix(Test.resp)), as.vector(Test.prob), levels=c("No", "Yes"))
text(.2, .2, paste("AUC=",round(auc.test, 3)))
layout(1)
c(auc.train, auc.test)7. Threshold Optimization
#cost.list = c(0,0,3,1)/4 # order of (TP, TN, FP, FN)
#cost.list = c(0,0,2,1)/3 # order of (TP, TN, FP, FN)
#cost.list = c(0,0,1,1)/2 # order of (TP, TN, FP, FN)
#cost.list = c(0,0,1,2)/3 # order of (TP, TN, FP, FN)
cost.list = c(0,0,1,3)/4 # order of (TP, TN, FP, FN)
threshold.list = seq(.01,.99,.01) # grid for threshold
#- Pick a model
svm.fit <- Tuned02$best.model
Test.prob0 = predict(svm.fit, Test.set, decision.values=TRUE)
Test.prob = attributes(Test.prob0)$decision.values
head(Test.pred) # Check col name for No/Yes or Yes/No
library(caret) # for confusionMatrix
cost=0
for (i in 1:length(threshold.list)){
threshold = threshold.list[i]
#- Check the training set accuracy
Test.pred = ifelse(Test.prob > threshold, "Yes", "No")
CM.test <- confusionMatrix(factor(Test.pred),
factor(as.matrix(Test.resp)),
positive="Yes")
TP = CM.test$table[2,2] # True Pos
TN = CM.test$table[1,1] # True Neg
FP = CM.test$table[2,1] # False Pos
FN = CM.test$table[1,2] # False Neg
cost[i] = sum(c(TP, TN, FP, FN) * cost.list)
}
plot(threshold.list, cost, xlab="threshold")
cost.list
which.min(cost)
min(cost)
threshold.list[which.min(cost)]CODE ONLY
##--------------------------------------
## 1. Load Data and Prepare
library(tidyverse)
Heart <- read_csv(file="https://nmimoto.github.io/datasets/heart.csv")
# Rename medv column as resp.
Heart2 <- Heart %>%
select(-"index") %>% # remove col named "index"
rename(resp=AHD) %>% # rename the column
relocate(resp) %>% # move "resp" to 1st column
mutate(resp=as.factor(resp), # Turn these columns to <factor> instead of <double>
ChestPain=as.factor(ChestPain),
Thal=as.factor(Thal),
Sex=as.factor(Sex),
Fbs=as.factor(Fbs),
RestECG=as.factor(RestECG),
ExAng=as.factor(ExAng))
print(Heart2, width=1000)
table(Heart2$resp) # note that first level is "No"
######################################################
###
### When we are doing SVM or NN, we need to turn all variables into numbers.
### The data will be scaled inside svm() automatically.
###
######################################################
print(Heart, width=1000)
# Rename medv column as resp.
Heart3 <- Heart %>%
select(-"index") %>% # remove col named "index"
rename(resp=AHD) %>% # rename the column
relocate(resp) %>% # move "resp" to 1st column
mutate(resp=as.factor(resp)) # Turn these columns to <factor> instead of <double>
print(Heart3, width=1000)
Heart3 <- Heart3 %>% select(-c("ChestPain", "Thal")) # remove columns
print(Heart3, width=1000)
##----------
## remove rows with NA.
Orig <- Heart3
# Check for N/A in data. Remove if there's any.
summary(Orig)
dim(Orig)
sum(is.na(Orig))
# If there is na in the data, run below
Orig <- Orig %>% na.omit()
dim(Orig)
##----------------------------------------
## 2. Data Separation
# Load the file from my website
source("https://nmimoto.github.io/R/ML-00.txt") # load CreateCV()
# Use this function to create CV dataset
my.seed = 131
CreateCV(Orig, numFolds=5, seed=my.seed)
# Output (all as Data Frame):
# Train.set Train.resp
# Test.set Test.resp
# CV.train[[[[k]]]] CV.train.resp[[[[k]]]]
# CV.valid[[[[k]]]] CV.valid.resp[[[[k]]]]
##----------------------------------------
## 3. SVM with Linear Kernel using In-House CV
library (e1071) # install.packages('e1071", repos='https://cran.case.edu/')
AUCs <- MSE.valid <- matrix(0, 5, 2)
colnames(AUCs) = c("Train AUC", "Valid AUC")
for (k in 1:5) {
Fit00 <- e1071::svm(resp~., data=CV.train[[k]], # <----- Change model here
kernel="linear", cost=5,
# kernel ="radial", gamma = .5, cost=1,
scale=TRUE)
summary(Fit00)
#- Training Confusion Matrix
Train.prob0 = predict(Fit00, CV.train[[k]], decision.values=TRUE)
Train.prob = attributes(Train.prob0)$decision.values
#- Predict in Validation Set
Valid.prob0 = predict(Fit00, CV.valid[[k]], decision.values=TRUE)
Valid.prob = attributes(Valid.prob0)$decision.values
#- Check the training set accuracy
library(caret)
library(pROC)
Train.prob1 = as.numeric(Train.prob)# (Train.prob-min(Train.prob))/(max(Train.prob)-min(Train.prob)))
Valid.prob1 = as.numeric(Train.prob) #(Valid.prob-min(Valid.prob))/(max(Valid.prob)-min(Valid.prob)))
AUCs[k,] <- round(c(auc(factor(as.matrix(CV.train.resp[[k]])), as.vector(Train.prob), levels=c("No", "Yes")),
auc(factor(as.matrix(CV.valid.resp[[k]])), as.vector(Valid.prob), levels=c("No", "Yes"))), 4)
}
AUCs
Av.AUCs = apply(AUCs, 2, mean)
names(Av.AUCs) = c("Av.Train AUC", "Av.Valid AUC")
Av.AUCs
###
### Make decision about best model based on Av Valid AUC.
### Better Cost parameter value must be found.
### After Best Cost value is found, final Train/Test fit must be performed.
###
##----------------------------------------
## 4. SVM with Linear Kernel using tune() on Train.set
library (e1071) # install.packages('e1071", repos='https://cran.case.edu/')
# Pick best Cost with automatic 5-fold CV
set.seed (my.seed)
Tuned01 = e1071::tune(svm, resp~., data=Train.set, kernel ="linear",
ranges=list(
cost=c(0.01, 0.1, 1, 5, 10, 100, 1000)), # <- try these cost parameters
scale=TRUE,
tunecontrol=tune.control(cross=5))
summary(Tuned01)
#--- Extract the model with best cost
summary(Tuned01$best.model)
# plot(resp~., Tuned01$best.model, data=Train.set)
###------------------
### Final Fit with Training Set vs Test Set using the best model
Final.fit <- Tuned01$best.model # result from the auto CV tuning
# Training Confusion Matrix
Train.prob0 = predict(Final.fit, Train.set, decision.values=TRUE)
Train.prob = attributes(Train.prob0)$decision.values
head(Train.prob) # Check col name for No/Yes or Yes/No
# Testing Confusion Matrix
Test.prob0 = predict(Final.fit, Test.set, decision.values=TRUE)
Test.prob = attributes(Test.prob0)$decision.values
head(Train.prob) # Check col name for No/Yes or Yes/No
# Pick a threshold
threshold <- .7
Train.pred = ifelse(Train.prob > threshold, "No", "Yes")
#Train.pred = ifelse(Train.prob > threshold, "Yes", "No") #-
Test.pred = ifelse(Test.prob > threshold, "No", "Yes")
#Test.pred = ifelse(Test.prob > threshold, "Yes", "No")
CM.train <- caret::confusionMatrix(factor(Train.pred), factor(as.matrix(Train.resp)), positive="Yes")
CM.train
CM.test <- caret::confusionMatrix(factor(Test.pred), factor(as.matrix(Test.resp)), positive="Yes")
CM.test
# ROC curve and AUC
layout(matrix(1:2,1,2))
pROC::plot.roc(factor(as.matrix(Train.resp)), as.vector(Train.prob), levels=c("No", "Yes"))
# point corresponding to CM.train
abline(h=CM.train[["byClass"]][["Sensitivity"]], v=1-CM.train[["byClass"]][["Specificity"]], col="red")
auc.train = pROC::auc(factor(as.matrix(Train.resp)), as.vector(Train.prob), levels=c("No", "Yes"))
text(.2, .2, paste("AUC=",round(auc.train, 3)))
# Testing
pROC::plot.roc(factor(as.matrix(Test.resp)),
as.vector(Test.prob), levels=c("No", "Yes"))
# point corresponding to CM.train
abline(h=CM.test[["byClass"]][["Sensitivity"]], v=1-CM.test[["byClass"]][["Specificity"]], col="red")
auc.test = pROC::auc(factor(as.matrix(Test.resp)), as.vector(Test.prob), levels=c("No", "Yes"))
text(.2, .2, paste("AUC=",round(auc.test, 3)))
layout(1)
c(auc.train, auc.test)
###
### AUC.test is the estimate for this model's prediction power. (with threshold = .5)
### ROC curve represents the overall capability of this model.
### By changing threshold, we can be at anywhere on the curve.
###
##----------------------------------------
## 5. SVM with Radial Kernel using Auto-CV on Train.set
library (e1071) # install.packages('e1071", repos='https://cran.case.edu/')
# Pick best Cost with automatic 5-fold CV
set.seed (my.seed)
Tuned02 = e1071::tune(svm, resp~., data=Train.set, kernel ="radial",
ranges=list(
gamma = 2^(-1:4),
cost=c(0.001, 0.01, 0.1, 1, 5, 10, 100, 1000)),
scale=TRUE,
tunecontrol=tune.control(cross=5))
summary(Tuned02)
#--- Extract the model with best cost
summary(Tuned02$best.model)
# plot(resp~., Tuned02$best.model, data=Train.set)
###----------
### Final Fit with Training Set vs Test Set using the best model
Final.fit <- Tuned02$best.model # results from CV tuning
# Training Confusion Matrix
Train.prob0 = predict(Final.fit, Train.set, decision.values=TRUE)
Train.prob = attributes(Train.prob0)$decision.values
head(Train.prob) # Check col name for No/Yes or Yes/No
# Testing Confusion Matrix
Test.prob0 = predict(Final.fit, Test.set, decision.values=TRUE)
Test.prob = attributes(Test.prob0)$decision.values
head(Train.prob) # Check col name for No/Yes or Yes/No
# Pick a threshold
threshold <- .7
Train.pred = ifelse(Train.prob > threshold, "No", "Yes")
#Train.pred = ifelse(Train.prob > threshold, "Yes", "No") #-
Test.pred = ifelse(Test.prob > threshold, "No", "Yes")
#Test.pred = ifelse(Test.prob > threshold, "Yes", "No")
CM.train <- caret::confusionMatrix(factor(Train.pred), factor(as.matrix(Train.resp)), positive="Yes")
CM.train
CM.test <- caret::confusionMatrix(factor(Test.pred), factor(as.matrix(Test.resp)), positive="Yes")
CM.test
#- ROC curve and AUC
layout(matrix(1:2,1,2))
pROC::plot.roc(factor(as.matrix(Train.resp)), as.vector(Train.prob), levels=c("No", "Yes"))
# point corresponding to CM.train
abline(h=CM.train[["byClass"]][["Sensitivity"]], v=1-CM.train[["byClass"]][["Specificity"]], col="red")
auc.train = pROC::auc(factor(as.matrix(Train.resp)), as.vector(Train.prob), levels=c("No", "Yes"))
text(.2, .2, paste("AUC=",round(auc.train, 3)))
#- Testing
pROC::plot.roc(factor(as.matrix(Test.resp)),
as.vector(Test.prob), levels=c("No", "Yes"))
# point corresponding to CM.train
abline(h=CM.test[["byClass"]][["Sensitivity"]], v=1-CM.test[["byClass"]][["Specificity"]], col="red")
auc.test = pROC::auc(factor(as.matrix(Test.resp)), as.vector(Test.prob), levels=c("No", "Yes"))
text(.2, .2, paste("AUC=",round(auc.test, 3)))
layout(1)
c(auc.train, auc.test)
##----------------------------------------
## 6. SVM with Polynomial Kernel using Auto-CV on Train.set
library (e1071) # install.packages('e1071", repos='https://cran.case.edu/')
# Pick best Cost with automatic 5-fold CV
set.seed (my.seed)
#Tuned03 = e1071::tune(svm, resp~., data=Train.set, kernel ="polynomial",
# ranges=list(
# degree = 2^(-1:3),
# gamma = 2^(2:5),
# coef0=c(-2, 0,2^(-1:2)),
# cost=c(0.0001, 0.001, 0.01, 0.1, 1, 5, 10, 100)),
# scale=TRUE,
# tunecontrol=tune.control(cross=5))
#summary(Tuned03)
#--- Extract the model with best cost
# plot(resp~., Tuned03$best.model, data=Train.set)
summary(Tuned03$best.model)
###----------
### Final Fit with Training Set vs Test Set using the best model
Final.fit <- Tuned03$best.model # results from CV tuning
# Training Confusion Matrix
Train.prob0 = predict(Final.fit, Train.set, decision.values=TRUE)
Train.prob = attributes(Train.prob0)$decision.values
head(Train.prob) # Check col name for No/Yes or Yes/No
# Testing Confusion Matrix
Test.prob0 = predict(Final.fit, Test.set, decision.values=TRUE)
Test.prob = attributes(Test.prob0)$decision.values
head(Train.prob) # Check col name for No/Yes or Yes/No
# Pick a threshold
threshold <- .7
Train.pred = ifelse(Train.prob > threshold, "No", "Yes")
#Train.pred = ifelse(Train.prob > threshold, "Yes", "No") #-
Test.pred = ifelse(Test.prob > threshold, "No", "Yes")
#Test.pred = ifelse(Test.prob > threshold, "Yes", "No")
CM.train <- caret::confusionMatrix(factor(Train.pred), factor(as.matrix(Train.resp)), positive="Yes")
CM.train
CM.test <- caret::confusionMatrix(factor(Test.pred), factor(as.matrix(Test.resp)), positive="Yes")
CM.test
#- ROC curve and AUC
layout(matrix(1:2,1,2))
pROC::plot.roc(factor(as.matrix(Train.resp)), as.vector(Train.prob), levels=c("No", "Yes"))
# point corresponding to CM.train
abline(h=CM.train[["byClass"]][["Sensitivity"]], v=1-CM.train[["byClass"]][["Specificity"]], col="red")
auc.train = pROC::auc(factor(as.matrix(Train.resp)), as.vector(Train.prob), levels=c("No", "Yes"))
text(.2, .2, paste("AUC=",round(auc.train, 3)))
#- Testing
pROC::plot.roc(factor(as.matrix(Test.resp)),
as.vector(Test.prob), levels=c("No", "Yes"))
# point corresponding to CM.train
abline(h=CM.test[["byClass"]][["Sensitivity"]], v=1-CM.test[["byClass"]][["Specificity"]], col="red")
auc.test = pROC::auc(factor(as.matrix(Test.resp)), as.vector(Test.prob), levels=c("No", "Yes"))
text(.2, .2, paste("AUC=",round(auc.test, 3)))
layout(1)
c(auc.train, auc.test)
##----------------------------------------
## 7. Threshold Optimization
#cost.list = c(0,0,3,1)/4 # order of (TP, TN, FP, FN)
#cost.list = c(0,0,2,1)/3 # order of (TP, TN, FP, FN)
#cost.list = c(0,0,1,1)/2 # order of (TP, TN, FP, FN)
#cost.list = c(0,0,1,2)/3 # order of (TP, TN, FP, FN)
cost.list = c(0,0,1,3)/4 # order of (TP, TN, FP, FN)
threshold.list = seq(.01,.99,.01) # grid for threshold
#- Pick a model
svm.fit <- Tuned02$best.model
Test.prob0 = predict(svm.fit, Test.set, decision.values=TRUE)
Test.prob = attributes(Test.prob0)$decision.values
head(Test.pred) # Check col name for No/Yes or Yes/No
library(caret) # for confusionMatrix
cost=0
for (i in 1:length(threshold.list)){
threshold = threshold.list[i]
#- Check the training set accuracy
Test.pred = ifelse(Test.prob > threshold, "Yes", "No")
CM.test <- confusionMatrix(factor(Test.pred),
factor(as.matrix(Test.resp)),
positive="Yes")
TP = CM.test$table[2,2] # True Pos
TN = CM.test$table[1,1] # True Neg
FP = CM.test$table[2,1] # False Pos
FN = CM.test$table[1,2] # False Neg
cost[i] = sum(c(TP, TN, FP, FN) * cost.list)
}
plot(threshold.list, cost, xlab="threshold")
cost.list
which.min(cost)
min(cost)
threshold.list[which.min(cost)]