Library

library(readxl) # read Excel File
library(dplyr) #data cleaning 
library(kableExtra) # Pretty Format
library(scorecard) # Scorecard builder
library(rsample) # split train test

Data Source

Dataset yang digunakan adalah data Taiwan Credit. sumber data: Source

Dataset memiliki 25 kolom dan 30,952 dengan penjelasan sebagai berikut:

• id = id debitur
• limit_bal = Besaran kredit limit yang diberikan dalam dolar NT
• sex = jenis kelamin
  • 1 = laki-laki
  • 2 = perempuan
• education = Pendidikan terakhir
  • 1 = pascasarjana (S2 & S3)
  • 2 = universitas (S1)
  • 3 = high school (SMA)
  • 4 = lain-lain
• marriage = Status pernikahan
  • 1 = menikah
  • 2 = lajang
  • 3 = lainnya
• age = Usia dalam tahun
• pay_* = Status pembayaran dalam bulan April (1) - September (6).
  • 0 = pembayaran tepat waktu
  • 1 = keterlambatan pembayaran satu bulan
  • 2 = keterlambatan pembayaran dua bulan
  • …
  • 8 = keterlambatan pembayaran delapan bulan atau lebih
• bill_amt* = Jumlah tagihan pada bulan April (1) - September (6) dalam dolar NT
• pay_amt* = Jumlah pembayaran/pengeluaran sebelumnya pada bulan April (1) - September(6) dalam dolar NT
• gb_flag = Flagging pembayaran default (gagal bayar) pada bulan berikutnya
  • 1 = default
  • 0 = not default

datasource <- read_excel("data_input/credit_taiwan.xlsx") %>% 
  select(-id) %>% 
  mutate(sex = as.factor(sex),
         education = as.factor(education),
         marriage = as.factor(marriage))
  

datasource %>% 
  head(2) %>% 
  kbl() %>%
  kable_styling()

Function

# EDA
create_eda_report <- function(data) {
  results <- list()
  
  for (column in colnames(data)) {
    column_name <- column
    column_type <- class(data[[column]])
    unique_count <- length(unique(data[[column]]))
    total <- nrow(data)
    null_count <- sum(is.na(data[[column]]))
    non_null_count <- total - null_count
    null_percentage <- sprintf("%.2f%%", (null_count / total * 100))
    
    if (column_type == "character" || column_type == "factor") {
      value_counts <- table(data[[column]])
      most_frequent <- names(sort(value_counts, decreasing = TRUE))[1]
      mode_percentage <- sprintf("%.2f%%", (sort(value_counts, decreasing = TRUE)[1] / total * 100))
      mean_or_top1 <- ""
      min_or_bottom <- names(sort(value_counts))[1]
      if (length(value_counts) >= 2) {
        min_or_bottom <- paste0(min_or_bottom, ": ", sprintf("%.2f%%", (sort(value_counts)[1] / total * 100)))
      }
      max_or_top1 <- names(sort(value_counts, decreasing = TRUE))[1]
      if (length(value_counts) >= 2) {
        max_or_top1 <- paste0(max_or_top1, ": ", sprintf("%.2f%%", (sort(value_counts, decreasing = TRUE)[1] / total * 100)))
      }
    } else {
      most_frequent <- names(sort(table(data[[column]]), decreasing = TRUE))[1]
      mode_percentage <- sprintf("%.2f%%", (sort(table(data[[column]]), decreasing = TRUE)[1] / total * 100))
      mean_or_top1 <- sprintf("%.3f", mean(data[[column]], na.rm = TRUE))
      min_or_bottom <- min(data[[column]], na.rm = TRUE)
      max_or_top1 <- max(data[[column]], na.rm = TRUE)
    }
    
    results[[length(results) + 1]] <- data.frame(
      column_name = column_name,
      type = column_type,
      unique_count = unique_count,
      total = total,
      null_count = null_count,
      non_null_count = non_null_count,
      null_percentage = null_percentage,
      most_frequent = most_frequent,
      mode_percentage = mode_percentage,
      mean = mean_or_top1,
      min_or_bottom = min_or_bottom,
      max_or_top1 = max_or_top1,
      stringsAsFactors = FALSE
    )
  }
  
  column_info <- do.call(rbind, results)
  
  return(column_info)
}

#Punctual & Overdue
count_consecutive_od <- function(row) {
  count <- 0
  max_count <- 0
  
  for (value in row) {
    if (value > 0) {
      count <- count + 1
      max_count <- max(max_count, count)
    } else {
      count <- 0
    }
  }
  
  return(max_count)
}

count_consecutive_punct <- function(row) {
  count <- 0
  max_count <- 0
  
  for (value in row) {
    if (value <= 0) {
      count <- count + 1
      max_count <- max(max_count, count)
    } else {
      count <- 0
    }
  }
  
  return(max_count)
}

# Approval Rate
approval_rate <- function(score, label, positive = 0){

  score_list <- list(data = score)
  label_list <- list(data = label)

  g <- gains_table(score = score_list, label = label_list, positive = positive)
  
  final_df <- g %>% 
    mutate(
      count_approved = max(cum_count) - cum_count,
      neg_approved = max(cum_neg) - cum_neg,
      neg_rate = round((neg_approved / count_approved), 4)
    ) %>% 
    replace(is.na(.), 0) %>% 
    select(bin, approval_rate, neg_rate, 
           count_approved, neg_approved,
           count, neg, pos)
  
  final_df
}

Derivative

Kali ini saya akan menambahkan beberapa variable derivative untuk menambah kombinasi variasi dari setiap parameter. Berikut adalah definisinya:

• der_sex_education = Penggabungan Jenis kelamin dan pendidikan
• der_sex_marriage = Penggabungan Jenis kelamin dan Status Perkawinan
• der_marriage_education = Penggabungan JStatus Perkawinan dan pendidikan
• der_max_pay = Max Overdue
• der_max_pay_first_3 = Max Overdue 3 bulan pertama
• der_max_pay_last_3 = Max Overdue 3 bulan terakhir
• der_cm_od_first_3 = berapa kali terlambat di 3 bulan pertama
• der_cm_od_last_3 = berapa kali terlambat di 3 bulan terakhir
• der_cm_od = berapa kali terlambat
• der_cm_punct_first_3 = berapa kali tepat waktu di 3 bulan pertama
• der_cm_punct_last_3 = berapa kali tepat waktu di 3 bulan terakhir
• der_cm_punct = berapa kali tepat waktu

datasource$der_sex_education <- paste(datasource$sex, datasource$education, sep = "-")
datasource$der_sex_marriage <- paste(datasource$sex, datasource$marriage, sep = "-")
datasource$der_marriage_education <- paste(datasource$marriage, datasource$education, sep = "-")
datasource$der_max_pay <- apply(datasource[, 6:11], 1, max)
datasource$der_max_pay_first_3 <- apply(datasource[, 6:9], 1, max)
datasource$der_max_pay_last_3 <- apply(datasource[, 9:11], 1, max)
datasource$der_cm_od_first_3 <- apply(datasource[, 6:8], 1, count_consecutive_od)
datasource$der_cm_od_last_3 <- apply(datasource[, 9:11], 1, count_consecutive_od)
datasource$der_cm_od <- apply(datasource[, 6:11], 1, count_consecutive_od)
datasource$der_cm_punct_first_3 <- apply(datasource[, 6:9], 1, count_consecutive_punct)
datasource$der_cm_punct_last_3 <- apply(datasource[, 9:11], 1, count_consecutive_punct)
datasource$der_cm_punct <- apply(datasource[, 6:11], 1, count_consecutive_punct)

Minor Data Cleaning

Untuk memastikan variable derifatif yang bersifat kategorik bertipe data factor

datasource <-datasource %>%  
  mutate(der_sex_education = as.factor(der_sex_education),
         der_sex_marriage = as.factor(der_sex_marriage),
         der_marriage_education = as.factor(der_marriage_education))

datasource %>% 
  head(2) %>% 
  kbl() %>%
  kable_styling()

EDA

datasource_eda <- create_eda_report(datasource) %>% 
  kbl() %>%
  kable_styling()

datasource_eda

Dari hasil EDA dapat disimpulkan bahwa dari dokumentasi variable terdapat null value yang di ganti dengan nilai 0 di marriage status yang akan diinterpretasikan sebagai debitur yang tidak menginput marriage status. Untuk Variable lain didapati 0% Missing value

Data Splitting

Metode splitting dengan mengambil 70% untuk training data dan 30% untuk validasi / Test data

RNGkind(sample.kind= "Rounding")
set.seed(52)

splitter <- initial_split(data = datasource,
                          prop = 0.7) 

train <- training(splitter)
test <- testing(splitter)

Binning

binning <- woebin(dt = train,
                  y = 'gb_flag',
                  positive = 0)

#> ✔ Binning on 21666 rows and 36 columns in 00:00:14

binning_df <- bind_rows(binning, .id = "variable")

binning_df %>% 
  kbl() %>%
  kable_styling()

Transform to WoE

Transformasi nilai asli dengan hasil WoE per masing masing bin variable untuk pemodelan scorecard

train_woe <- woebin_ply(dt = train,
                        bins = binning)

#> ✔ Woe transformating on 21666 rows and 35 columns in 00:00:12

test_woe <- woebin_ply(dt = test,
                        bins = binning)

#> ✔ Woe transformating on 9286 rows and 35 columns in 00:00:02

Information Value

Karena memang dikhususkan untuk existing customer variable derivatif terkait jumlah berapa kali pembayaran terlambat mendapat IV tertinggi, kemungkinan dikarenakan dalam penentuan target variable diambil dari riwayat pembiayaan.

info_val <- iv(dt = train_woe,
   y ='gb_flag',
   positive = 0)



info_val

#>                       variable    info_value
#>  1:              der_cm_od_woe 2.52036183299
#>  2:           der_cm_punct_woe 2.44561980495
#>  3:   der_cm_punct_first_3_woe 1.61550587939
#>  4:      der_cm_od_first_3_woe 1.40762584814
#>  5:    der_cm_punct_last_3_woe 1.30669920771
#>  6:            der_max_pay_woe 1.24363616181
#>  7:       der_cm_od_last_3_woe 1.21682171977
#>  8:    der_max_pay_first_3_woe 1.14345582141
#>  9:     der_max_pay_last_3_woe 0.95236305086
#> 10:                  pay_1_woe 0.90156463822
#> 11:                  pay_4_woe 0.89215545424
#> 12:                  pay_5_woe 0.86446845529
#> 13:                  pay_2_woe 0.74473009488
#> 14:                  pay_3_woe 0.73165916222
#> 15:                  pay_6_woe 0.71573380240
#> 16:              limit_bal_woe 0.68600256213
#> 17:              bill_amt6_woe 0.36505153390
#> 18:              bill_amt4_woe 0.34552057008
#> 19:              bill_amt5_woe 0.34278407367
#> 20:              bill_amt3_woe 0.30144297443
#> 21:              bill_amt2_woe 0.27862402949
#> 22:              bill_amt1_woe 0.24178320934
#> 23:      der_sex_education_woe 0.10661170429
#> 24:               pay_amt2_woe 0.09651588574
#> 25: der_marriage_education_woe 0.09379099179
#> 26:               pay_amt1_woe 0.09246700895
#> 27:              education_woe 0.07869044349
#> 28:               pay_amt6_woe 0.06885583223
#> 29:               pay_amt3_woe 0.06603215570
#> 30:               pay_amt5_woe 0.05397319591
#> 31:                    age_woe 0.04981323815
#> 32:               pay_amt4_woe 0.04418106082
#> 33:       der_sex_marriage_woe 0.01720165935
#> 34:                    sex_woe 0.01646469643
#> 35:               marriage_woe 0.00006328743
#>                       variable    info_value

Feature selection

Setelah di transform ke WoE maka dilakukan feature selection dengan mengeliminasi variable yang saling berkorelasi agar tidak terkena double punishment dan apabila saling berkorelasi maka dipilih yang IV nya lebih tinggi

cor_matrix <- cor(train_woe %>% select(-gb_flag))
rounded_abs_cor_matrix <- round(abs(cor_matrix), 2)

formatted_matrix <- rounded_abs_cor_matrix %>%
  as.data.frame() %>%
  mutate_all(~ cell_spec(., "html", color = "white", background = ifelse(. > 0.4, "red", "white")))


formatted_matrix %>%
  kbl(escape = FALSE) %>%
  kable_styling()

Binning Plot

Mari highlight variable der_cm_od_woe, dapat dilihat bahwa semakin sering seseorang tepat waktu maka akan semakin kecil probability of default nya. Antar Bin sudah risk rank dan merepresentasikan bin yang baik. Hasil IV bisa tinggi kemungkinan karena dalam penentuan badrate dilakukan analisa seperti Roll rate analysis yang diambil dari riwayat pembiayaan juga

plot <- woebin_plot(bins = binning)
plot$der_cm_od

Final Selection

pada tahap akhir setelah mendapatkan nilai IV tertinggi dan variable yang tidak saling berkorelasi, didapatkan variable berikut: der_cm_od_woe,limit_bal_woe,bill_amt6_woe,der_sex_education_woe,pay_amt2_woe,age_woe

train_woe_final<- train_woe %>% 
  select(c(der_cm_od_woe,limit_bal_woe,bill_amt6_woe,der_sex_education_woe,pay_amt2_woe,age_woe,gb_flag))

test_woe_final<- test_woe %>% 
  select(c(der_cm_od_woe,limit_bal_woe,bill_amt6_woe,der_sex_education_woe,pay_amt2_woe,age_woe,gb_flag))

train_woe_final %>% 
  head(2) %>% 
  kbl() %>%
  kable_styling()

Model Building

Semua variable yang digunakan dalam pembuatan model signifikan terhadap target kecuali variable umur

model <- glm(formula = gb_flag ~ .,
             data = train_woe_final,
             family = "binomial")

summary(model) 

#> 
#> Call:
#> glm(formula = gb_flag ~ ., family = "binomial", data = train_woe_final)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -3.9496  -0.6513  -0.2721   0.6772   2.4296  
#> 
#> Coefficients:
#>                       Estimate Std. Error z value             Pr(>|z|)    
#> (Intercept)           -0.15801    0.01857  -8.507 < 0.0000000000000002 ***
#> der_cm_od_woe         -0.86343    0.01793 -48.152 < 0.0000000000000002 ***
#> limit_bal_woe         -0.84729    0.02421 -35.000 < 0.0000000000000002 ***
#> bill_amt6_woe         -1.11035    0.03081 -36.037 < 0.0000000000000002 ***
#> der_sex_education_woe -0.32730    0.05701  -5.741   0.0000000094126076 ***
#> pay_amt2_woe          -0.47269    0.06219  -7.601   0.0000000000000293 ***
#> age_woe               -0.08688    0.08181  -1.062                0.288    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 29911  on 21665  degrees of freedom
#> Residual deviance: 18937  on 21659  degrees of freedom
#> AIC: 18951
#> 
#> Number of Fisher Scoring iterations: 6

Vif

Untuk membuktikan korelasi antar variable yang digunakan dalam model

vif(model)

#>                 variable     gvif
#> 1:         der_cm_od_woe 1.025403
#> 2:         limit_bal_woe 1.213765
#> 3:         bill_amt6_woe 1.092026
#> 4: der_sex_education_woe 1.053359
#> 5:          pay_amt2_woe 1.037778
#> 6:               age_woe 1.093182

Prediction

test_woe_final$pred_risk <- predict(object = model,
                                    newdata = test_woe_final,
                                    type = "response")
test_woe_final %>% 
  head(2) %>% 
  kbl() %>%
  kable_styling()

Model Performance

Jika dilihat nilai KS sebesar 0,57 yang berarti model mampu membedakan kelas negatif dan positif dengan jarak 57% (cukup besar kemampuan nya dalam membedakan kelas positif dan negatif, sehingga meminimalisir prediksi yang salah)

list_pred <- list(test = test_woe_final$pred_risk)
list_label <- list(test = test_woe_final$gb_flag)

perf_eva(pred = list_pred,
         label = list_label,
         confusion_matrix = TRUE,
         threshold = 0.5,
         show_plot = c("ks", "roc"))

#> $binomial_metric
#> $binomial_metric$test
#>          MSE      RMSE   LogLoss        R2        KS       AUC      Gini
#> 1: 0.1466613 0.3829638 0.4400123 0.4092079 0.5712766 0.8743167 0.7486334
#> 
#> 
#> $confusion_matrix
#> $confusion_matrix$test
#>    label pred_0 pred_1     error
#> 1:     0   4242    790 0.1569952
#> 2:     1   1237   3017 0.2907851
#> 3: total   5479   3807 0.2182856
#> 
#> 
#> $pic
#> TableGrob (1 x 2) "arrange": 2 grobs
#>   z     cells    name           grob
#> 1 1 (1-1,1-1) arrange gtable[layout]
#> 2 2 (1-1,2-2) arrange gtable[layout]

Pembuatan Scorecard

score_card <- scorecard(bins = binning, 
                        model = model,
                        odds0 = 1/19,
                        points0 = 600,
                        pdo = 20)

scorecard_df <- bind_rows(score_card, .id = "variable")

scorecard_df %>% 
  kbl() %>%
  kable_styling()

score_train <- scorecard_ply(dt = train,
                       card = score_card,
                       only_total_score = F)

score_test <- scorecard_ply(dt = test,
                       card = score_card,
                       only_total_score = F)

score_train %>% head() %>% 
  kbl() %>%
  kable_styling()

PSI

Dapat dismpulkan bahwa PSI antara train dan test sebesar 0.001937147 yang berarti tidak ada pergeseran signifikan antar variable yang digunakan dalam pembuatan scoring sehingga dapat dikatakan variable di data train dan test stabil

score_list <- list(train = score_train$score, 
                   test = score_test$score)

label_list <- list(train = train_woe_final$gb_flag, 
                   test = test_woe_final$gb_flag)

psi <- perf_psi(score = score_list, 
                label = label_list, 
                positive = 0)
psi$psi # psi data frame

#>    variable    dataset         psi
#> 1:     pred train_test 0.001937147

head(score_test,2) %>% 
    kbl() %>%
  kable_styling()

Approval Rate

approval_rate_table <- approval_rate(score = score_test$score,
              label = test_woe_final$gb_flag,
              positive = 0)

approval_rate_table %>% 
  kbl() %>%
  kable_styling()

Kesimpulan

Cut-off score yang dapat diambil >=426 dengan menghasilkan approval rate sekitar 80% dan menyisakan badrate sekitar 34%