Data Exploration

Sebagai ilustrasi, akan digunakan data Default yang berisi data simulasi sebanyak 10,000 konsumen kartu kredit.

library(ISLR)
library(tidyverse)
library(DataExplorer)

data(Default)
help(Default)

## starting httpd help server ... done

glimpse(Default)

## Rows: 10,000
## Columns: 4
## $ default <fct> No, 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, N…
## $ balance <dbl> 729.5265, 817.1804, 1073.5492, 529.2506, 785.6559, 919.5885, 8…
## $ income  <dbl> 44361.625, 12106.135, 31767.139, 35704.494, 38463.496, 7491.55…

Terlihat pada output di atas bahwa terdapat 4 peubah pada data tersebut, yaitu default, student, balance, dan income. Peubah default dan student merupakan peubah kategorik yang terdiri dari dua kategori, yaitu Yes atau No. Selanjutnya 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(Default)

plot_bar(Default)

Perhatikan bahwa kecenderungan pola hubungan antara peubah numerik dan kategorik dapat lebih mudah diidentifikasi menggunakan boxplot seperti yang dilakukan pada syntax berikut ini.

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.

tab = table (Default$student, Default$default)
tab[,2]/apply (tab,1,sum)

##         No        Yes 
## 0.02919501 0.04313859

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

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

Perhatikan bahwa setelah kita menelaah hal yang sama dengan mempertimbangkan apakah nasabah tersebut merupakan siswa atau bukan, apakah informasinya masih sama?

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(method = "glm", method.args = list(family = "binomial")) +
  ggtitle("Logistic regression model fit") +
  xlab("Balance") +
  ylab("Probability of Default")

## `geom_smooth()` using formula = 'y ~ x'

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

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.

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

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.

Model Evaluation

Pseudo $R^2$

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

## fitting null model for pseudo-r2

##  McFadden 
## 0.4533916

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

## fitting null model for pseudo-r2

##  McFadden 
## 0.4619194

Residual Assessment

plot(full.mod, which = 4, id.n = 5)

Testing logistic model out of sample

Pada ilustrasi ini akan digunakan data Lending Club Data. 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, dan income menunjukkan pendapatan. Sebagai catatan, peubah income bernilai NA jika sumber pendapatan tidak terverifikasi.

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

## Rows: 42535 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): good, purpose
## dbl (4): fico, dti, loan_amnt, income
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

loan

## # A tibble: 42,535 × 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
## # ℹ 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

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

table(loan$good)

## 
##   bad  good 
##  6371 36164

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


loan %>%
  rename()

## # A tibble: 42,535 × 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
## # ℹ 42,525 more rows

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

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

## 
## Call:  glm(formula = good ~ fico + dti + loan_amnt + purpose, family = "binomial", 
##     data = train)
## 
## Coefficients:
##             (Intercept)                     fico                      dti  
##              -7.206e+00                1.333e-02               -9.424e-03  
##               loan_amnt       purposeeducational  purposehome_improvement  
##              -2.301e-05               -4.904e-01               -1.662e-01  
##   purposemajor_purchase           purposemedical             purposeother  
##               3.395e-02               -3.516e-01               -3.974e-01  
##   purposesmall_business  purposevacation_wedding  
##              -1.004e+00                8.197e-02  
## 
## Degrees of Freedom: 34027 Total (i.e. Null);  34017 Residual
## Null Deviance:       28830 
## Residual Deviance: 27640     AIC: 27660

Model Evaluation

VIF

car::vif(logit4)

##               GVIF Df GVIF^(1/(2*Df))
## fico      1.070237  1        1.034523
## dti       1.062019  1        1.030543
## loan_amnt 1.130559  1        1.063278
## purpose   1.154023  7        1.010285

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. Spesifisitas juga dapat dihitung sebagai 1 - Sensitivitas.

library(caret)
test$pred <- predict(logit4, test, type="response")

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

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  bad good
##       bad   380 1270
##       good  869 5988
##                                           
##                Accuracy : 0.7486          
##                  95% CI : (0.7392, 0.7578)
##     No Information Rate : 0.8532          
##     P-Value [Acc > NIR] : 1               
##                                           
##                   Kappa : 0.1141          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.8250          
##             Specificity : 0.3042          
##          Pos Pred Value : 0.8733          
##          Neg Pred Value : 0.2303          
##              Prevalence : 0.8532          
##          Detection Rate : 0.7039          
##    Detection Prevalence : 0.8060          
##       Balanced Accuracy : 0.5646          
##                                           
##        'Positive' Class : good            
##

broom::tidy(conf.mat)

## # A tibble: 14 × 6
##    term                 class estimate conf.low conf.high   p.value
##    <chr>                <chr>    <dbl>    <dbl>     <dbl>     <dbl>
##  1 accuracy             <NA>     0.749    0.739     0.758  1   e+ 0
##  2 kappa                <NA>     0.114   NA        NA     NA       
##  3 mcnemar              <NA>    NA       NA        NA      5.21e-18
##  4 sensitivity          good     0.825   NA        NA     NA       
##  5 specificity          good     0.304   NA        NA     NA       
##  6 pos_pred_value       good     0.873   NA        NA     NA       
##  7 neg_pred_value       good     0.230   NA        NA     NA       
##  8 precision            good     0.873   NA        NA     NA       
##  9 recall               good     0.825   NA        NA     NA       
## 10 f1                   good     0.848   NA        NA     NA       
## 11 prevalence           good     0.853   NA        NA     NA       
## 12 detection_rate       good     0.704   NA        NA     NA       
## 13 detection_prevalence good     0.806   NA        NA     NA       
## 14 balanced_accuracy    good     0.565   NA        NA     NA

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.

library(pROC)
(roc_score=roc(test$good, test$pred)) #AUC score

## 
## Call:
## roc.default(response = test$good, predictor = test$pred)
## 
## Data: test$pred in 1249 controls (test$good bad) < 7258 cases (test$good good).
## Area under the curve: 0.6384

plot(roc_score ,main ="ROC curve -- Logistic Regression ")

References

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

Program R untuk Analisis Regresi Logistik

Rahma Anisa

2023-05-05

Data Exploration

The Model, Estimation and Inference

Interpreting coefficients

Prediction

Categorical Independent Variable

Multiple Logistic Regression

Model Evaluation

Pseudo \(R^2\)

Residual Assessment

Testing logistic model out of sample

Model Evaluation

VIF

Confusion Matrix

ROC Curve

References