Logistics Regression
G1501202053 Ikhlasul Amalia Rahmi
9/29/2020
Data Exploration Sebagai ilustrasi, akan digunakan data Default yang berisi data simulasi sebanyak 10,000 konsumen kartu kredit.
library(ISLR)## Warning: package 'ISLR' was built under R version 4.0.5library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.0.5## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.6 v dplyr 1.0.8
## v tidyr 1.2.0 v stringr 1.4.0
## v readr 2.1.2 v forcats 0.5.1## Warning: package 'ggplot2' was built under R version 4.0.5## Warning: package 'tibble' was built under R version 4.0.5## Warning: package 'tidyr' was built under R version 4.0.5## Warning: package 'readr' was built under R version 4.0.5## Warning: package 'dplyr' was built under R version 4.0.5## Warning: package 'forcats' was built under R version 4.0.4## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(DataExplorer)
data(Default)
help(Default)## starting httpd help server ...## doneglimpse(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~view(Default)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. :73554plot_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.04313859Default %>% 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: 8Default %>%
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.005498917Seperti 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, misalnya odds kejadian kredit macet terhadap kredit lancar.
exp(reglog$coefficients)## (Intercept) balance
## 2.366933e-05 1.005514e+00Pada 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.005943365Prediction 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.005752145Artinya, 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: 6predict(reglog2, newdata = data.frame(student="Yes"), type="res")## 1
## 0.04313859predict(reglog2, newdata = data.frame(student="No"), type="res")## 1
## 0.02919501Multiple 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: 8Perhatikan bahwa pola pengaruh peubah penjelas tidak sama seperti yang dihasilkan pada model sebelumnya yang hanya menggunakan satu peubah penjelas.
Model Evaluation Pseudo R2
library(pscl)## Classes and Methods for R developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University
## Simon Jackman
## hurdle and zeroinfl functions by Achim Zeileispscl::pR2(reglog)["McFadden"]## fitting null model for pseudo-r2## McFadden
## 0.4533916pscl::pR2(full.mod)["McFadden"]## fitting null model for pseudo-r2## McFadden
## 0.4619194Residual 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
##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.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 rowssummary(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 :18758Pada ilustrasi ini, peubah fico, dti, loan_amnt, dan purpose akan digunakan sebagai peubah penjelas.
table(loan$good)##
## bad good
## 6371 36164loan$good<-as.factor(loan$good)
loan %>%
rename()## # 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 rowsset.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: 27660Model 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.010285Confusion 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)## Warning: package 'caret' was built under R version 4.0.5## Loading required package: lattice## Warning: package 'lattice' was built under R version 4.0.5##
## Attaching package: 'caret'## The following object is masked from 'package:purrr':
##
## lifttest$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 x 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 NAROC 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(InformationValue)##
## Attaching package: 'InformationValue'## The following objects are masked from 'package:caret':
##
## confusionMatrix, precision, sensitivity, specificityplotROC(test$good=="good", test$pred)