STA581-05 - Logistics Regression

A. Pendahuluan

A.1 Package

Berikut adalah packages yang akan digunakan dalam praktikum ini, yaitu :

Nama Package	Kegunaan
`broom`	untuk merapikan tampilan data
`caret`	akan digunakan beberapa fungsi untuk membentuk Confusion Matrix
`DataExplorer`	akan digunakan beberapa fungsi seperti `plot_bar`, `plot_histogram`, `plot_boxplot`, untuk membantu visualisasi data
`grid`	diperlukan agar user defined function `arrange` bisa digunakan
`InformationValue`	akan digunakan fungsi `plotROC` untuk menampilkan Receiver Operating Characteristics (ROC) curve
`ISLR`	Sumber Datasets `default`
`pscl`	akan digunakan fungsi `pR2` (Pseudo $R^2$) untuk Evaluasi Model
`tidyverse`

Untuk melakukan Instalasi, Anda bisa menggunakan syntax berikut ini :

install.packages(c("broom","caret","DataExplorer","grid","InformationValue","ISLR","pscl","tidyverse"))

Jangan lupa untuk load package tersebut.

library(broom)
library(caret)
library(DataExplorer)
library(grid)
library(InformationValue)
library(ISLR)
library(pscl)
library(tidyverse)

Sebenarnya Anda tidak perlu untuk load keseluruhan pakcage, Anda bisa menggunakan syntax namapackage::fungsi untuk menggunakan suatu fungsi dari suatu package di R, selama package tersebut sudah Anda Install.

A.2 Data

Data yang akan digunakan dalam praktikum ini ada 2, yaitu :
- Default Datasets yang berasal dari package ISLR
- Lending Club Datasets

A.3 Tambahan

Untuk mempermudah visualisasi beberapa bagian dalam Praktikum ini, digunakan fungsi arrange yang dibuat oleh Stephen Turner

vp.layout <- function(x, y) viewport(layout.pos.row=x, layout.pos.col=y)
arrange <- function(..., nrow=NULL, ncol=NULL, as.table=FALSE) {
  dots <- list(...)
  n <- length(dots)
  if(is.null(nrow) & is.null(ncol)) { nrow = floor(n/2) ; ncol = ceiling(n/nrow)}
  if(is.null(nrow)) { nrow = ceiling(n/ncol)}
  if(is.null(ncol)) { ncol = ceiling(n/nrow)}
  ## NOTE see n2mfrow in grDevices for possible alternative
  grid.newpage()
  pushViewport(viewport(layout=grid.layout(nrow,ncol) ) )
  ii.p <- 1
  for(ii.row in seq(1, nrow)){
    ii.table.row <- ii.row
    if(as.table) {ii.table.row <- nrow - ii.table.row + 1}
    for(ii.col in seq(1, ncol)){
      ii.table <- ii.p
      if(ii.p > n) break
      print(dots[[ii.table]], vp=vp.layout(ii.table.row, ii.col))
      ii.p <- ii.p + 1
    }
  }
}

B. `Default` Datasets

B.1 Eksplorasi Data

data(Default)

Untuk mengetahui lebih lanjut mengenai Default, Anda bisa menggunakan syntax berikut :

help(Default)

Deskripsi

A simulated data set containing information on ten thousand customers. The aim here is to predict which customers will default on their credit card debt.

Format

A data frame with 10000 observations on the following 4 variables.

default : A factor with levels No and Yes indicating whether the customer defaulted on their debt
student : A factor with levels No and Yes indicating whether the customer is a student
balance : The average balance that the customer has remaining on their credit card after making their monthly payment
income : Income of customer

glimpse

Anda bisa melihatnya dengan glimpse dari dplyr yang merupakan bagian dari tidyverse

glimpse(Default)

## Rows: 10,000
## Columns: 4
## $ default <fct> No, No, No, No, No, No, No, No, No, No, No, No, No, No, No,...
## $ student <fct> No, Yes, No, No, No, Yes, No, Yes, No, No, Yes, Yes, No, No...
## $ balance <dbl> 729.5265, 817.1804, 1073.5492, 529.2506, 785.6559, 919.5885...
## $ income  <dbl> 44361.625, 12106.135, 31767.139, 35704.494, 38463.496, 7491...

summary

Seandainya kita ingin mengetahui sebaran dan kecenderungan pola hubungan antar peubah, salah satu yang dapat dilakukan adalah sebagai berikut.

summary(Default)

##  default    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

plot_histogram

plot_histogram(Default)

plot_histogram ini menampilkan peubah dengan tipe numerik.

plot_bar

plot_bar(Default)

plot_bar ini menampilkan peubah dengan tipe kategorik.

plot_boxplot

Untuk melihat pola hubungan antara peubah numerik dan kategorik dapat lebih mudah diidentifikasi menggunakan plot_boxplot.

plot_boxplot(Default,by="default")

Output di atas memperlihatkan bahwa median pendapatan antara nasabah yang mengalami kredit macet atau tidak macet cederung mirip. Sebaliknya, terdapat perbedaan jumlah tagihan kartu kredit pada nasabah dengan kredit macet dan yang tidak. Pada kasus ini, orang yang memiliki tagihan lebih banyak, cenderung lebih mengalami kredit macet.

student - balance

a1<-Default %>% filter(student=="Yes") %>%
  select(balance,default) %>%  
  plot_boxplot(by="default", title="Default Occurrence on Student")

a2<-Default %>% filter(student=="No") %>%
  select(balance,default) %>%  
  plot_boxplot(by="default", title="Default Occurrence on Non-Student")

arrange(a1,a2)

## $page_1
## 
## $page_1

student - income

b1<-Default %>% filter(student=="Yes") %>%
  select(income,default) %>%  
  plot_boxplot(by="default", title="Default Occurrence on Student")

b2<-Default %>% filter(student=="No") %>%
  select(income,default) %>%  
  plot_boxplot(by="default", title="Default Occurrence on Non-Student")

arrange(b1,b2)

## $page_1
## 
## $page_1

B.2 The Model, Estimation and Inference

Sebagai awal dari ilustrasi yang akan dibahas pada modul kali ini, akan diperlihatkan pemodelan regresi logistik dengan menggunakan satu peubah penjelas terlebih dulu.

reglog<-glm(default~balance,data=Default,family=binomial)
summary(reglog)

## 
## Call:
## glm(formula = default ~ balance, family = binomial, data = Default)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2697  -0.1465  -0.0589  -0.0221   3.7589  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.065e+01  3.612e-01  -29.49   <2e-16 ***
## balance      5.499e-03  2.204e-04   24.95   <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: 2920.6  on 9999  degrees of freedom
## Residual deviance: 1596.5  on 9998  degrees of freedom
## AIC: 1600.5
## 
## Number of Fisher Scoring iterations: 8

Default %>%
  mutate(prob = ifelse(default == "Yes", 1, 0)) %>%
  ggplot(aes(balance, prob)) +
  geom_point(alpha = .15) +
  geom_smooth(formula=y~x, method = "glm", method.args = list(family = "binomial")) +
  ggtitle("Logistic regression model fit") +
  xlab("Balance") +
  ylab("Probability of Default")+
  theme_minimal()

B.2.1 Interpreting coefficients

reglog$coefficients

##   (Intercept)       balance 
## -10.651330614   0.005498917

Seperti yang telah dibahas pada kelas perkuliahan,

hasil $β_1$=0.0055 (positif) berarti membesarnya tagihan akan meningkatkan peluang gagal bayar. Pengguna kartu kredit yang lebih besar tagihannya memiliki risiko 0.0055 kali lebih besar untuk gagal bayar. Lebih tepat lagi, peningkatan satu unit tagihan (dollar) berhubungan dengan meningkatnya log-odds gagal bayar sebesar 0.0055 unit.

Untuk memudahkan interpretasi, seringkali digunakan rasio odds, yaitu rasio antara dua odds.

exp(reglog$coefficients)

##  (Intercept)      balance 
## 2.366933e-05 1.005514e+00

Pada ilustrasi ini, terlihat bahwa nilai rasio odds untuk peubah balance adalah 1.00551. Hal ini dapat diartikan bahwa ketika tagihan kartu kredit meningkat $ 1 maka kemungkinan kejadian kredit macet meningkat sebesar 1.00551 kali dibadingkan tagihannya tetap

confint(reglog)

## Waiting for profiling to be done...

##                     2.5 %       97.5 %
## (Intercept) -11.383288936 -9.966565064
## balance       0.005078926  0.005943365

B.2.2 Prediction

Seandainya kita ingin memprediksi peluang kredit macet jika besarnya tagihan seorang nasabah adalah sebesar $ 1,000

predict(reglog, newdata = data.frame(balance=1000), type="response")

##           1 
## 0.005752145

Artinya, sesorang yang memiliki jumlah tagihan $ 1,000 diprediksi memiliki peluang kredit macet sebesar 0.005752145.

B.2.3 Categorical Independent Variable

reglog2<-glm(default~student,data=Default,family=binomial)
summary(reglog2)

## 
## Call:
## glm(formula = default ~ student, family = binomial, data = Default)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.2970  -0.2970  -0.2434  -0.2434   2.6585  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.50413    0.07071  -49.55  < 2e-16 ***
## studentYes   0.40489    0.11502    3.52 0.000431 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2920.6  on 9999  degrees of freedom
## Residual deviance: 2908.7  on 9998  degrees of freedom
## AIC: 2912.7
## 
## Number of Fisher Scoring iterations: 6

predict(reglog2, newdata = data.frame(student="Yes"), type="res")

##          1 
## 0.04313859

predict(reglog2, newdata = data.frame(student="No"), type="res")

##          1 
## 0.02919501

B.2.4 Multiple Logistic Regression

full.mod = glm(default~., data=Default,family=binomial)
summary(full.mod)

## 
## Call:
## glm(formula = default ~ ., family = binomial, data = Default)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4691  -0.1418  -0.0557  -0.0203   3.7383  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.087e+01  4.923e-01 -22.080  < 2e-16 ***
## studentYes  -6.468e-01  2.363e-01  -2.738  0.00619 ** 
## balance      5.737e-03  2.319e-04  24.738  < 2e-16 ***
## income       3.033e-06  8.203e-06   0.370  0.71152    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2920.6  on 9999  degrees of freedom
## Residual deviance: 1571.5  on 9996  degrees of freedom
## AIC: 1579.5
## 
## Number of Fisher Scoring iterations: 8

Perhatikan bahwa pola pengaruh peubah penjelas tidak sama seperti yang dihasilkan pada model sebelumnya yang hanya menggunakan satu peubah penjelas. Dimana sebelumnya studentYes bernilai positif.

B.3 Model Evaluation

B.3.1 Pseudo $R^2$

pscl::pR2(reglog)["McFadden"]

## fitting null model for pseudo-r2

##  McFadden 
## 0.4533916

pscl::pR2(reglog2)["McFadden"]

## fitting null model for pseudo-r2

##    McFadden 
## 0.004097255

pscl::pR2(full.mod)["McFadden"]

## fitting null model for pseudo-r2

##  McFadden 
## 0.4619194

C. `Lending Club` Datasets

C.1 Eksplorasi Data

Pada ilustrasi ini akan digunakan data Lending Club Data.

loan <- read_csv("https://www.dropbox.com/s/89g1yyhwpcqwjn9/lending_club_cleaned.csv?raw=1")

## 
## -- Column specification --------------------------------------------------------
## cols(
##   good = col_character(),
##   purpose = col_character(),
##   fico = col_double(),
##   dti = col_double(),
##   loan_amnt = col_double(),
##   income = col_double()
## )

Peubah good merupakan indikator apakah pinjaman tersebut termasuk baik (good) atau buruk (bad);
Peubah fico merupakan skor kredit dimana skor yang lebih tinggi adalah yang lebih baik;
Peubah dti menunjukkan rasio debt-to-income,
Peubah loan_amnt menunjukkan jumlah pinjaman dan
Peubah income menunjukkan pendapatan. Sebagai catatan, peubah income bernilai NA jika sumber pendapatan tidak terverifikasi.

loan

## # A tibble: 42,535 x 6
##    good  purpose             fico   dti loan_amnt income
##    <chr> <chr>              <dbl> <dbl>     <dbl>  <dbl>
##  1 good  debt_consolidation   737 27.6       5000  24000
##  2 bad   major_purchase       742  1         2500  30000
##  3 good  small_business       737  8.72      2400     NA
##  4 good  other                692 20        10000  49200
##  5 good  other                697 17.9       3000  80000
##  6 good  vacation_wedding     732 11.2       5000  36000
##  7 good  debt_consolidation   692 23.5       7000     NA
##  8 good  major_purchase       662  5.35      3000  48000
##  9 bad   small_business       677  5.55      5600  40000
## 10 bad   other                727 18.1       5375  15000
## # ... with 42,525 more rows

summary(loan)

##      good             purpose               fico            dti       
##  Length:42535       Length:42535       Min.   :612.0   Min.   : 0.00  
##  Class :character   Class :character   1st Qu.:687.0   1st Qu.: 8.20  
##  Mode  :character   Mode  :character   Median :712.0   Median :13.47  
##                                        Mean   :715.1   Mean   :13.37  
##                                        3rd Qu.:742.0   3rd Qu.:18.68  
##                                        Max.   :827.0   Max.   :29.99  
##                                                                       
##    loan_amnt         income       
##  Min.   :  500   Min.   :   4800  
##  1st Qu.: 5200   1st Qu.:  44995  
##  Median : 9700   Median :  63000  
##  Mean   :11090   Mean   :  75186  
##  3rd Qu.:15000   3rd Qu.:  90000  
##  Max.   :35000   Max.   :6000000  
##                  NA's   :18758

table(loan$good)

## 
##   bad  good 
##  6371 36164

loan$good<-as.factor(loan$good)
loan

## # A tibble: 42,535 x 6
##    good  purpose             fico   dti loan_amnt income
##    <fct> <chr>              <dbl> <dbl>     <dbl>  <dbl>
##  1 good  debt_consolidation   737 27.6       5000  24000
##  2 bad   major_purchase       742  1         2500  30000
##  3 good  small_business       737  8.72      2400     NA
##  4 good  other                692 20        10000  49200
##  5 good  other                697 17.9       3000  80000
##  6 good  vacation_wedding     732 11.2       5000  36000
##  7 good  debt_consolidation   692 23.5       7000     NA
##  8 good  major_purchase       662  5.35      3000  48000
##  9 bad   small_business       677  5.55      5600  40000
## 10 bad   other                727 18.1       5375  15000
## # ... with 42,525 more rows

C.2 Membagi Data

Untuk kepentingan prediksi model, data akan dibagi dua menjadi data train dan data test. Karena pembagian data ini menggunakan sampling, maka perlu dilakukan pengaturan seed untuk randomisasi.

set.seed(5815)
sample <- sample(nrow(loan),floor(nrow(loan)*0.8))
train <- loan[sample,]
test <- loan[-sample,]

Sehingga, kita mempunyai dua data, train untuk membentuk model, dan test untuk menguji model yang telah dibuat.

C.3 Model

Pada ilustrasi ini, peubah fico, dti, loan_amnt, dan purpose akan digunakan sebagai peubah penjelas

logit <- glm(good ~ fico + dti+ loan_amnt + purpose, data = train, family = "binomial")
summary(logit)

## 
## Call:
## glm(formula = good ~ fico + dti + loan_amnt + purpose, family = "binomial", 
##     data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.7320   0.3863   0.5134   0.6182   1.2736  
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -7.209e+00  3.483e-01 -20.701  < 2e-16 ***
## fico                     1.336e-02  4.929e-04  27.098  < 2e-16 ***
## dti                     -1.007e-02  2.420e-03  -4.161 3.17e-05 ***
## loan_amnt               -2.353e-05  2.131e-06 -11.042  < 2e-16 ***
## purposeeducational      -5.343e-01  1.404e-01  -3.805 0.000142 ***
## purposehome_improvement -1.807e-01  5.879e-02  -3.075 0.002108 ** 
## purposemajor_purchase    5.330e-02  6.384e-02   0.835 0.403748    
## purposemedical          -4.120e-01  1.116e-01  -3.693 0.000222 ***
## purposeother            -3.539e-01  4.797e-02  -7.377 1.62e-13 ***
## purposesmall_business   -9.250e-01  6.222e-02 -14.867  < 2e-16 ***
## purposevacation_wedding  8.354e-02  9.607e-02   0.870 0.384502    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 28698  on 34027  degrees of freedom
## Residual deviance: 27547  on 34017  degrees of freedom
## AIC: 27569
## 
## Number of Fisher Scoring iterations: 5

C.4 Evaluasi Model

C.4.1 Confusion Matrix

Sensitivitas (atau True Positive Rate) adalah persentase pengamatan (aktual) yang diprediksi dengan benar oleh model, sedangkan spesifisitas adalah persentase dari 0 (aktual) yang diprediksi dengan benar.

test$pred <- predict(logit, test, type="response")

test$good_pred <- ifelse(test$pred > 0.80, "good", "bad")
test$good_pred <- as.factor(test$good_pred)
(conf.mat<-caret::confusionMatrix(test$good_pred, test$good, positive="good"))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  bad good
##       bad   418 1272
##       good  869 5948
##                                          
##                Accuracy : 0.7483         
##                  95% CI : (0.739, 0.7575)
##     No Information Rate : 0.8487         
##     P-Value [Acc > NIR] : 1              
##                                          
##                   Kappa : 0.1317         
##                                          
##  Mcnemar's Test P-Value : <2e-16         
##                                          
##             Sensitivity : 0.8238         
##             Specificity : 0.3248         
##          Pos Pred Value : 0.8725         
##          Neg Pred Value : 0.2473         
##              Prevalence : 0.8487         
##          Detection Rate : 0.6992         
##    Detection Prevalence : 0.8013         
##       Balanced Accuracy : 0.5743         
##                                          
##        'Positive' Class : good           
##

broom::tidy(conf.mat)

## # A tibble: 14 x 6
##    term                 class estimate conf.low conf.high   p.value
##    <chr>                <chr>    <dbl>    <dbl>     <dbl>     <dbl>
##  1 accuracy             <NA>     0.748    0.739     0.758  1.00e+ 0
##  2 kappa                <NA>     0.132   NA        NA     NA       
##  3 mcnemar              <NA>    NA       NA        NA      3.69e-18
##  4 sensitivity          good     0.824   NA        NA     NA       
##  5 specificity          good     0.325   NA        NA     NA       
##  6 pos_pred_value       good     0.873   NA        NA     NA       
##  7 neg_pred_value       good     0.247   NA        NA     NA       
##  8 precision            good     0.873   NA        NA     NA       
##  9 recall               good     0.824   NA        NA     NA       
## 10 f1                   good     0.847   NA        NA     NA       
## 11 prevalence           good     0.849   NA        NA     NA       
## 12 detection_rate       good     0.699   NA        NA     NA       
## 13 detection_prevalence good     0.801   NA        NA     NA       
## 14 balanced_accuracy    good     0.574   NA        NA     NA

C.4.2 ROC Curve

Receiver Operating Characteristics (ROC) curve melacak persentase true positive saat cut-off peluang prediksi diturunkan dari 1 menjadi 0. Model yang baik akan memperlihat kurva yang lebih curam, artinya True Positive Rate meningkat lebih cepat dibandingkan dengan False Positive Rate ketika cut-off menurun. Dengan kata lain, semakin besar luas area di bawah kurva ROC maka kemampuan prediksi yang dihasilkan oleh model semakin baik.

plotROC(test$good=="good", test$pred)

Referensi

Air Force Institute of Technology. (n.d.). Logistic regression · AFIT data science lab R programming guide. AFIT Data Science Lab R Programming Guide. Retrieved from: https://afit-r.github.io/logistic_regression
James, G., Witten, D., Hastie, T., Tibshirani, R. (2014). An Introduction to Statistical Learning with Applications in R. Springer.
Prabhakaran, S. (n.d.). Logistic regression with R. Tutorials on Advanced Stats and Machine Learning With R. Retrieved from: https://r-statistics.co/Logistic-Regression-With-R.html
Tomas, D. (n.d.). Lab 9 logistic regression. RPubs. Retrieved from: https://rpubs.com/dvorakt/151334
Turner S. 2010. Arrange multiple ggplot2 plots in the same image window [3 Maret 2021]. Retrieved from: https://gettinggeneticsdone.blogspot.com/2010/03/arrange-multiple-ggplot2-plots-in-same.html
Zhu X. (n.d.). Logistic regression, prediction and ROC. Xiaorui (Jeremy) Zhu 朱浚铭. Retrieved from: https://xiaoruizhu.github.io/Data-Mining-R/lecture/4_LogisticReg.html

IPB University, rahmaanisa@apps.ipb.ac.id︎↩︎
IPB University, gdito@apps.ipb.ac.id ↩︎
Badan Informasi Geospasial, abdul.aziz@big.go.id ↩︎

STA581-05 - Logistics Regression