Default - 4. Logistic Regression Analysis
Nao Mimoto – Class Web Page
1. Load data from ISLR package
As we did in Ch2-Lab.
library(ISLR) # install.packages('ISLR', repos='https://cran.case.edu/')
head(Default) # see the data## default student balance income
## 1 No No 729.5265 44361.625
## 2 No Yes 817.1804 12106.135
## 3 No No 1073.5492 31767.139
## 4 No No 529.2506 35704.494
## 5 No No 785.6559 38463.496
## 6 No Yes 919.5885 7491.559## [1] 10000 4# ?Default # see explanation for variables
# A data frame with 10000 observations on the following 4 variables.
#
# 1 default (factor) A factor with levels No and Yes indicating whether the customer defaulted on their debt
# 2 student (factor) A factor with levels No and Yes indicating whether the customer is a student
# 3 balance (numeric) The average balance that the customer has remaining on their credit card after making their monthly payment
# 4 income (numeric) Income of customer
#
# column "default" is the respoinse variable.
#--- Turn the data into tibble ---
library(tidyverse) # install.packages("tidyverse")
Default <- as_tibble(Default)
Default## # A tibble: 10,000 × 4
## default student balance income
## <fct> <fct> <dbl> <dbl>
## 1 No No 730. 44362.
## 2 No Yes 817. 12106.
## 3 No No 1074. 31767.
## 4 No No 529. 35704.
## 5 No No 786. 38463.
## 6 No Yes 920. 7492.
## 7 No No 826. 24905.
## 8 No Yes 809. 17600.
## 9 No No 1161. 37469.
## 10 No No 0 29275.
## # ℹ 9,990 more rowsDefault2 <- Default %>%
rename(resp=default) %>% # Rename "default" column as "resp"
relocate(resp) # move "resp" columnm to 1st
Default2## # A tibble: 10,000 × 4
## resp student balance income
## <fct> <fct> <dbl> <dbl>
## 1 No No 730. 44362.
## 2 No Yes 817. 12106.
## 3 No No 1074. 31767.
## 4 No No 529. 35704.
## 5 No No 786. 38463.
## 6 No Yes 920. 7492.
## 7 No No 826. 24905.
## 8 No Yes 809. 17600.
## 9 No No 1161. 37469.
## 10 No No 0 29275.
## # ℹ 9,990 more rows## resp student balance income
## No :9667 No :7056 Min. : 0.0 Min. : 772
## Yes: 333 Yes:2944 1st Qu.: 481.7 1st Qu.:21340
## Median : 823.6 Median :34553
## Mean : 835.4 Mean :33517
## 3rd Qu.:1166.3 3rd Qu.:43808
## Max. :2654.3 Max. :73554## [1] 0## [1] 10000 4## [1] 10000 42. Data Separation (Copied from Ch3-Lab)
# 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: 10000
## training data: 1666 x 5 = 8330
## size of test data: 1670
##
## 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. Logistic Regression on Training Set
3.1 Model1 Just balance
##
## Call:
## glm(formula = resp ~ balance, family = binomial, data = Train.set)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.071e+01 4.013e-01 -26.69 <2e-16 ***
## balance 5.528e-03 2.445e-04 22.61 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2423.5 on 8329 degrees of freedom
## Residual deviance: 1328.3 on 8328 degrees of freedom
## AIC: 1332.3
##
## Number of Fisher Scoring iterations: 8 # Plot the result for visualization
plot(Train.set$balance, Train.set$resp=="Yes", ylab="Resp", xlab="Balance")
lines(Train.set$balance, Fit01$fitted, lwd=2, col="red", type="p")
3.2 Model2 All three
##
## Call:
## glm(formula = resp ~ ., family = binomial, data = Train.set)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.062e+01 5.320e-01 -19.971 < 2e-16 ***
## studentYes -7.591e-01 2.554e-01 -2.973 0.00295 **
## balance 5.759e-03 2.572e-04 22.393 < 2e-16 ***
## income -4.756e-06 8.992e-06 -0.529 0.59684
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2423.5 on 8329 degrees of freedom
## Residual deviance: 1310.6 on 8326 degrees of freedom
## AIC: 1318.6
##
## Number of Fisher Scoring iterations: 8## (Intercept) studentYes balance income
## -1.062473e+01 -7.591343e-01 5.759184e-03 -4.756460e-06## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.062473e+01 5.320173e-01 -19.9706430 9.916855e-89
## studentYes -7.591343e-01 2.553536e-01 -2.9728749 2.950247e-03
## balance 5.759184e-03 2.571868e-04 22.3929968 4.605582e-111
## income -4.756460e-06 8.992349e-06 -0.5289453 5.968434e-01## (Intercept) studentYes balance income
## 9.916855e-89 2.950247e-03 4.605582e-111 5.968434e-01 #--- Plot the result for visualization
plot(Train.set$balance, Train.set$resp=="Yes", ylab="Resp", xlab="Balance")
lines(Train.set$balance, Fit02$fitted, lwd=2, col="red", type="p")
3.3 Model3 Remove income
##
## Call:
## glm(formula = resp ~ balance + student, family = binomial, data = Train.set)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.081e+01 4.099e-01 -26.366 < 2e-16 ***
## balance 5.754e-03 2.569e-04 22.396 < 2e-16 ***
## studentYes -6.535e-01 1.600e-01 -4.084 4.42e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2423.5 on 8329 degrees of freedom
## Residual deviance: 1310.9 on 8327 degrees of freedom
## AIC: 1316.9
##
## Number of Fisher Scoring iterations: 8 #--- Plot the result for visualization
plot(Train.set$balance, Train.set$resp=="Yes", ylab="Resp", xlab="Balance")
lines(Train.set$balance, Fit03$fitted, lwd=2, col="red", type="p")
# Two curves are for Student= Yes/No
# Turn Student=Yes into Blue
ix = (Train.set$student=="Yes")
lines(Train.set$balance[ix], Fit03$fitted[ix], lwd=2, col="blue", type="p")
###----------
Fit00 = Fit03
# Extract fitted response (training)
Train.prob =predict(Fit00, type ="response")
# Predict in Test Set
Test.prob = predict(Fit00, newdata=Test.set, type="response")
###----------
### 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.9493302 0.9506749 ###----------
### 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 8023 192
## Yes 31 84
##
## Accuracy : 0.9732
## 95% CI : (0.9695, 0.9766)
## No Information Rate : 0.9669
## P-Value [Acc > NIR] : 0.0004646
##
## Kappa : 0.4183
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.30435
## Specificity : 0.99615
## Pos Pred Value : 0.73043
## Neg Pred Value : 0.97663
## Prevalence : 0.03313
## Detection Rate : 0.01008
## Detection Prevalence : 0.01381
## Balanced Accuracy : 0.65025
##
## 'Positive' Class : Yes
## ## Confusion Matrix and Statistics
##
## Reference
## Prediction No Yes
## No 1605 36
## Yes 8 21
##
## Accuracy : 0.9737
## 95% CI : (0.9648, 0.9808)
## No Information Rate : 0.9659
## P-Value [Acc > NIR] : 0.04196
##
## Kappa : 0.4763
##
## Mcnemar's Test P-Value : 4.693e-05
##
## Sensitivity : 0.36842
## Specificity : 0.99504
## Pos Pred Value : 0.72414
## Neg Pred Value : 0.97806
## Prevalence : 0.03413
## Detection Rate : 0.01257
## Detection Prevalence : 0.01737
## Balanced Accuracy : 0.68173
##
## 'Positive' Class : Yes
## ## Reference
## Prediction No Yes
## No 8023 192
## Yes 31 84## Reference
## Prediction No Yes
## No 1605 36
## Yes 8 21## [1] 0.3684211## [1] 0.9950403 ## 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.96586826 0.03413174## No Yes
## 0.98263473 0.01736527 ###----------
### 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.9493302 0.95067494. Search for the Best Model using CV
AUCs <- MSE.valid <- 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)
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.9496 0.9477
## [2,] 0.9480 0.9555
## [3,] 0.9496 0.9468
## [4,] 0.9470 0.9572
## [5,] 0.9525 0.9377## Av.Train AUC Av.Valid AUC
## 0.94934 0.94898###
### Make decision about best model based on Av Valid AUC
### (no tuning of hyperparameter with logistic regression)
###5. Final Training/Test fit using best model
### Best model
Fit03 <- glm(resp ~ balance + student, family=binomial, data=Train.set)
summary(Fit03)##
## Call:
## glm(formula = resp ~ balance + student, family = binomial, data = Train.set)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.081e+01 4.099e-01 -26.366 < 2e-16 ***
## balance 5.754e-03 2.569e-04 22.396 < 2e-16 ***
## studentYes -6.535e-01 1.600e-01 -4.084 4.42e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2423.5 on 8329 degrees of freedom
## Residual deviance: 1310.9 on 8327 degrees of freedom
## AIC: 1316.9
##
## Number of Fisher Scoring iterations: 8###----------
Fit00 = Fit03
# Extract fitted response (training)
Train.prob =predict(Fit00, type ="response")
# Predict in Test Set
Test.prob = predict(Fit00, newdata=Test.set, type="response")
###----------
### 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.9493302 0.9506749 ###----------
### 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 8023 192
## Yes 31 84
##
## Accuracy : 0.9732
## 95% CI : (0.9695, 0.9766)
## No Information Rate : 0.9669
## P-Value [Acc > NIR] : 0.0004646
##
## Kappa : 0.4183
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.30435
## Specificity : 0.99615
## Pos Pred Value : 0.73043
## Neg Pred Value : 0.97663
## Prevalence : 0.03313
## Detection Rate : 0.01008
## Detection Prevalence : 0.01381
## Balanced Accuracy : 0.65025
##
## 'Positive' Class : Yes
## ## Confusion Matrix and Statistics
##
## Reference
## Prediction No Yes
## No 1605 36
## Yes 8 21
##
## Accuracy : 0.9737
## 95% CI : (0.9648, 0.9808)
## No Information Rate : 0.9659
## P-Value [Acc > NIR] : 0.04196
##
## Kappa : 0.4763
##
## Mcnemar's Test P-Value : 4.693e-05
##
## Sensitivity : 0.36842
## Specificity : 0.99504
## Pos Pred Value : 0.72414
## Neg Pred Value : 0.97806
## Prevalence : 0.03413
## Detection Rate : 0.01257
## Detection Prevalence : 0.01737
## Balanced Accuracy : 0.68173
##
## 'Positive' Class : Yes
## ## Reference
## Prediction No Yes
## No 8023 192
## Yes 31 84## Reference
## Prediction No Yes
## No 1605 36
## Yes 8 21## [1] 0.3684211## [1] 0.9950403 ## 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.96586826 0.03413174## No Yes
## 0.98263473 0.01736527 ###----------
### 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.9493302 0.9506749 ###----------
### 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.9496 0.9477
## [2,] 0.9480 0.9555
## [3,] 0.9496 0.9468
## [4,] 0.9470 0.9572
## [5,] 0.9525 0.9377## Av.Train AUC Av.Valid AUC
## 0.94934 0.94898## [1] 0.9493302 0.95067496. Determine the Best Threshold
# Decide on Threshould penalty
#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.
## 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.
## [1] 0.00 0.00 0.75 0.25## [1] 72## [1] 12.75## [1] 0.72CODE ONLY
##----------------------------------------
## 1. Load data from ISLR package
library(ISLR) # install.packages('ISLR', repos='https://cran.case.edu/')
head(Default) # see the data
dim(Default)
# ?Default # see explanation for variables
# A data frame with 10000 observations on the following 4 variables.
#
# 1 default (factor) A factor with levels No and Yes indicating whether the customer defaulted on their debt
# 2 student (factor) A factor with levels No and Yes indicating whether the customer is a student
# 3 balance (numeric) The average balance that the customer has remaining on their credit card after making their monthly payment
# 4 income (numeric) Income of customer
#
# column "default" is the respoinse variable.
#--- Turn the data into tibble ---
library(tidyverse) # install.packages("tidyverse")
Default <- as_tibble(Default)
Default
Default2 <- Default %>%
rename(resp=default) %>% # Rename "default" column as "resp"
relocate(resp) # move "resp" columnm to 1st
Default2
Orig <- Default2
#- 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 (Copied from Ch3-Lab)
# 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]]]]
##----------------------------------------
## 3. Logistic Regression on Training Set
###----------
### 3.1 Model1 Just balance
Fit01 <- glm(resp ~ balance, family=binomial, data=Train.set)
summary(Fit01)
# Plot the result for visualization
plot(Train.set$balance, Train.set$resp=="Yes", ylab="Resp", xlab="Balance")
lines(Train.set$balance, Fit01$fitted, lwd=2, col="red", type="p")
###----------
### 3.2 Model2 All three
Fit02 <- glm(resp ~ ., family=binomial, data=Train.set)
summary(Fit02)
coef(Fit02)
summary(Fit02)$coef # extract estimated parameters
summary(Fit02)$coef[,4] # extract 4th column (P-value) of above
#--- Plot the result for visualization
plot(Train.set$balance, Train.set$resp=="Yes", ylab="Resp", xlab="Balance")
lines(Train.set$balance, Fit02$fitted, lwd=2, col="red", type="p")
###----------
### 3.3 Model3 Remove income
Fit03 <- glm(resp ~ balance + student, family=binomial, data=Train.set)
summary(Fit03)
#--- Plot the result for visualization
plot(Train.set$balance, Train.set$resp=="Yes", ylab="Resp", xlab="Balance")
lines(Train.set$balance, Fit03$fitted, lwd=2, col="red", type="p")
# Two curves are for Student= Yes/No
# Turn Student=Yes into Blue
ix = (Train.set$student=="Yes")
lines(Train.set$balance[ix], Fit03$fitted[ix], lwd=2, col="blue", type="p")
#- Extract fitted response (training)
Train.prob =predict(Fit03, type ="response")
head(Train.prob)
#- Predict in Test Set
Test.prob = predict(Fit03, 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)
##----------------------------------------
## 4. Logistic Regression with CV
AUCs <- MSE.valid <- 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)
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
###
### Make decision about best model based on Av Valid AUC
### (no tuning of hyperparameter here)
###
##----------------------------------------
## 5. Final Training/Test fit using best model
### Best model
Fit03 <- glm(resp ~ balance + student, family=binomial, data=Train.set)
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)
##----------------------------------------
## 6. Determine the Best Threshold
# Decide on Threshould penalty
#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)]