heart - 4. Logistic Regression
Nao Mimoto – Class Web Page
1. Load data
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))
# Low p-value from Chi-sq test of association (see heart 0)
# 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"
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## resp Age Sex ChestPain RestBP
## No :164 Min. :29.00 0: 97 asymptomatic:144 Min. : 94.0
## Yes:139 1st Qu.:48.00 1:206 nonanginal : 86 1st Qu.:120.0
## Median :56.00 nontypical : 50 Median :130.0
## Mean :54.44 typical : 23 Mean :131.7
## 3rd Qu.:61.00 3rd Qu.:140.0
## Max. :77.00 Max. :200.0
##
## Chol Fbs RestECG MaxHR ExAng Oldpeak
## Min. :126.0 0:258 0:151 Min. : 71.0 0:204 Min. :0.00
## 1st Qu.:211.0 1: 45 1: 4 1st Qu.:133.5 1: 99 1st Qu.:0.00
## Median :241.0 2:148 Median :153.0 Median :0.80
## Mean :246.7 Mean :149.6 Mean :1.04
## 3rd Qu.:275.0 3rd Qu.:166.0 3rd Qu.:1.60
## Max. :564.0 Max. :202.0 Max. :6.20
##
## Slope Ca Thal
## Min. :1.000 Min. :0.0000 fixed : 18
## 1st Qu.:1.000 1st Qu.:0.0000 normal :166
## Median :2.000 Median :0.0000 reversable:117
## Mean :1.601 Mean :0.6722 NA's : 2
## 3rd Qu.:2.000 3rd Qu.:1.0000
## Max. :3.000 Max. :3.0000
## NA's :4## [1] 303 14## [1] 6## [1] 303 14## [1] 297 142. 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
CreateCV(Orig, numFolds=5, seed=5480)##
## number of folds: 5
## number of obs: 297
## training data: 49 x 5 = 245
## size of test data: 52
##
## 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]]]]
# Visualization doesn't help much any more
ix_no = (Train.set$resp=="No")
ix_yes = (Train.set$resp=="Yes")
plot(Train.set$RestBP[ix_no], Train.set$Chol[ix_no], xlab="RestBp", ylab="Chol")
lines(Train.set$RestBP[ix_yes], Train.set$Chol[ix_yes], col="red", type="p")
3. Logistic Regression with CV
AUCs <- matrix(0, 5, 2)
colnames(AUCs) = c("Train AUC", "Valid AUC")
for (k in 1:5) {
Fit00 <- glm(resp ~. , family=binomial, data=CV.train[[k]]) # <----- Change model here
#- Extract fitted response (training)
Train.prob =predict(Fit00, type ="response") # fitted responses
#- Predict in Validation Set
Valid.prob = predict(Fit00, newdata=CV.valid[[k]], type="response")
#- Check the training set accuracy
library(caret)
library(pROC)
# Train.pred = ifelse(Train.prob > threshold, "Yes", "No") # Turn the fitted values to Up/Down using threshold of .5
# Valid.pred = ifelse(Valid.prob > threshold, "Yes", "No")
# CM.train <- confusionMatrix(factor(Train.pred), factor(as.matrix(CV.train.resp[[k]])), positive="Yes")
# CM.valid <- confusionMatrix(factor(Valid.pred), factor(as.matrix(CV.valid.resp[[k]])), positive="Yes")
AUCs[k,] <- round(c(auc(factor(as.matrix(CV.train.resp[[k]])), Train.prob, levels=c("No", "Yes")),
auc(factor(as.matrix(CV.valid.resp[[k]])), Valid.prob, levels=c("No", "Yes"))), 4)
}
AUCs## Train AUC Valid AUC
## [1,] 0.9235 0.9044
## [2,] 0.9288 0.9211
## [3,] 0.9159 0.9398
## [4,] 0.9451 0.7525
## [5,] 0.9329 0.8655## Av.Train AUC Av.Valid AUC
## 0.92924 0.87666### Run this section with different models for Fit00.
### Make decision about the best model based on Av Valid AUC
###4. Final Training/Test fit using best model
### Best model (of your choice)
Fit03 <- glm(resp ~ ., family=binomial, data=Train.set) # <===== needs to be changed to your model
summary(Fit03)##
## Call:
## glm(formula = resp ~ ., family = binomial, data = Train.set)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.459914 3.282489 -1.359 0.17424
## Age -0.008370 0.028472 -0.294 0.76877
## Sex1 1.176484 0.550751 2.136 0.03267 *
## ChestPainnonanginal -1.457325 0.527059 -2.765 0.00569 **
## ChestPainnontypical -0.731788 0.593656 -1.233 0.21770
## ChestPaintypical -1.640239 0.749599 -2.188 0.02866 *
## RestBP 0.019652 0.012005 1.637 0.10163
## Chol 0.005199 0.004091 1.271 0.20386
## Fbs1 -0.553422 0.630513 -0.878 0.38009
## RestECG1 0.788465 2.368506 0.333 0.73921
## RestECG2 0.480117 0.414701 1.158 0.24697
## MaxHR -0.016947 0.013010 -1.303 0.19271
## ExAng1 0.896079 0.489790 1.830 0.06732 .
## Oldpeak 0.444619 0.243525 1.826 0.06789 .
## Slope 0.652920 0.390366 1.673 0.09441 .
## Ca 1.484943 0.323065 4.596 4.3e-06 ***
## Thalnormal 0.011247 0.837676 0.013 0.98929
## Thalreversable 1.152122 0.811619 1.420 0.15574
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 334.62 on 244 degrees of freedom
## Residual deviance: 167.94 on 227 degrees of freedom
## AIC: 203.94
##
## Number of Fisher Scoring iterations: 6###----------
Fit00 = Fit03
# Extract fitted response (training)
Train.prob =predict(Fit00, type ="response")
head(Train.prob)## 1 2 3 4 5 6
## 0.26420804 0.99885545 0.99584106 0.30174141 0.05237129 0.96995924## 1 2 3 4 5 6
## 0.04445500 0.24093499 0.99582351 0.92737561 0.22136791 0.03733577 ###----------
### a. Output ROC curve and AUC without threshold
library(pROC)
layout(matrix(1:2, 1, 2))
plot.roc(factor(as.matrix(Train.resp)), Train.prob, levels=c("No", "Yes"))
auc.train = auc(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"))
auc.test = auc(factor(as.matrix(Test.resp)), Test.prob, levels=c("No", "Yes"))
text(.2, .2, paste("Train AUC=",round(auc.test, 3)))
## [1] 0.9268027 0.9468750 ###----------
### b. Output Confusion Matrix for given threshold
threshold = .5 # pick a threshold
#- Check the training set accuracy
library(caret)
Train.pred = ifelse(Train.prob > threshold, "Yes", "No") # Turn the fitted values to Up/Down using threshold of .5
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")
##===== Training CM =====
CM.train # Training set result## Confusion Matrix and Statistics
##
## Reference
## Prediction No Yes
## No 130 23
## Yes 10 82
##
## Accuracy : 0.8653
## 95% CI : (0.8161, 0.9054)
## No Information Rate : 0.5714
## P-Value [Acc > NIR] : < 2e-16
##
## Kappa : 0.7207
##
## Mcnemar's Test P-Value : 0.03671
##
## Sensitivity : 0.7810
## Specificity : 0.9286
## Pos Pred Value : 0.8913
## Neg Pred Value : 0.8497
## Prevalence : 0.4286
## Detection Rate : 0.3347
## Detection Prevalence : 0.3755
## Balanced Accuracy : 0.8548
##
## 'Positive' Class : Yes
## ## Confusion Matrix and Statistics
##
## Reference
## Prediction No Yes
## No 16 3
## Yes 4 29
##
## Accuracy : 0.8654
## 95% CI : (0.7421, 0.9441)
## No Information Rate : 0.6154
## P-Value [Acc > NIR] : 7.047e-05
##
## Kappa : 0.7129
##
## Mcnemar's Test P-Value : 1
##
## Sensitivity : 0.9062
## Specificity : 0.8000
## Pos Pred Value : 0.8788
## Neg Pred Value : 0.8421
## Prevalence : 0.6154
## Detection Rate : 0.5577
## Detection Prevalence : 0.6346
## Balanced Accuracy : 0.8531
##
## 'Positive' Class : Yes
## ## Reference
## Prediction No Yes
## No 130 23
## Yes 10 82## Reference
## Prediction No Yes
## No 16 3
## Yes 4 29## [1] 0.90625## [1] 0.8 ## Note: Confusion Matrix
## Reference
## Prediction No Yes
## No 1437 62 [Specificity][ ] = TrueNeg / sum of col
## Yes 0 1 [ ][Sensitivity] = TruePos / sum of col
## High threshold -> Specificity goes to 1. Very conservative. Everyone is flagged neg.
## Low threshold -> Sensitivity goes to 1. Very agressive. Eveyone is flagged pos.
colSums(CM.test$table) / sum(colSums(CM.test$table)) # % of Actual Yes/No## No Yes
## 0.3846154 0.6153846## No Yes
## 0.3653846 0.6346154 ###----------
### c. ROC and corresponding point for the threshold chosen
library(pROC)
layout(matrix(1:2, 1, 2))
threshold # current value of threshold## [1] 0.5 #- 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")
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")
text(.2, .2, paste("Test AUC=",round(auc.test, 3)))
## [1] 0.9268027 0.9468750 ###----------
### d. Comparison of Validation AUC and Test AUC
plot(AUCs[,2], col="red", ylim=c(.5,1),
main="Av. CV training/valid (blk/red) AUC and Test AUC (line)")
lines(AUCs[,1], type="p")
abline(h=auc.test)
## Train AUC Valid AUC
## [1,] 0.9235 0.9044
## [2,] 0.9288 0.9211
## [3,] 0.9159 0.9398
## [4,] 0.9451 0.7525
## [5,] 0.9329 0.8655## Av.Train AUC Av.Valid AUC
## 0.92924 0.87666## [1] 0.9268027 0.94687505. Determine the Best Threshold
## Decide on penalty for each cell of confusion matrix
#cost.list = c(1,1,1,1)/4 # order of (TP, TN, FP, FN)
cost.list = c(0,0,3,1)/4 # 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
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.
## [1] 0.00 0.00 0.75 0.25## [1] 69## [1] 2.5## [1] 0.69CODE ONLY
##----------------------------------------
## 1. Load data
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))
# Low p-value from Chi-sq test of association (see heart 0)
# 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"
print(Heart2, width=1000)
table(Heart2$resp) # note that first level is "No"
Orig <- Heart2
#- 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
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
CreateCV(Orig, numFolds=5, seed=5480)
# 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]]]]
# Visualization doesn't help much any more
ix_no = (Train.set$resp=="No")
ix_yes = (Train.set$resp=="Yes")
plot(Train.set$RestBP[ix_no], Train.set$Chol[ix_no], xlab="RestBp", ylab="Chol")
lines(Train.set$RestBP[ix_yes], Train.set$Chol[ix_yes], col="red", type="p")
##----------------------------------------
## 3. Logistic Regression with CV
AUCs <- matrix(0, 5, 2)
colnames(AUCs) = c("Train AUC", "Valid AUC")
for (k in 1:5) {
Fit00 <- glm(resp ~. , family=binomial, data=CV.train[[k]]) # <----- Change model here
#- Extract fitted response (training)
Train.prob =predict(Fit00, type ="response") # fitted responses
#- Predict in Validation Set
Valid.prob = predict(Fit00, newdata=CV.valid[[k]], type="response")
#- Check the training set accuracy
library(caret)
library(pROC)
# Train.pred = ifelse(Train.prob > threshold, "Yes", "No") # Turn the fitted values to Up/Down using threshold of .5
# Valid.pred = ifelse(Valid.prob > threshold, "Yes", "No")
# CM.train <- confusionMatrix(factor(Train.pred), factor(as.matrix(CV.train.resp[[k]])), positive="Yes")
# CM.valid <- confusionMatrix(factor(Valid.pred), factor(as.matrix(CV.valid.resp[[k]])), positive="Yes")
AUCs[k,] <- round(c(auc(factor(as.matrix(CV.train.resp[[k]])), Train.prob, levels=c("No", "Yes")),
auc(factor(as.matrix(CV.valid.resp[[k]])), 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
### Run this section with different models for Fit00.
### Make decision about the best model based on Av Valid AUC
###
##----------------------------------------
## 4. Final Training/Test fit using best model
### Best model (of your choice)
Fit03 <- glm(resp ~ ., family=binomial, data=Train.set) # <===== needs to be changed to your model
summary(Fit03)
###----------
Fit00 = Fit03
# Extract fitted response (training)
Train.prob =predict(Fit00, type ="response")
head(Train.prob)
# Predict in Test Set
Test.prob = predict(Fit00, newdata=Test.set, type="response")
head(Test.prob)
###----------
### a. Output ROC curve and AUC without threshold
library(pROC)
layout(matrix(1:2, 1, 2))
plot.roc(factor(as.matrix(Train.resp)), Train.prob, levels=c("No", "Yes"))
auc.train = auc(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"))
auc.test = auc(factor(as.matrix(Test.resp)), Test.prob, levels=c("No", "Yes"))
text(.2, .2, paste("Train AUC=",round(auc.test, 3)))
layout(1)
c(auc.train, auc.test)
###----------
### b. Output Confusion Matrix for given threshold
threshold = .5 # pick a threshold
#- Check the training set accuracy
library(caret)
Train.pred = ifelse(Train.prob > threshold, "Yes", "No") # Turn the fitted values to Up/Down using threshold of .5
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")
##===== Training CM =====
CM.train # Training set result
##===== Test CM =====
CM.test # Testing set
##----- Just CMs -----
CM.train$table # train CM
CM.test$table # test CM
## If you want to extract specific metric
CM.test[["byClass"]][["Sensitivity"]]
CM.test[["byClass"]][["Specificity"]]
## Note: Confusion Matrix
## Reference
## Prediction No Yes
## No 1437 62 [Specificity][ ] = TrueNeg / sum of col
## Yes 0 1 [ ][Sensitivity] = TruePos / sum of col
## High threshold -> Specificity goes to 1. Very conservative. Everyone is flagged neg.
## Low threshold -> Sensitivity goes to 1. Very agressive. Eveyone is flagged pos.
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. ROC and corresponding point for the threshold chosen
library(pROC)
layout(matrix(1:2, 1, 2))
threshold # current value of threshold
#- 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")
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")
text(.2, .2, paste("Test AUC=",round(auc.test, 3)))
layout(1)
c(auc.train, auc.test)
###----------
### d. Comparison of Validation AUC and Test AUC
plot(AUCs[,2], col="red", ylim=c(.5,1),
main="Av. CV training/valid (blk/red) AUC and Test AUC (line)")
lines(AUCs[,1], type="p")
abline(h=auc.test)
AUCs
Av.AUCs
c(auc.train, auc.test)
##----------------------------------------
## 5. Determine the Best Threshold
## Decide on penalty for each cell of confusion matrix
#cost.list = c(1,1,1,1)/4 # order of (TP, TN, FP, FN)
cost.list = c(0,0,3,1)/4 # 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
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)]