1. Load Data and Prepare

(top)

library(tidyverse)
Heart <- read_csv(file="https://nmimoto.github.io/datasets/heart.csv")

# Column manipulation
Heart2 <- Heart %>%
    select(-"index")  %>%            # remove col named "index"
    rename(resp=AHD) %>%             # rename "AHD" as "resp"
    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
table(Heart2$resp)    # note that first level is "No"
## 
##  No Yes 
## 164 139
# Low p-value from Chi-sq test of association
#  3 4 8 10 11 12 13 14
#  "Sex" "ChestPain" "RestECG" "ExAng" "Oldpeak" "Slope" "Ca" "Thal"
# High p-value from Chi-sq test of association
#  7
#  "Fbs"
# need to remove rows with NA.
Orig <- Heart2

#- Check for N/A in data. Remove if there's any.
  summary(Orig)
  sum(is.na(Orig))
  # If there is na in the data, run below
  dim(Orig)
  Orig <- Orig %>% na.omit()
  dim(Orig)


2. Data Separation

(top)

# 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=2452
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. Demo of How Decision Tree has high variance

(top)

library(tree)          # install.packages('tree', repos='https://cran.case.edu/')
tree31 <- tree(resp ~., CV.train[[1]])
summary(tree31)
## 
## Classification tree:
## tree(formula = resp ~ ., data = CV.train[[1]])
## Variables actually used in tree construction:
## [1] "ChestPain" "Age"       "Sex"       "Chol"      "Ca"        "Thal"     
## [7] "ExAng"     "MaxHR"     "RestBP"   
## Number of terminal nodes:  17 
## Residual mean deviance:  0.4116 = 73.67 / 179 
## Misclassification error rate: 0.09694 = 19 / 196
plot(tree31)
text(tree31, pretty=1, cex=.7)

#---------
tree32 <- tree(resp ~., CV.train[[2]])
summary(tree32)
## 
## Classification tree:
## tree(formula = resp ~ ., data = CV.train[[2]])
## Variables actually used in tree construction:
## [1] "Thal"      "Age"       "ChestPain" "Oldpeak"   "RestBP"    "Ca"       
## [7] "Chol"      "MaxHR"    
## Number of terminal nodes:  16 
## Residual mean deviance:  0.5564 = 100.1 / 180 
## Misclassification error rate: 0.1276 = 25 / 196
plot(tree32)
text(tree32, pretty=0, cex=.7)


4. Decision Tree Growing using rpart(). It has color plotting.

(top)

### (Train set fitting only)

library(rpart)         # install.packages('rpart', repos='https://cran.case.edu/')
library(rpart.plot)    # install.packages('rpart.plot', repos='https://cran.case.edu/')
tree00<- rpart(resp~., data=Train.set)
summary(tree00)

rpart.plot(tree00)     # now you can plot in color

plot(tree00)           # old way of plotting
text(tree00)


5. Decision Tree (Grow and Prune using tree())

(top)

### a. Growing the tree
library(tree)
tree1 = tree(resp~., Train.set)
summary(tree1)
## 
## Classification tree:
## tree(formula = resp ~ ., data = Train.set)
## Variables actually used in tree construction:
##  [1] "Thal"      "Ca"        "Age"       "MaxHR"     "RestBP"    "ChestPain"
##  [7] "Chol"      "Sex"       "Oldpeak"   "RestECG"  
## Number of terminal nodes:  17 
## Residual mean deviance:  0.5516 = 125.8 / 228 
## Misclassification error rate: 0.1347 = 33 / 245
tree1
## node), split, n, deviance, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 245 339.400 No ( 0.51429 0.48571 )  
##    2) Thal: normal 134 151.800 No ( 0.74627 0.25373 )  
##      4) Ca < 0.5 91  70.970 No ( 0.86813 0.13187 )  
##        8) Age < 54.5 55  17.180 No ( 0.96364 0.03636 )  
##         16) MaxHR < 152.5 14  11.480 No ( 0.85714 0.14286 ) *
##         17) MaxHR > 152.5 41   0.000 No ( 1.00000 0.00000 ) *
##        9) Age > 54.5 36  42.540 No ( 0.72222 0.27778 )  
##         18) RestBP < 157 30  30.020 No ( 0.80000 0.20000 ) *
##         19) RestBP > 157 6   7.638 Yes ( 0.33333 0.66667 ) *
##      5) Ca > 0.5 43  59.590 Yes ( 0.48837 0.51163 )  
##       10) ChestPain: nonanginal,nontypical,typical 24  26.990 No ( 0.75000 0.25000 )  
##         20) Chol < 236.5 7   9.561 Yes ( 0.42857 0.57143 ) *
##         21) Chol > 236.5 17  12.320 No ( 0.88235 0.11765 )  
##           42) ChestPain: nonanginal 11   0.000 No ( 1.00000 0.00000 ) *
##           43) ChestPain: nontypical,typical 6   7.638 No ( 0.66667 0.33333 ) *
##       11) ChestPain: asymptomatic 19  16.570 Yes ( 0.15789 0.84211 )  
##         22) Sex: 0 6   8.318 No ( 0.50000 0.50000 ) *
##         23) Sex: 1 13   0.000 Yes ( 0.00000 1.00000 ) *
##    3) Thal: fixed,reversable 111 120.800 Yes ( 0.23423 0.76577 )  
##      6) Oldpeak < 0.55 30  41.050 No ( 0.56667 0.43333 )  
##       12) Chol < 237.5 19  21.900 No ( 0.73684 0.26316 )  
##         24) MaxHR < 149 5   0.000 No ( 1.00000 0.00000 ) *
##         25) MaxHR > 149 14  18.250 No ( 0.64286 0.35714 ) *
##       13) Chol > 237.5 11  12.890 Yes ( 0.27273 0.72727 )  
##         26) RestECG: 0 6   8.318 No ( 0.50000 0.50000 ) *
##         27) RestECG: 2 5   0.000 Yes ( 0.00000 1.00000 ) *
##      7) Oldpeak > 0.55 81  56.510 Yes ( 0.11111 0.88889 )  
##       14) ChestPain: nonanginal,typical 21  27.910 Yes ( 0.38095 0.61905 )  
##         28) MaxHR < 143.5 7   0.000 Yes ( 0.00000 1.00000 ) *
##         29) MaxHR > 143.5 14  19.120 No ( 0.57143 0.42857 ) *
##       15) ChestPain: asymptomatic,nontypical 60  10.170 Yes ( 0.01667 0.98333 )  
##         30) RestBP < 112 6   5.407 Yes ( 0.16667 0.83333 ) *
##         31) RestBP > 112 54   0.000 Yes ( 0.00000 1.00000 ) *
plot(tree1)
text(tree1, pretty=0, cex=1)

# Check the training fit and test prediction
#Train.fitted = predict(tree1, type="vector")
#Test.pred    = predict(tree1, Test.set, type="vector")

threshold = .5

Train.prob = predict(tree1, type="vector")[,"Yes"]
Train.pred = ifelse(Train.prob > threshold, "Yes", "No")
Test.prob  = predict(tree1, Test.set, type="vector")[,"Yes"]
Test.pred  = ifelse(Test.prob > threshold, "Yes", "No")

    ###----------
    ### a. Pruning the tree
    set.seed(my.seed)
    cv.for.pruning = cv.tree(tree1, FUN=prune.misclass, K=5)   #automatic 5-fold CV
         #use FUN= prune.tree if you are doing regression tree
    names(cv.for.pruning)
## [1] "size"   "dev"    "k"      "method"
    plot(cv.for.pruning$size, cv.for.pruning$dev, type="b")

    plot(cv.for.pruning$k,    cv.for.pruning$dev, type="b")

                                            # size is the number of terminal nodes
                                            # dev is the Av. CV error rate
                                            # k is the pruning parameter (alpha)
    cv.for.pruning
## $size
## [1] 17 11  6  4  2  1
## 
## $dev
## [1]  69  67  66  71  77 119
## 
## $k
## [1] -Inf  0.0  1.0  4.5  6.5 59.0
## 
## $method
## [1] "misclass"
## 
## attr(,"class")
## [1] "prune"         "tree.sequence"
    # We conclude that it's better to purne at 6 nodes.
    pruned1 = prune.tree(tree1, best=6)  # prune to 6-node tree
    plot(pruned1)
    text(pruned1, pretty=0, cex=1)

    ###----------
    ### b. Pick a threshold and visualize fit
    Chosen.model <- pruned1  # tree1 (unpruned) or pruned1
    threshold = .6   # pick a threshold

    #- Check the training set accuracy
    library(caret)     # install.packages('caret', repos='https://cran.case.edu/')
    Train.prob = predict(Chosen.model, type="vector")[,"Yes"]
    Train.pred = ifelse(Train.prob > threshold, "Yes", "No")  # Turn the fitted values to Up/Down using the threshold
    Test.prob  = predict(Chosen.model, Test.set, type="vector")[,"Yes"]
    Test.pred  = ifelse(Test.prob > threshold, "Yes", "No")
    CM.train <- caret::confusionMatrix(factor(Train.pred), factor(as.matrix(Train.resp)), positive="Yes")
    CM.test <- caret::confusionMatrix(factor(Test.pred), factor(as.matrix(Test.resp)), positive="Yes")

    CM.train            # Training set result
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  No Yes
##        No  114  31
##        Yes  12  88
##                                        
##                Accuracy : 0.8245       
##                  95% CI : (0.771, 0.87)
##     No Information Rate : 0.5143       
##     P-Value [Acc > NIR] : < 2.2e-16    
##                                        
##                   Kappa : 0.6471       
##                                        
##  Mcnemar's Test P-Value : 0.006052     
##                                        
##             Sensitivity : 0.7395       
##             Specificity : 0.9048       
##          Pos Pred Value : 0.8800       
##          Neg Pred Value : 0.7862       
##              Prevalence : 0.4857       
##          Detection Rate : 0.3592       
##    Detection Prevalence : 0.4082       
##       Balanced Accuracy : 0.8221       
##                                        
##        'Positive' Class : Yes          
## 
    CM.train$table      # output just the table
##           Reference
## Prediction  No Yes
##        No  114  31
##        Yes  12  88
    CM.train[["byClass"]][["Sensitivity"]]
## [1] 0.7394958
    CM.train[["byClass"]][["Specificity"]]
## [1] 0.9047619
    CM.test             # Testing set
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction No Yes
##        No  29   7
##        Yes  5  11
##                                           
##                Accuracy : 0.7692          
##                  95% CI : (0.6316, 0.8747)
##     No Information Rate : 0.6538          
##     P-Value [Acc > NIR] : 0.05117         
##                                           
##                   Kappa : 0.4765          
##                                           
##  Mcnemar's Test P-Value : 0.77283         
##                                           
##             Sensitivity : 0.6111          
##             Specificity : 0.8529          
##          Pos Pred Value : 0.6875          
##          Neg Pred Value : 0.8056          
##              Prevalence : 0.3462          
##          Detection Rate : 0.2115          
##    Detection Prevalence : 0.3077          
##       Balanced Accuracy : 0.7320          
##                                           
##        'Positive' Class : Yes             
## 
    CM.test$table      # output just the table
##           Reference
## Prediction No Yes
##        No  29   7
##        Yes  5  11
    colSums(CM.test$table) / sum(colSums(CM.test$table))    # % of Actual Yes/No
##        No       Yes 
## 0.6538462 0.3461538
    rowSums(CM.test$table) / sum(rowSums(CM.test$table))    # % of predicted Yes/No
##        No       Yes 
## 0.6923077 0.3076923
    ###----------
    ### c.  Output ROC curve and AUC for all threshold
    library(pROC)
    #- Training Set
    plot.roc(factor(as.matrix(Train.resp)),  Train.prob, levels=c("No", "Yes"))
    # point corresponding to CM.train
    abline(h=CM.train[["byClass"]][["Sensitivity"]], v=CM.train[["byClass"]][["Specificity"]], col="red")
    auc.train = auc(factor(as.matrix(Train.resp)), Train.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Train AUC=",round(auc.train, 3)))

    #- Test Set
    plot.roc(factor(as.matrix(Test.resp)),  Test.prob, levels=c("No", "Yes"))
    # point corresponding to CM.test
    abline(h=CM.test[["byClass"]][["Sensitivity"]], v=CM.test[["byClass"]][["Specificity"]], col="red")
    auc.test = auc(factor(as.matrix(Test.resp)), Test.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Test AUC=",round(auc.test, 3)))

    GrowAndPrune = c(auc.train, auc.test)
    GrowAndPrune
## [1] 0.8862212 0.8856209


6. Bagging (Train and Test only)

(top)

library(ipred)    #  install.packages('ipred', repos='https://cran.case.edu/')
set.seed(my.seed)
treeBag01 = bagging(resp~., data=Train.set, nbagg=500,
                    coob = TRUE,
                    control = rpart.control(minsplit = 2, cp = 0))
treeBag01
## 
## Bagging classification trees with 500 bootstrap replications 
## 
## Call: bagging.data.frame(formula = resp ~ ., data = Train.set, nbagg = 500, 
##     coob = TRUE, control = rpart.control(minsplit = 2, cp = 0))
## 
## Out-of-bag estimate of misclassification error:  0.1918
###----------
Chosen.model <- treeBag01

    ###----------
    ### a. Pick a threshold and visualize fit
    threshold = .5    # pick a threshold

    #- Check the training set accuracy
    library(caret)     # install.packages('caret', repos='https://cran.case.edu/')
    Train.prob = predict(Chosen.model, type="prob")[,"Yes"]
    Train.pred = ifelse(Train.prob > threshold, "Yes", "No")  # Turn the fitted values to Up/Down using the threshold
    Test.prob = predict(Chosen.model, newdata=Test.set, type="prob")[,"Yes"]
    Test.pred = ifelse(Test.prob > threshold, "Yes", "No")  # Turn the fitted values to Up/Down using the threshold
    CM.train <- caret::confusionMatrix(factor(Train.pred), factor(as.matrix(Train.resp)), positive="Yes")
    CM.test  <- caret::confusionMatrix(factor(Test.pred),  factor(as.matrix(Test.resp)), positive="Yes")

    CM.train            # Training set result
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  No Yes
##        No  104  25
##        Yes  22  94
##                                           
##                Accuracy : 0.8082          
##                  95% CI : (0.7532, 0.8555)
##     No Information Rate : 0.5143          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.6157          
##                                           
##  Mcnemar's Test P-Value : 0.7705          
##                                           
##             Sensitivity : 0.7899          
##             Specificity : 0.8254          
##          Pos Pred Value : 0.8103          
##          Neg Pred Value : 0.8062          
##              Prevalence : 0.4857          
##          Detection Rate : 0.3837          
##    Detection Prevalence : 0.4735          
##       Balanced Accuracy : 0.8077          
##                                           
##        'Positive' Class : Yes             
## 
    CM.train$table      # output just the table
##           Reference
## Prediction  No Yes
##        No  104  25
##        Yes  22  94
    CM.train[["byClass"]][["Sensitivity"]]
## [1] 0.789916
    CM.train[["byClass"]][["Specificity"]]
## [1] 0.8253968
    CM.test             # Testing set
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction No Yes
##        No  29   5
##        Yes  5  13
##                                           
##                Accuracy : 0.8077          
##                  95% CI : (0.6747, 0.9037)
##     No Information Rate : 0.6538          
##     P-Value [Acc > NIR] : 0.01172         
##                                           
##                   Kappa : 0.5752          
##                                           
##  Mcnemar's Test P-Value : 1.00000         
##                                           
##             Sensitivity : 0.7222          
##             Specificity : 0.8529          
##          Pos Pred Value : 0.7222          
##          Neg Pred Value : 0.8529          
##              Prevalence : 0.3462          
##          Detection Rate : 0.2500          
##    Detection Prevalence : 0.3462          
##       Balanced Accuracy : 0.7876          
##                                           
##        'Positive' Class : Yes             
## 
    CM.test$table      # output just the table
##           Reference
## Prediction No Yes
##        No  29   5
##        Yes  5  13
    colSums(CM.test$table) / sum(colSums(CM.test$table))    # % of Actual Yes/No
##        No       Yes 
## 0.6538462 0.3461538
    rowSums(CM.test$table) / sum(rowSums(CM.test$table))    # % of predicted Yes/No
##        No       Yes 
## 0.6538462 0.3461538
    ###----------
    ### b. Output ROC curve and AUC for all threshold
    library(pROC)
    #- Training Set
    plot.roc(factor(as.matrix(Train.resp)),  Train.prob, levels=c("No", "Yes"))
    # point corresponding to CM.train
    abline(h=CM.train[["byClass"]][["Sensitivity"]], v=CM.train[["byClass"]][["Specificity"]], col="red")
    auc.train = auc(factor(as.matrix(Train.resp)), Train.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Train AUC=",round(auc.train, 3)))

    #- Test Set
    plot.roc(factor(as.matrix(Test.resp)),  Test.prob, levels=c("No", "Yes"))
    # point corresponding to CM.test
    abline(h=CM.test[["byClass"]][["Sensitivity"]], v=CM.test[["byClass"]][["Specificity"]], col="red")
    auc.test = auc(factor(as.matrix(Test.resp)), Test.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Test AUC=",round(auc.test, 3)))

    Bagging = c(auc.train, auc.test)
    Bagging
## [1] 0.8768841 0.9240196
    # Training ROC and Test ROC side by side
    layout(matrix(1:2, 1, 2))
    plot.roc(factor(as.matrix(Train.resp)),  Train.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Train AUC=",round(auc.train, 3)))
    plot.roc(factor(as.matrix(Test.resp)),  Test.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Test AUC=",round(auc.test, 3)))

    layout(1)

    Bagging
## [1] 0.8768841 0.9240196


7. Random Forest (Train and Test only)

(top)

library(randomForest)
set.seed(my.seed)
treeBag01 = randomForest(resp~., data=Train.set, mtry=13, ntree=500, importance=TRUE)
treeBag01
## 
## Call:
##  randomForest(formula = resp ~ ., data = Train.set, mtry = 13,      ntree = 500, importance = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 13
## 
##         OOB estimate of  error rate: 20.82%
## Confusion matrix:
##      No Yes class.error
## No  103  23   0.1825397
## Yes  28  91   0.2352941
importance (treeBag01)
##                   No         Yes MeanDecreaseAccuracy MeanDecreaseGini
## Age        3.2933590  4.47584534            5.3289205        9.6417137
## Sex        7.9477240  3.62846731            8.3540488        2.8661468
## ChestPain 11.6628905 16.01873751           18.6622589       17.1371934
## RestBP     3.5721179  2.67930078            4.1059576        9.1056263
## Chol       1.4974958 -5.00455871           -2.3398021       10.7054969
## Fbs        0.8137165 -0.46957564            0.1068307        0.2788583
## RestECG    3.5313202 -1.31657393            1.2208067        1.3675410
## MaxHR      3.3809587 -0.01992428            2.7277163        9.5447071
## ExAng      5.8807054  3.20243753            6.6196601        4.7205487
## Oldpeak   17.6547766 13.74079052           21.7354998       15.8952784
## Slope      3.2229843  8.51488644            9.1738665        4.5061035
## Ca        19.4679635 14.15465801           22.5078493       15.1270350
## Thal      19.5375340 14.42202065           22.9233971       20.9736938
varImpPlot (treeBag01)

###----------
Chosen.model <- treeBag01

    ###----------
    ### a. Pick a threshold and visualize fit
    threshold = .5    # pick a threshold

    #- Check the training set accuracy
    library(caret)     # install.packages('caret', repos='https://cran.case.edu/')
    Train.prob = predict(Chosen.model, type="prob")[,"Yes"]
    Train.pred = ifelse(Train.prob > threshold, "Yes", "No")  # Turn the fitted values to Up/Down using the threshold
    Test.prob = predict(Chosen.model, newdata=Test.set, type="prob")[,"Yes"]
    Test.pred = ifelse(Test.prob > threshold, "Yes", "No")  # Turn the fitted values to Up/Down using the threshold
    CM.train <- caret::confusionMatrix(factor(Train.pred), factor(as.matrix(Train.resp)), positive="Yes")
    CM.test  <- caret::confusionMatrix(factor(Test.pred),  factor(as.matrix(Test.resp)), positive="Yes")

    CM.train            # Training set result
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  No Yes
##        No  103  28
##        Yes  23  91
##                                           
##                Accuracy : 0.7918          
##                  95% CI : (0.7356, 0.8409)
##     No Information Rate : 0.5143          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.5828          
##                                           
##  Mcnemar's Test P-Value : 0.5754          
##                                           
##             Sensitivity : 0.7647          
##             Specificity : 0.8175          
##          Pos Pred Value : 0.7982          
##          Neg Pred Value : 0.7863          
##              Prevalence : 0.4857          
##          Detection Rate : 0.3714          
##    Detection Prevalence : 0.4653          
##       Balanced Accuracy : 0.7911          
##                                           
##        'Positive' Class : Yes             
## 
    CM.train$table      # output just the table
##           Reference
## Prediction  No Yes
##        No  103  28
##        Yes  23  91
    CM.train[["byClass"]][["Sensitivity"]]
## [1] 0.7647059
    CM.train[["byClass"]][["Specificity"]]
## [1] 0.8174603
    CM.test             # Testing set
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction No Yes
##        No  28   5
##        Yes  6  13
##                                          
##                Accuracy : 0.7885         
##                  95% CI : (0.653, 0.8894)
##     No Information Rate : 0.6538         
##     P-Value [Acc > NIR] : 0.02577        
##                                          
##                   Kappa : 0.5387         
##                                          
##  Mcnemar's Test P-Value : 1.00000        
##                                          
##             Sensitivity : 0.7222         
##             Specificity : 0.8235         
##          Pos Pred Value : 0.6842         
##          Neg Pred Value : 0.8485         
##              Prevalence : 0.3462         
##          Detection Rate : 0.2500         
##    Detection Prevalence : 0.3654         
##       Balanced Accuracy : 0.7729         
##                                          
##        'Positive' Class : Yes            
## 
    CM.test$table      # output just the table
##           Reference
## Prediction No Yes
##        No  28   5
##        Yes  6  13
    colSums(CM.test$table) / sum(colSums(CM.test$table))    # % of Actual Yes/No
##        No       Yes 
## 0.6538462 0.3461538
    rowSums(CM.test$table) / sum(rowSums(CM.test$table))    # % of predicted Yes/No
##        No       Yes 
## 0.6346154 0.3653846
    ###----------
    ### b. Output ROC curve and AUC for all threshold
    library(pROC)
    #- Training Set
    plot.roc(factor(as.matrix(Train.resp)),  Train.prob, levels=c("No", "Yes"))
    # point corresponding to CM.train
    abline(h=CM.train[["byClass"]][["Sensitivity"]], v=CM.train[["byClass"]][["Specificity"]], col="red")
    auc.train = auc(factor(as.matrix(Train.resp)), Train.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Train AUC=",round(auc.train, 3)))

    #- Test Set
    plot.roc(factor(as.matrix(Test.resp)),  Test.prob, levels=c("No", "Yes"))
    # point corresponding to CM.test
    abline(h=CM.test[["byClass"]][["Sensitivity"]], v=CM.test[["byClass"]][["Specificity"]], col="red")
    auc.test = auc(factor(as.matrix(Test.resp)), Test.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Test AUC=",round(auc.test, 3)))

    RandomForest = c(auc.train, auc.test)
    RandomForest
## [1] 0.8644124 0.9019608
    # Training ROC and Test ROC side by side
    layout(matrix(1:2, 1, 2))
    plot.roc(factor(as.matrix(Train.resp)),  Train.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Train AUC=",round(auc.train, 3)))
    plot.roc(factor(as.matrix(Test.resp)),  Test.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Test AUC=",round(auc.test, 3)))

    layout(1)

    RandomForest
## [1] 0.8644124 0.9019608


8. Boosting (Train and Test only)

(top)

library(gbm)     # install.packages('gbm', repos='https://cran.case.edu/')
set.seed(my.seed)
treeBT01 = gbm::gbm(as.numeric(resp=="Yes") ~., data=Train.set,
                    distribution="bernoulli",   # for binary classification
                    n.trees=5000,
                    interaction.depth=4)

summary(treeBT01)

##                 var   rel.inf
## ChestPain ChestPain 13.480058
## Thal           Thal 12.904968
## Age             Age 11.883537
## Oldpeak     Oldpeak 11.838874
## MaxHR         MaxHR 11.199102
## Chol           Chol 10.427892
## Ca               Ca 10.159034
## RestBP       RestBP  7.376253
## ExAng         ExAng  3.632848
## Sex             Sex  3.319270
## Slope         Slope  1.840794
## RestECG     RestECG  1.784650
## Fbs             Fbs  0.152720
plot(treeBT01, i="Chol")

plot(treeBT01, i="MaxHR")

###----------
Chosen.model <- treeBT01

    ###----------
    ### a. Pick a threshold and visualize fit
    threshold = .5    # pick a threshold

    #- Check the training set accuracy
    library(caret)     # install.packages('caret', repos='https://cran.case.edu/')
    Train.prob = predict(Chosen.model, type="response", n.trees=1000)
    Train.pred = ifelse(Train.prob > threshold, "Yes", "No")  # Turn the fitted values to Up/Down using the threshold
    Test.prob = predict(Chosen.model, newdata=Test.set, type="response", n.trees=1000)
    Test.pred = ifelse(Test.prob > threshold, "Yes", "No")  # Turn the fitted values to Up/Down using the threshold
    CM.train <- caret::confusionMatrix(factor(Train.pred), factor(as.matrix(Train.resp)), positive="Yes")
    CM.test  <- caret::confusionMatrix(factor(Test.pred),  factor(as.matrix(Test.resp)), positive="Yes")

    CM.train            # Training set result
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  No Yes
##        No  126   0
##        Yes   0 119
##                                      
##                Accuracy : 1          
##                  95% CI : (0.9851, 1)
##     No Information Rate : 0.5143     
##     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.4857     
##          Detection Rate : 0.4857     
##    Detection Prevalence : 0.4857     
##       Balanced Accuracy : 1.0000     
##                                      
##        'Positive' Class : Yes        
## 
    CM.train$table      # output just the table
##           Reference
## Prediction  No Yes
##        No  126   0
##        Yes   0 119
    CM.train[["byClass"]][["Sensitivity"]]
## [1] 1
    CM.train[["byClass"]][["Specificity"]]
## [1] 1
    CM.test             # Testing set
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction No Yes
##        No  30   7
##        Yes  4  11
##                                          
##                Accuracy : 0.7885         
##                  95% CI : (0.653, 0.8894)
##     No Information Rate : 0.6538         
##     P-Value [Acc > NIR] : 0.02577        
##                                          
##                   Kappa : 0.5136         
##                                          
##  Mcnemar's Test P-Value : 0.54649        
##                                          
##             Sensitivity : 0.6111         
##             Specificity : 0.8824         
##          Pos Pred Value : 0.7333         
##          Neg Pred Value : 0.8108         
##              Prevalence : 0.3462         
##          Detection Rate : 0.2115         
##    Detection Prevalence : 0.2885         
##       Balanced Accuracy : 0.7467         
##                                          
##        'Positive' Class : Yes            
## 
    CM.test$table      # output just the table
##           Reference
## Prediction No Yes
##        No  30   7
##        Yes  4  11
    colSums(CM.test$table) / sum(colSums(CM.test$table))    # % of Actual Yes/No
##        No       Yes 
## 0.6538462 0.3461538
    rowSums(CM.test$table) / sum(rowSums(CM.test$table))    # % of predicted Yes/No
##        No       Yes 
## 0.7115385 0.2884615
    ###----------
    ### b. Output ROC curve and AUC for all threshold
    library(pROC)
    #- Training Set
    plot.roc(factor(as.matrix(Train.resp)),  Train.prob, levels=c("No", "Yes"))
    # point corresponding to CM.train
    abline(h=CM.train[["byClass"]][["Sensitivity"]], v=CM.train[["byClass"]][["Specificity"]], col="red")
    auc.train = auc(factor(as.matrix(Train.resp)), Train.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Train AUC=",round(auc.train, 3)))

    #- Test Set
    plot.roc(factor(as.matrix(Test.resp)),  Test.prob, levels=c("No", "Yes"))
    # point corresponding to CM.test
    abline(h=CM.test[["byClass"]][["Sensitivity"]], v=CM.test[["byClass"]][["Specificity"]], col="red")
    auc.test = auc(factor(as.matrix(Test.resp)), Test.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Test AUC=",round(auc.test, 3)))

    Boosting = c(auc.train, auc.test)
    Boosting
## [1] 1.0000000 0.9166667
    # Training ROC and Test ROC side by side
    layout(matrix(1:2, 1, 2))
    plot.roc(factor(as.matrix(Train.resp)),  Train.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Train AUC=",round(auc.train, 3)))
    plot.roc(factor(as.matrix(Test.resp)),  Test.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Test AUC=",round(auc.test, 3)))

    layout(1)

    Boosting
## [1] 1.0000000 0.9166667


9. Logistic Regression for comparison (Train vs Test only)

(top)

Fit03 <- glm(resp ~ .,  family=binomial, data=Train.set)
summary(Fit03)
## 
## Call:
## glm(formula = resp ~ ., family = binomial, data = Train.set)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -4.488742   3.118359  -1.439  0.15002    
## Age                 -0.013622   0.026098  -0.522  0.60169    
## Sex1                 1.450860   0.572178   2.536  0.01122 *  
## ChestPainnonanginal -1.577396   0.528432  -2.985  0.00284 ** 
## ChestPainnontypical -1.228959   0.661458  -1.858  0.06318 .  
## ChestPaintypical    -1.750393   0.667910  -2.621  0.00877 ** 
## RestBP               0.021695   0.012304   1.763  0.07785 .  
## Chol                 0.007615   0.004872   1.563  0.11803    
## Fbs1                -0.432408   0.641247  -0.674  0.50011    
## RestECG1             1.077491   2.439620   0.442  0.65873    
## RestECG2             0.494664   0.415305   1.191  0.23362    
## MaxHR               -0.015378   0.012298  -1.250  0.21114    
## ExAng1               0.884633   0.478870   1.847  0.06470 .  
## Oldpeak              0.392126   0.244580   1.603  0.10888    
## Slope                0.511958   0.410723   1.246  0.21259    
## Ca                   1.087811   0.271716   4.003 6.24e-05 ***
## Thalnormal          -0.415896   0.863582  -0.482  0.63009    
## Thalreversable       0.817687   0.837796   0.976  0.32907    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 339.44  on 244  degrees of freedom
## Residual deviance: 168.82  on 227  degrees of freedom
## AIC: 204.82
## 
## Number of Fisher Scoring iterations: 6
###----------
Chosen.model <- Fit03

    ###----------
    #- a Output result for given threshold value
    threshold = .5 # pick a threshold

    #- Check the training set accuracy
    library(caret)
    Train.prob =predict(Chosen.model, type ="response")
    Train.pred = ifelse(Train.prob > threshold, "Yes", "No")  # Turn the fitted values to Up/Down using threshold of .5
    Test.prob = predict(Chosen.model, newdata=Test.set, type="response")
    Test.pred  = ifelse(Test.prob  > threshold, "Yes", "No")
    CM.train <- confusionMatrix(factor(Train.pred), factor(as.matrix(Train.resp)), positive="Yes")
    CM.test  <- confusionMatrix(factor(Test.pred),  factor(as.matrix(Test.resp)),  positive="Yes")

    CM.train            # Training set result
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  No Yes
##        No  112  23
##        Yes  14  96
##                                           
##                Accuracy : 0.849           
##                  95% CI : (0.7979, 0.8914)
##     No Information Rate : 0.5143          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.6971          
##                                           
##  Mcnemar's Test P-Value : 0.1884          
##                                           
##             Sensitivity : 0.8067          
##             Specificity : 0.8889          
##          Pos Pred Value : 0.8727          
##          Neg Pred Value : 0.8296          
##              Prevalence : 0.4857          
##          Detection Rate : 0.3918          
##    Detection Prevalence : 0.4490          
##       Balanced Accuracy : 0.8478          
##                                           
##        'Positive' Class : Yes             
## 
    CM.train$table      # output just the table
##           Reference
## Prediction  No Yes
##        No  112  23
##        Yes  14  96
    CM.train[["byClass"]][["Sensitivity"]]
## [1] 0.8067227
    CM.train[["byClass"]][["Specificity"]]
## [1] 0.8888889
    CM.test             # Testing set
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction No Yes
##        No  32   3
##        Yes  2  15
##                                          
##                Accuracy : 0.9038         
##                  95% CI : (0.7897, 0.968)
##     No Information Rate : 0.6538         
##     P-Value [Acc > NIR] : 3.378e-05      
##                                          
##                   Kappa : 0.7848         
##                                          
##  Mcnemar's Test P-Value : 1              
##                                          
##             Sensitivity : 0.8333         
##             Specificity : 0.9412         
##          Pos Pred Value : 0.8824         
##          Neg Pred Value : 0.9143         
##              Prevalence : 0.3462         
##          Detection Rate : 0.2885         
##    Detection Prevalence : 0.3269         
##       Balanced Accuracy : 0.8873         
##                                          
##        'Positive' Class : Yes            
## 
    CM.test$table      # output just the table
##           Reference
## Prediction No Yes
##        No  32   3
##        Yes  2  15
    ###----------
    #- b Output ROC curve and AUC for all threshold
    library(pROC)
    #- Training Set
    plot.roc(factor(as.matrix(Train.resp)),  Train.prob, levels=c("No", "Yes"))
    # point corresponding to CM.train
    abline(h=CM.train[["byClass"]][["Sensitivity"]], v=CM.train[["byClass"]][["Specificity"]], col="red")
    auc.train = auc(factor(as.matrix(Train.resp)), Train.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Train AUC=",round(auc.train, 3)))

    #- Test Set
    plot.roc(factor(as.matrix(Test.resp)),  Test.prob, levels=c("No", "Yes"))
    # point corresponding to CM.test
    abline(h=CM.test[["byClass"]][["Sensitivity"]], v=CM.test[["byClass"]][["Specificity"]], col="red")
    auc.test = auc(factor(as.matrix(Test.resp)), Test.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Test AUC=",round(auc.test, 3)))

    LogiReg <- data.frame(
        tr.AUC   = auc.train,
        test.AUC = auc.test
    )
    LogiReg
##      tr.AUC  test.AUC
## 1 0.9264372 0.9509804
    layout(matrix(1:2, 1, 2))
    plot.roc(factor(as.matrix(Train.resp)),  Train.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Train AUC=",round(auc.train, 3)))
    plot.roc(factor(as.matrix(Test.resp)),  Test.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Test AUC=",round(auc.test, 3)))

    layout(1)

    LogiReg
##      tr.AUC  test.AUC
## 1 0.9264372 0.9509804


Compare Performances

(top)

##  With my.seed=7211
##  Compare      (auc.train, auc.test)
##  LogiReg       0.9437945  0.8964497
##  GrowAndPrune  0.8898077  0.8853550
##  Bagging       0.8768657  0.9215976
##  RandomForest  0.8659742  0.9260355
##  Boosting      1.0000000  0.8727811


## Compare (auc.train, auc.test)
LogiReg
##      tr.AUC  test.AUC
## 1 0.9264372 0.9509804
GrowAndPrune
## [1] 0.8862212 0.8856209
Bagging
## [1] 0.8768841 0.9240196
RandomForest
## [1] 0.8644124 0.9019608
Boosting
## [1] 1.0000000 0.9166667


10. Threshold Picker

(top)

#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(0.01,.99,.01)    # grid for threshold
cost=0
library(caret)      # for confusionMatrix

#  Have to be switched depending on the model
Chosen.model <- pruned1;    Test.prob  = predict(Chosen.model, Test.set, type="vector")[,"Yes"]
# Chosen.model <- treeBag01;  Test.prob = predict(Chosen.model, newdata=Test.set, type="prob")[,"Yes"]
# Chosen.model <- treeBT01;   Test.prob = predict(Chosen.model, newdata=Test.set, type="response", n.trees=1000)


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)
}
## Warning in confusionMatrix.default(factor(Test.pred),
## factor(as.matrix(Test.resp)), : Levels are not in the same order for reference
## and data. Refactoring data to match.
## Warning in confusionMatrix.default(factor(Test.pred),
## factor(as.matrix(Test.resp)), : Levels are not in the same order for reference
## and data. Refactoring data to match.
## Warning in confusionMatrix.default(factor(Test.pred),
## factor(as.matrix(Test.resp)), : Levels are not in the same order for reference
## and data. Refactoring data to match.
## Warning in confusionMatrix.default(factor(Test.pred),
## factor(as.matrix(Test.resp)), : Levels are not in the same order for reference
## and data. Refactoring data to match.
## Warning in confusionMatrix.default(factor(Test.pred),
## factor(as.matrix(Test.resp)), : Levels are not in the same order for reference
## and data. Refactoring data to match.
## Warning in confusionMatrix.default(factor(Test.pred),
## factor(as.matrix(Test.resp)), : Levels are not in the same order for reference
## and data. Refactoring data to match.
## Warning in confusionMatrix.default(factor(Test.pred),
## factor(as.matrix(Test.resp)), : Levels are not in the same order for reference
## and data. Refactoring data to match.
## Warning in confusionMatrix.default(factor(Test.pred),
## factor(as.matrix(Test.resp)), : Levels are not in the same order for reference
## and data. Refactoring data to match.
## Warning in confusionMatrix.default(factor(Test.pred),
## factor(as.matrix(Test.resp)), : Levels are not in the same order for reference
## and data. Refactoring data to match.
## Warning in confusionMatrix.default(factor(Test.pred),
## factor(as.matrix(Test.resp)), : Levels are not in the same order for reference
## and data. Refactoring data to match.
## Warning in confusionMatrix.default(factor(Test.pred),
## factor(as.matrix(Test.resp)), : Levels are not in the same order for reference
## and data. Refactoring data to match.
## Warning in confusionMatrix.default(factor(Test.pred),
## factor(as.matrix(Test.resp)), : Levels are not in the same order for reference
## and data. Refactoring data to match.
## Warning in confusionMatrix.default(factor(Test.pred),
## factor(as.matrix(Test.resp)), : Levels are not in the same order for reference
## and data. Refactoring data to match.
## Warning in confusionMatrix.default(factor(Test.pred),
## factor(as.matrix(Test.resp)), : Levels are not in the same order for reference
## and data. Refactoring data to match.
plot(threshold.list, cost, xlab="threshold")

cost.list
## [1] 0.0000000 0.0000000 0.6666667 0.3333333
which.min(cost)
## [1] 62
min(cost)
## [1] 3.333333
threshold.list[which.min(cost)]
## [1] 0.62


CODE ONLY

(top)

##--------------------------------------
## 1. Load Data and Prepare

library(tidyverse)
Heart <- read_csv(file="https://nmimoto.github.io/datasets/heart.csv")

# Column manipulation
Heart2 <- Heart %>%
    select(-"index")  %>%            # remove col named "index"
    rename(resp=AHD) %>%             # rename "AHD" as "resp"
    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"

# Low p-value from Chi-sq test of association
#  3 4 8 10 11 12 13 14
#  "Sex" "ChestPain" "RestECG" "ExAng" "Oldpeak" "Slope" "Ca" "Thal"
# High p-value from Chi-sq test of association
#  7
#  "Fbs"


# need to remove rows with NA.
Orig <- Heart2

#- Check for N/A in data. Remove if there's any.
  summary(Orig)
  sum(is.na(Orig))
  # If there is na in the data, run below
  dim(Orig)
  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=2452
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. Demo of How Decision Tree has high variance

library(tree)          # install.packages('tree', repos='https://cran.case.edu/')
tree31 <- tree(resp ~., CV.train[[1]])
summary(tree31)

plot(tree31)
text(tree31, pretty=1, cex=.7)

#---------
tree32 <- tree(resp ~., CV.train[[2]])
summary(tree32)

plot(tree32)
text(tree32, pretty=0, cex=.7)


##----------------------------------------
## 4. Decision Tree Growing using rpart().  It has color plotting.

### (Train set fitting only)

library(rpart)         # install.packages('rpart', repos='https://cran.case.edu/')
library(rpart.plot)    # install.packages('rpart.plot', repos='https://cran.case.edu/')
tree00<- rpart(resp~., data=Train.set)
summary(tree00)

rpart.plot(tree00)     # now you can plot in color

plot(tree00)           # old way of plotting
text(tree00)


##----------------------------------------
## 5. Decision Tree (Grow and Prune using tree())

### a. Growing the tree
library(tree)
tree1 = tree(resp~., Train.set)
summary(tree1)
tree1

plot(tree1)
text(tree1, pretty=0, cex=1)

# Check the training fit and test prediction
#Train.fitted = predict(tree1, type="vector")
#Test.pred    = predict(tree1, Test.set, type="vector")

threshold = .5

Train.prob = predict(tree1, type="vector")[,"Yes"]
Train.pred = ifelse(Train.prob > threshold, "Yes", "No")
Test.prob  = predict(tree1, Test.set, type="vector")[,"Yes"]
Test.pred  = ifelse(Test.prob > threshold, "Yes", "No")

    ###----------
    ### a. Pruning the tree
    set.seed(my.seed)
    cv.for.pruning = cv.tree(tree1, FUN=prune.misclass, K=5)   #automatic 5-fold CV
         #use FUN= prune.tree if you are doing regression tree
    names(cv.for.pruning)

    plot(cv.for.pruning$size, cv.for.pruning$dev, type="b")
    plot(cv.for.pruning$k,    cv.for.pruning$dev, type="b")
                                            # size is the number of terminal nodes
                                            # dev is the Av. CV error rate
                                            # k is the pruning parameter (alpha)
    cv.for.pruning

    # We conclude that it's better to purne at 6 nodes.
    pruned1 = prune.tree(tree1, best=6)  # prune to 6-node tree
    plot(pruned1)
    text(pruned1, pretty=0, cex=1)

    ###----------
    ### b. Pick a threshold and visualize fit
    Chosen.model <- pruned1  # tree1 (unpruned) or pruned1
    threshold = .6   # pick a threshold

    #- Check the training set accuracy
    library(caret)     # install.packages('caret', repos='https://cran.case.edu/')
    Train.prob = predict(Chosen.model, type="vector")[,"Yes"]
    Train.pred = ifelse(Train.prob > threshold, "Yes", "No")  # Turn the fitted values to Up/Down using the threshold
    Test.prob  = predict(Chosen.model, Test.set, type="vector")[,"Yes"]
    Test.pred  = ifelse(Test.prob > threshold, "Yes", "No")
    CM.train <- caret::confusionMatrix(factor(Train.pred), factor(as.matrix(Train.resp)), positive="Yes")
    CM.test <- caret::confusionMatrix(factor(Test.pred), factor(as.matrix(Test.resp)), positive="Yes")

    CM.train            # Training set result
    CM.train$table      # output just the table

    CM.train[["byClass"]][["Sensitivity"]]
    CM.train[["byClass"]][["Specificity"]]

    CM.test             # Testing set
    CM.test$table      # output just the table

    colSums(CM.test$table) / sum(colSums(CM.test$table))    # % of Actual Yes/No
    rowSums(CM.test$table) / sum(rowSums(CM.test$table))    # % of predicted Yes/No

    ###----------
    ### c.  Output ROC curve and AUC for all threshold
    library(pROC)
    #- Training Set
    plot.roc(factor(as.matrix(Train.resp)),  Train.prob, levels=c("No", "Yes"))
    # point corresponding to CM.train
    abline(h=CM.train[["byClass"]][["Sensitivity"]], v=CM.train[["byClass"]][["Specificity"]], col="red")
    auc.train = auc(factor(as.matrix(Train.resp)), Train.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Train AUC=",round(auc.train, 3)))

    #- Test Set
    plot.roc(factor(as.matrix(Test.resp)),  Test.prob, levels=c("No", "Yes"))
    # point corresponding to CM.test
    abline(h=CM.test[["byClass"]][["Sensitivity"]], v=CM.test[["byClass"]][["Specificity"]], col="red")
    auc.test = auc(factor(as.matrix(Test.resp)), Test.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Test AUC=",round(auc.test, 3)))

    GrowAndPrune = c(auc.train, auc.test)
    GrowAndPrune


##----------------------------------------
## 6. Bagging (Train and Test only)


library(ipred)    #  install.packages('ipred', repos='https://cran.case.edu/')
set.seed(my.seed)
treeBag01 = bagging(resp~., data=Train.set, nbagg=500,
                    coob = TRUE,
                    control = rpart.control(minsplit = 2, cp = 0))
treeBag01


###----------
Chosen.model <- treeBag01

    ###----------
    ### a. Pick a threshold and visualize fit
    threshold = .5    # pick a threshold

    #- Check the training set accuracy
    library(caret)     # install.packages('caret', repos='https://cran.case.edu/')
    Train.prob = predict(Chosen.model, type="prob")[,"Yes"]
    Train.pred = ifelse(Train.prob > threshold, "Yes", "No")  # Turn the fitted values to Up/Down using the threshold
    Test.prob = predict(Chosen.model, newdata=Test.set, type="prob")[,"Yes"]
    Test.pred = ifelse(Test.prob > threshold, "Yes", "No")  # Turn the fitted values to Up/Down using the threshold
    CM.train <- caret::confusionMatrix(factor(Train.pred), factor(as.matrix(Train.resp)), positive="Yes")
    CM.test  <- caret::confusionMatrix(factor(Test.pred),  factor(as.matrix(Test.resp)), positive="Yes")

    CM.train            # Training set result
    CM.train$table      # output just the table

    CM.train[["byClass"]][["Sensitivity"]]
    CM.train[["byClass"]][["Specificity"]]

    CM.test             # Testing set
    CM.test$table      # output just the table

    colSums(CM.test$table) / sum(colSums(CM.test$table))    # % of Actual Yes/No
    rowSums(CM.test$table) / sum(rowSums(CM.test$table))    # % of predicted Yes/No

    ###----------
    ### b. Output ROC curve and AUC for all threshold
    library(pROC)
    #- Training Set
    plot.roc(factor(as.matrix(Train.resp)),  Train.prob, levels=c("No", "Yes"))
    # point corresponding to CM.train
    abline(h=CM.train[["byClass"]][["Sensitivity"]], v=CM.train[["byClass"]][["Specificity"]], col="red")
    auc.train = auc(factor(as.matrix(Train.resp)), Train.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Train AUC=",round(auc.train, 3)))

    #- Test Set
    plot.roc(factor(as.matrix(Test.resp)),  Test.prob, levels=c("No", "Yes"))
    # point corresponding to CM.test
    abline(h=CM.test[["byClass"]][["Sensitivity"]], v=CM.test[["byClass"]][["Specificity"]], col="red")
    auc.test = auc(factor(as.matrix(Test.resp)), Test.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Test AUC=",round(auc.test, 3)))

    Bagging = c(auc.train, auc.test)
    Bagging

    # Training ROC and Test ROC side by side
    layout(matrix(1:2, 1, 2))
    plot.roc(factor(as.matrix(Train.resp)),  Train.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Train AUC=",round(auc.train, 3)))
    plot.roc(factor(as.matrix(Test.resp)),  Test.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Test AUC=",round(auc.test, 3)))
    layout(1)

    Bagging


##----------------------------------------
## 7. Random Forest (Train and Test only)

library(randomForest)
set.seed(my.seed)
treeBag01 = randomForest(resp~., data=Train.set, mtry=13, ntree=500, importance=TRUE)
treeBag01

importance (treeBag01)
varImpPlot (treeBag01)


###----------
Chosen.model <- treeBag01

    ###----------
    ### a. Pick a threshold and visualize fit
    threshold = .5    # pick a threshold

    #- Check the training set accuracy
    library(caret)     # install.packages('caret', repos='https://cran.case.edu/')
    Train.prob = predict(Chosen.model, type="prob")[,"Yes"]
    Train.pred = ifelse(Train.prob > threshold, "Yes", "No")  # Turn the fitted values to Up/Down using the threshold
    Test.prob = predict(Chosen.model, newdata=Test.set, type="prob")[,"Yes"]
    Test.pred = ifelse(Test.prob > threshold, "Yes", "No")  # Turn the fitted values to Up/Down using the threshold
    CM.train <- caret::confusionMatrix(factor(Train.pred), factor(as.matrix(Train.resp)), positive="Yes")
    CM.test  <- caret::confusionMatrix(factor(Test.pred),  factor(as.matrix(Test.resp)), positive="Yes")

    CM.train            # Training set result
    CM.train$table      # output just the table

    CM.train[["byClass"]][["Sensitivity"]]
    CM.train[["byClass"]][["Specificity"]]

    CM.test             # Testing set
    CM.test$table      # output just the table

    colSums(CM.test$table) / sum(colSums(CM.test$table))    # % of Actual Yes/No
    rowSums(CM.test$table) / sum(rowSums(CM.test$table))    # % of predicted Yes/No

    ###----------
    ### b. Output ROC curve and AUC for all threshold
    library(pROC)
    #- Training Set
    plot.roc(factor(as.matrix(Train.resp)),  Train.prob, levels=c("No", "Yes"))
    # point corresponding to CM.train
    abline(h=CM.train[["byClass"]][["Sensitivity"]], v=CM.train[["byClass"]][["Specificity"]], col="red")
    auc.train = auc(factor(as.matrix(Train.resp)), Train.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Train AUC=",round(auc.train, 3)))

    #- Test Set
    plot.roc(factor(as.matrix(Test.resp)),  Test.prob, levels=c("No", "Yes"))
    # point corresponding to CM.test
    abline(h=CM.test[["byClass"]][["Sensitivity"]], v=CM.test[["byClass"]][["Specificity"]], col="red")
    auc.test = auc(factor(as.matrix(Test.resp)), Test.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Test AUC=",round(auc.test, 3)))

    RandomForest = c(auc.train, auc.test)
    RandomForest

    # Training ROC and Test ROC side by side
    layout(matrix(1:2, 1, 2))
    plot.roc(factor(as.matrix(Train.resp)),  Train.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Train AUC=",round(auc.train, 3)))
    plot.roc(factor(as.matrix(Test.resp)),  Test.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Test AUC=",round(auc.test, 3)))
    layout(1)

    RandomForest


##----------------------------------------
## 8. Boosting (Train and Test only)

library(gbm)     # install.packages('gbm', repos='https://cran.case.edu/')
set.seed(my.seed)
treeBT01 = gbm::gbm(as.numeric(resp=="Yes") ~., data=Train.set,
                    distribution="bernoulli",   # for binary classification
                    n.trees=5000,
                    interaction.depth=4)

summary(treeBT01)

plot(treeBT01, i="Chol")
plot(treeBT01, i="MaxHR")


###----------
Chosen.model <- treeBT01

    ###----------
    ### a. Pick a threshold and visualize fit
    threshold = .5    # pick a threshold

    #- Check the training set accuracy
    library(caret)     # install.packages('caret', repos='https://cran.case.edu/')
    Train.prob = predict(Chosen.model, type="response", n.trees=1000)
    Train.pred = ifelse(Train.prob > threshold, "Yes", "No")  # Turn the fitted values to Up/Down using the threshold
    Test.prob = predict(Chosen.model, newdata=Test.set, type="response", n.trees=1000)
    Test.pred = ifelse(Test.prob > threshold, "Yes", "No")  # Turn the fitted values to Up/Down using the threshold
    CM.train <- caret::confusionMatrix(factor(Train.pred), factor(as.matrix(Train.resp)), positive="Yes")
    CM.test  <- caret::confusionMatrix(factor(Test.pred),  factor(as.matrix(Test.resp)), positive="Yes")

    CM.train            # Training set result
    CM.train$table      # output just the table

    CM.train[["byClass"]][["Sensitivity"]]
    CM.train[["byClass"]][["Specificity"]]

    CM.test             # Testing set
    CM.test$table      # output just the table

    colSums(CM.test$table) / sum(colSums(CM.test$table))    # % of Actual Yes/No
    rowSums(CM.test$table) / sum(rowSums(CM.test$table))    # % of predicted Yes/No

    ###----------
    ### b. Output ROC curve and AUC for all threshold
    library(pROC)
    #- Training Set
    plot.roc(factor(as.matrix(Train.resp)),  Train.prob, levels=c("No", "Yes"))
    # point corresponding to CM.train
    abline(h=CM.train[["byClass"]][["Sensitivity"]], v=CM.train[["byClass"]][["Specificity"]], col="red")
    auc.train = auc(factor(as.matrix(Train.resp)), Train.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Train AUC=",round(auc.train, 3)))

    #- Test Set
    plot.roc(factor(as.matrix(Test.resp)),  Test.prob, levels=c("No", "Yes"))
    # point corresponding to CM.test
    abline(h=CM.test[["byClass"]][["Sensitivity"]], v=CM.test[["byClass"]][["Specificity"]], col="red")
    auc.test = auc(factor(as.matrix(Test.resp)), Test.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Test AUC=",round(auc.test, 3)))

    Boosting = c(auc.train, auc.test)
    Boosting

    # Training ROC and Test ROC side by side
    layout(matrix(1:2, 1, 2))
    plot.roc(factor(as.matrix(Train.resp)),  Train.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Train AUC=",round(auc.train, 3)))
    plot.roc(factor(as.matrix(Test.resp)),  Test.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Test AUC=",round(auc.test, 3)))
    layout(1)

    Boosting


##----------------------------------------
## 9. Logistic Regression for comparison (Train vs Test only)

Fit03 <- glm(resp ~ .,  family=binomial, data=Train.set)
summary(Fit03)


###----------
Chosen.model <- Fit03

    ###----------
    #- a Output result for given threshold value
    threshold = .5 # pick a threshold

    #- Check the training set accuracy
    library(caret)
    Train.prob =predict(Chosen.model, type ="response")
    Train.pred = ifelse(Train.prob > threshold, "Yes", "No")  # Turn the fitted values to Up/Down using threshold of .5
    Test.prob = predict(Chosen.model, newdata=Test.set, type="response")
    Test.pred  = ifelse(Test.prob  > threshold, "Yes", "No")
    CM.train <- confusionMatrix(factor(Train.pred), factor(as.matrix(Train.resp)), positive="Yes")
    CM.test  <- confusionMatrix(factor(Test.pred),  factor(as.matrix(Test.resp)),  positive="Yes")

    CM.train            # Training set result
    CM.train$table      # output just the table

    CM.train[["byClass"]][["Sensitivity"]]
    CM.train[["byClass"]][["Specificity"]]

    CM.test             # Testing set
    CM.test$table      # output just the table

    ###----------
    #- b Output ROC curve and AUC for all threshold
    library(pROC)
    #- Training Set
    plot.roc(factor(as.matrix(Train.resp)),  Train.prob, levels=c("No", "Yes"))
    # point corresponding to CM.train
    abline(h=CM.train[["byClass"]][["Sensitivity"]], v=CM.train[["byClass"]][["Specificity"]], col="red")
    auc.train = auc(factor(as.matrix(Train.resp)), Train.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Train AUC=",round(auc.train, 3)))

    #- Test Set
    plot.roc(factor(as.matrix(Test.resp)),  Test.prob, levels=c("No", "Yes"))
    # point corresponding to CM.test
    abline(h=CM.test[["byClass"]][["Sensitivity"]], v=CM.test[["byClass"]][["Specificity"]], col="red")
    auc.test = auc(factor(as.matrix(Test.resp)), Test.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Test AUC=",round(auc.test, 3)))

    LogiReg = c(auc.train, auc.test)
    LogiReg

    layout(matrix(1:2, 1, 2))
    plot.roc(factor(as.matrix(Train.resp)),  Train.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Train AUC=",round(auc.train, 3)))
    plot.roc(factor(as.matrix(Test.resp)),  Test.prob, levels=c("No", "Yes"))
    text(.2, .2, paste("Test AUC=",round(auc.test, 3)))
    layout(1)

    LogiReg


##—————————————- ## Compare Performances

With my.seed=7211

Compare (auc.train, auc.test)

LogiReg 0.9437945 0.8964497

GrowAndPrune 0.8898077 0.8853550

Bagging 0.8768657 0.9215976

RandomForest 0.8659742 0.9260355

Boosting 1.0000000 0.8727811

Compare (auc.train, auc.test)

LogiReg GrowAndPrune Bagging RandomForest Boosting

##—————————————- ## 10. Threshold Picker

#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(0.01,.99,.01) # grid for threshold cost=0 library(caret) # for confusionMatrix

Have to be switched depending on the model

Chosen.model <- pruned1; Test.prob = predict(Chosen.model, Test.set, type=“vector”)[,“Yes”] # Chosen.model <- treeBag01; Test.prob = predict(Chosen.model, newdata=Test.set, type=“prob”)[,“Yes”] # Chosen.model <- treeBT01; Test.prob = predict(Chosen.model, newdata=Test.set, type=“response”, n.trees=1000)

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)] ```