#=================================================================================
# Using Accuracy Criterion for Credit Classification Can Cause Disaster for Banks. 
#=================================================================================
 

# Load R packages: 
rm(list = ls())
library(cutpointr)
library(tidyverse)
library(caret)

# Load GermanCredit Data: 
data("GermanCredit")
df <- GermanCredit 

# Split data: 
set.seed(1)
id <- createDataPartition(y = df$Class, p = 0.7, list = FALSE)
train <- df[id, ]
test <- df[-id, ]

# Set conditions for training Logistic Model: 

number <- 5
repeats <- 5
n <- number*repeats

set.seed(1)
train.control <- trainControl(method = "repeatedcv", 
                              number = number,
                              repeats = repeats, 
                              classProbs = TRUE,
                              allowParallel = TRUE, 
                              summaryFunction = twoClassSummary)


# Trani Logistic Model with train data: 
my_logit <- train(Class ~., 
                  data = train, 
                  method = "glm", 
                  trControl = train.control)


# Use model for predicting PD: 
pd <- predict(my_logit, test, type = "prob") %>% pull(Bad)


# Calculate optimal cutoff by sensitivity and accuracy criterion: 
m <- cutpointr(x = pd, class = test$Class, metric = sens_constrain, pos_class = "Bad")
n <- cutpointr(x = pd, class = test$Class, metric = accuracy, pos_class = "Bad")

# Optimal cutoff: 
t1 <- m$optimal_cutpoint # maximize sensitivity. 
t2 <- n$optimal_cutpoint # maximize accuracy. 

# Plot ROC curve with optimal cutoff: 
gridExtra::grid.arrange(plot(m), plot(n))

# X is unchanged regardless of method selected: 

m$AUC

## [1] 0.7504762

n$AUC

## [1] 0.7504762

# Function for labelling credit applications: 

label_predicted <- function(cutoff) {
  y <- case_when(pd >= cutoff ~ "Bad", TRUE ~ "Good") %>% as.factor()
  return(y)
}

# Confution maxtrix: 
confusionMatrix(label_predicted(t1), test$Class, positive = "Bad") # If t1 is selected for classification. 

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Bad Good
##       Bad   76  102
##       Good  14  108
##                                           
##                Accuracy : 0.6133          
##                  95% CI : (0.5557, 0.6687)
##     No Information Rate : 0.7             
##     P-Value [Acc > NIR] : 0.9995          
##                                           
##                   Kappa : 0.2804          
##                                           
##  Mcnemar's Test P-Value : 6.597e-16       
##                                           
##             Sensitivity : 0.8444          
##             Specificity : 0.5143          
##          Pos Pred Value : 0.4270          
##          Neg Pred Value : 0.8852          
##              Prevalence : 0.3000          
##          Detection Rate : 0.2533          
##    Detection Prevalence : 0.5933          
##       Balanced Accuracy : 0.6794          
##                                           
##        'Positive' Class : Bad             
## 

confusionMatrix(label_predicted(t2), test$Class, positive = "Bad") # If t2 is selected for classification. 

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Bad Good
##       Bad   38   23
##       Good  52  187
##                                         
##                Accuracy : 0.75          
##                  95% CI : (0.697, 0.798)
##     No Information Rate : 0.7           
##     P-Value [Acc > NIR] : 0.032245      
##                                         
##                   Kappa : 0.3444        
##                                         
##  Mcnemar's Test P-Value : 0.001224      
##                                         
##             Sensitivity : 0.4222        
##             Specificity : 0.8905        
##          Pos Pred Value : 0.6230        
##          Neg Pred Value : 0.7824        
##              Prevalence : 0.3000        
##          Detection Rate : 0.1267        
##    Detection Prevalence : 0.2033        
##       Balanced Accuracy : 0.6563        
##                                         
##        'Positive' Class : Bad           
## 

# Results of Classification by the two methods: 

test %>% 
  mutate(Class_max_sen = label_predicted(t1),
         Class_max_accuracy = label_predicted(t2)) -> test_t1_t2

test_t1_t2 %>% 
  filter(Class_max_sen == "Good", Class == "Good") %>% 
  select(Class, Class_max_sen, Amount) -> df_Good_max_sen

test_t1_t2 %>% 
  filter(Class_max_sen == "Good", Class == "Bad") %>% 
  select(Class, Class_max_sen, Amount) -> df_Bad_max_sen

test_t1_t2 %>% 
  filter(Class_max_accuracy == "Good", Class == "Good") %>% 
  select(Class, Class_max_accuracy, Amount) -> df_Good_max_acc

test_t1_t2 %>% 
  filter(Class_max_accuracy == "Good", Class == "Bad") %>% 
  select(Class, Class_max_accuracy, Amount) -> df_Bad_max_acc


# Calculate profit if interest rate = 30%: 
ir <- 0.3
sum(ir*df_Good_max_sen$Amount) - sum(df_Bad_max_sen$Amount)

## [1] 42215.5

sum(ir*df_Good_max_acc$Amount) - sum(df_Bad_max_acc$Amount)

## [1] -21770.7

# Function for calculating profit: 

profit <- function(ir) {
  pro1 <- sum(ir*df_Good_max_sen$Amount) - sum(df_Bad_max_sen$Amount)
  pro2 <- sum(ir*df_Good_max_acc$Amount) - sum(df_Bad_max_acc$Amount)
  return(data.frame(IR = ir, Pro = c(pro1, pro2), Method = c("MaxSen", "MaxAcc")))
}


# Compare profit: 

lapply(seq(0.15, 0.5, by = 0.01), profit) -> profit_list 

do.call("rbind", profit_list) %>% 
  mutate(Pro = Pro / 1000) %>% 
  ggplot(aes(IR, Pro, color = Method)) + 
  geom_line() + 
  geom_point() + 
  scale_x_continuous(breaks = seq(0.15, 0.5, by = 0.05), labels = scales::percent) + 
  labs(x = "Interest Rate", y = "Profit", 
       title = "Profit Comparision by Classification Method") + 
  hrbrthemes::theme_modern_rc()

LS0tDQp0aXRsZTogIlVzaW5nIEFjY3VyYWN5IENyaXRlcmlvbiBmb3IgQ3JlZGl0IENsYXNzaWZpY2F0aW9uIENhbiBDYXVzZSBEaXNhc3RlciBmb3IgQmFua3MiDQphdXRob3I6ICJOZ3V5ZW4gQ2hpIER1bmciDQpzdWJ0aXRsZTogIlIgZm9yIFBsZWFzdXJlIg0Kb3V0cHV0Og0KICBodG1sX2RvY3VtZW50Og0KICAgIGNvZGVfZG93bmxvYWQ6IHllcw0KICAgIGNvZGVfZm9sZGluZzogaGlkZQ0KICAgIGhpZ2hsaWdodDogemVuYnVybg0KICAgIHRoZW1lOiBmbGF0bHkNCiAgICB0b2M6IHllcw0KICAgIHRvY19mbG9hdDogeWVzDQogIHdvcmRfZG9jdW1lbnQ6DQogICAgdG9jOiB5ZXMNCi0tLQ0KDQpgYGB7ciBzZXR1cCxpbmNsdWRlPUZBTFNFfQ0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFLCB3YXJuaW5nID0gRkFMU0UsIG1lc3NhZ2UgPSBGQUxTRSwgZmlnLnJldGluYT0yKQ0KYGBgDQoNCiMgV2h5IEFjY3VyYWN5IElzIGEgQmFkIE1lYXN1cmUgZm9yIENyZWRpdCBDbGFzc2lmaWNhdGlvbj8gDQoNCkRldGFpbGVkIGV4cGxhbmF0aW9ucyB5b3UgY2FuIHJlYWQgW2hlcmVdKGh0dHA6Ly9ycHVicy5jb20vY2hpZHVuZ2t0LzI5NzgyNSkuIEFuIGVtcGlyaWNhbCBldmlkZW5jZSB5b3UgY2FuIGZpbmQgW2hlcmVdKGh0dHA6Ly9ycHVicy5jb20vY2hpZHVuZ2t0LzQ4NzkxMikuIA0KDQojIE1heGltaXplIFRocmVzaG9sZCBmb3IgQWNjdXJhY3kgYW5kIFNlbnNpdGl2aXR5DQoNCg0KYGBge3J9DQojPT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09DQojIFVzaW5nIEFjY3VyYWN5IENyaXRlcmlvbiBmb3IgQ3JlZGl0IENsYXNzaWZpY2F0aW9uIENhbiBDYXVzZSBEaXNhc3RlciBmb3IgQmFua3MuIA0KIz09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PQ0KIA0KDQojIExvYWQgUiBwYWNrYWdlczogDQpybShsaXN0ID0gbHMoKSkNCmxpYnJhcnkoY3V0cG9pbnRyKQ0KbGlicmFyeSh0aWR5dmVyc2UpDQpsaWJyYXJ5KGNhcmV0KQ0KDQojIExvYWQgR2VybWFuQ3JlZGl0IERhdGE6IA0KZGF0YSgiR2VybWFuQ3JlZGl0IikNCmRmIDwtIEdlcm1hbkNyZWRpdCANCg0KIyBTcGxpdCBkYXRhOiANCnNldC5zZWVkKDEpDQppZCA8LSBjcmVhdGVEYXRhUGFydGl0aW9uKHkgPSBkZiRDbGFzcywgcCA9IDAuNywgbGlzdCA9IEZBTFNFKQ0KdHJhaW4gPC0gZGZbaWQsIF0NCnRlc3QgPC0gZGZbLWlkLCBdDQoNCiMgU2V0IGNvbmRpdGlvbnMgZm9yIHRyYWluaW5nIExvZ2lzdGljIE1vZGVsOiANCg0KbnVtYmVyIDwtIDUNCnJlcGVhdHMgPC0gNQ0KbiA8LSBudW1iZXIqcmVwZWF0cw0KDQpzZXQuc2VlZCgxKQ0KdHJhaW4uY29udHJvbCA8LSB0cmFpbkNvbnRyb2wobWV0aG9kID0gInJlcGVhdGVkY3YiLCANCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIG51bWJlciA9IG51bWJlciwNCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHJlcGVhdHMgPSByZXBlYXRzLCANCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGNsYXNzUHJvYnMgPSBUUlVFLA0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgYWxsb3dQYXJhbGxlbCA9IFRSVUUsIA0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgc3VtbWFyeUZ1bmN0aW9uID0gdHdvQ2xhc3NTdW1tYXJ5KQ0KDQoNCiMgVHJhbmkgTG9naXN0aWMgTW9kZWwgd2l0aCB0cmFpbiBkYXRhOiANCm15X2xvZ2l0IDwtIHRyYWluKENsYXNzIH4uLCANCiAgICAgICAgICAgICAgICAgIGRhdGEgPSB0cmFpbiwgDQogICAgICAgICAgICAgICAgICBtZXRob2QgPSAiZ2xtIiwgDQogICAgICAgICAgICAgICAgICB0ckNvbnRyb2wgPSB0cmFpbi5jb250cm9sKQ0KDQoNCiMgVXNlIG1vZGVsIGZvciBwcmVkaWN0aW5nIFBEOiANCnBkIDwtIHByZWRpY3QobXlfbG9naXQsIHRlc3QsIHR5cGUgPSAicHJvYiIpICU+JSBwdWxsKEJhZCkNCg0KDQojIENhbGN1bGF0ZSBvcHRpbWFsIGN1dG9mZiBieSBzZW5zaXRpdml0eSBhbmQgYWNjdXJhY3kgY3JpdGVyaW9uOiANCm0gPC0gY3V0cG9pbnRyKHggPSBwZCwgY2xhc3MgPSB0ZXN0JENsYXNzLCBtZXRyaWMgPSBzZW5zX2NvbnN0cmFpbiwgcG9zX2NsYXNzID0gIkJhZCIpDQpuIDwtIGN1dHBvaW50cih4ID0gcGQsIGNsYXNzID0gdGVzdCRDbGFzcywgbWV0cmljID0gYWNjdXJhY3ksIHBvc19jbGFzcyA9ICJCYWQiKQ0KDQojIE9wdGltYWwgY3V0b2ZmOiANCnQxIDwtIG0kb3B0aW1hbF9jdXRwb2ludCAjIG1heGltaXplIHNlbnNpdGl2aXR5LiANCnQyIDwtIG4kb3B0aW1hbF9jdXRwb2ludCAjIG1heGltaXplIGFjY3VyYWN5LiANCg0KIyBQbG90IFJPQyBjdXJ2ZSB3aXRoIG9wdGltYWwgY3V0b2ZmOiANCmdyaWRFeHRyYTo6Z3JpZC5hcnJhbmdlKHBsb3QobSksIHBsb3QobikpDQoNCiMgWCBpcyB1bmNoYW5nZWQgcmVnYXJkbGVzcyBvZiBtZXRob2Qgc2VsZWN0ZWQ6IA0KDQptJEFVQw0KbiRBVUMNCmBgYA0KDQoNCiMgRWZmZWN0cyBvbiBQcm9maXQNCg0KDQpgYGB7cn0NCiMgRnVuY3Rpb24gZm9yIGxhYmVsbGluZyBjcmVkaXQgYXBwbGljYXRpb25zOiANCg0KbGFiZWxfcHJlZGljdGVkIDwtIGZ1bmN0aW9uKGN1dG9mZikgew0KICB5IDwtIGNhc2Vfd2hlbihwZCA+PSBjdXRvZmYgfiAiQmFkIiwgVFJVRSB+ICJHb29kIikgJT4lIGFzLmZhY3RvcigpDQogIHJldHVybih5KQ0KfQ0KDQojIENvbmZ1dGlvbiBtYXh0cml4OiANCmNvbmZ1c2lvbk1hdHJpeChsYWJlbF9wcmVkaWN0ZWQodDEpLCB0ZXN0JENsYXNzLCBwb3NpdGl2ZSA9ICJCYWQiKSAjIElmIHQxIGlzIHNlbGVjdGVkIGZvciBjbGFzc2lmaWNhdGlvbi4gDQpjb25mdXNpb25NYXRyaXgobGFiZWxfcHJlZGljdGVkKHQyKSwgdGVzdCRDbGFzcywgcG9zaXRpdmUgPSAiQmFkIikgIyBJZiB0MiBpcyBzZWxlY3RlZCBmb3IgY2xhc3NpZmljYXRpb24uIA0KDQojIFJlc3VsdHMgb2YgQ2xhc3NpZmljYXRpb24gYnkgdGhlIHR3byBtZXRob2RzOiANCg0KdGVzdCAlPiUgDQogIG11dGF0ZShDbGFzc19tYXhfc2VuID0gbGFiZWxfcHJlZGljdGVkKHQxKSwNCiAgICAgICAgIENsYXNzX21heF9hY2N1cmFjeSA9IGxhYmVsX3ByZWRpY3RlZCh0MikpIC0+IHRlc3RfdDFfdDINCg0KdGVzdF90MV90MiAlPiUgDQogIGZpbHRlcihDbGFzc19tYXhfc2VuID09ICJHb29kIiwgQ2xhc3MgPT0gIkdvb2QiKSAlPiUgDQogIHNlbGVjdChDbGFzcywgQ2xhc3NfbWF4X3NlbiwgQW1vdW50KSAtPiBkZl9Hb29kX21heF9zZW4NCg0KdGVzdF90MV90MiAlPiUgDQogIGZpbHRlcihDbGFzc19tYXhfc2VuID09ICJHb29kIiwgQ2xhc3MgPT0gIkJhZCIpICU+JSANCiAgc2VsZWN0KENsYXNzLCBDbGFzc19tYXhfc2VuLCBBbW91bnQpIC0+IGRmX0JhZF9tYXhfc2VuDQoNCnRlc3RfdDFfdDIgJT4lIA0KICBmaWx0ZXIoQ2xhc3NfbWF4X2FjY3VyYWN5ID09ICJHb29kIiwgQ2xhc3MgPT0gIkdvb2QiKSAlPiUgDQogIHNlbGVjdChDbGFzcywgQ2xhc3NfbWF4X2FjY3VyYWN5LCBBbW91bnQpIC0+IGRmX0dvb2RfbWF4X2FjYw0KDQp0ZXN0X3QxX3QyICU+JSANCiAgZmlsdGVyKENsYXNzX21heF9hY2N1cmFjeSA9PSAiR29vZCIsIENsYXNzID09ICJCYWQiKSAlPiUgDQogIHNlbGVjdChDbGFzcywgQ2xhc3NfbWF4X2FjY3VyYWN5LCBBbW91bnQpIC0+IGRmX0JhZF9tYXhfYWNjDQoNCg0KIyBDYWxjdWxhdGUgcHJvZml0IGlmIGludGVyZXN0IHJhdGUgPSAzMCU6IA0KaXIgPC0gMC4zDQpzdW0oaXIqZGZfR29vZF9tYXhfc2VuJEFtb3VudCkgLSBzdW0oZGZfQmFkX21heF9zZW4kQW1vdW50KQ0Kc3VtKGlyKmRmX0dvb2RfbWF4X2FjYyRBbW91bnQpIC0gc3VtKGRmX0JhZF9tYXhfYWNjJEFtb3VudCkNCg0KDQoNCiMgRnVuY3Rpb24gZm9yIGNhbGN1bGF0aW5nIHByb2ZpdDogDQoNCnByb2ZpdCA8LSBmdW5jdGlvbihpcikgew0KICBwcm8xIDwtIHN1bShpcipkZl9Hb29kX21heF9zZW4kQW1vdW50KSAtIHN1bShkZl9CYWRfbWF4X3NlbiRBbW91bnQpDQogIHBybzIgPC0gc3VtKGlyKmRmX0dvb2RfbWF4X2FjYyRBbW91bnQpIC0gc3VtKGRmX0JhZF9tYXhfYWNjJEFtb3VudCkNCiAgcmV0dXJuKGRhdGEuZnJhbWUoSVIgPSBpciwgUHJvID0gYyhwcm8xLCBwcm8yKSwgTWV0aG9kID0gYygiTWF4U2VuIiwgIk1heEFjYyIpKSkNCn0NCg0KDQojIENvbXBhcmUgcHJvZml0OiANCg0KbGFwcGx5KHNlcSgwLjE1LCAwLjUsIGJ5ID0gMC4wMSksIHByb2ZpdCkgLT4gcHJvZml0X2xpc3QgDQoNCmRvLmNhbGwoInJiaW5kIiwgcHJvZml0X2xpc3QpICU+JSANCiAgbXV0YXRlKFBybyA9IFBybyAvIDEwMDApICU+JSANCiAgZ2dwbG90KGFlcyhJUiwgUHJvLCBjb2xvciA9IE1ldGhvZCkpICsgDQogIGdlb21fbGluZSgpICsgDQogIGdlb21fcG9pbnQoKSArIA0KICBzY2FsZV94X2NvbnRpbnVvdXMoYnJlYWtzID0gc2VxKDAuMTUsIDAuNSwgYnkgPSAwLjA1KSwgbGFiZWxzID0gc2NhbGVzOjpwZXJjZW50KSArIA0KICBsYWJzKHggPSAiSW50ZXJlc3QgUmF0ZSIsIHkgPSAiUHJvZml0IiwgDQogICAgICAgdGl0bGUgPSAiUHJvZml0IENvbXBhcmlzaW9uIGJ5IENsYXNzaWZpY2F0aW9uIE1ldGhvZCIpICsgDQogIGhyYnJ0aGVtZXM6OnRoZW1lX21vZGVybl9yYygpDQoNCmBgYA0KDQo=