Model logitowy
Marcin Bukowski
02 czerwiec 2023
Wprowadzenie
Model regresji logistycznej wykorzystywany jest do objaśniania zmiennych jakościowych, których wartości są z zakresu od 0 do 1. W literaturze przedmiotu model ten powszechnie znany jest jako logitowy. Model regresji logistycznej oparty jest na funkcji logistycznej. Funkcję tę określa się wzorem: \[ f(z) = \frac{e^{z}}{1+e^{z}} = \frac{1}{1+e^{-z}} \] gdzie: \[ z = \beta_{0}+\sum_{i = 1}^{n}\beta_{i}x_{i} \]
Zastosowanie modelu regresji logistycznej pozwala określić zarówno siłę, jak i kierunek zależności między czynnikiem jakościowym (typu klasowego) lub ilościowym (typu dyskretnego lub ciągłego) a dychotomiczną zmienną objaśnianą.
Przykładowy model z jedną zmieną ilościową i jakościową
Dla lepszego zrozumienia czym jest model logitowy i na czym polega
obliczenie wartości na jego podstawie można zbudować przykładowy model
okreslający jakość kredytu w zalezności od stanu konta czekowego
(zmienna jakościowa) i kwoty kredytu (zmienna ilościowa). Do budowy
modelu zostanie wykorzystana funkcja glm z arumentem
family=binomial().
library(ROCR)
load("C:/Users/mbuko/OneDrive/Dokumenty/credit.RData")
credit2<-credit[,c(1, 5, 21)]
model<-glm(jakosc ~ .,
data=credit2,family=binomial())Podsumowanie modelu może zostać wywołane za pomocą funkcji
summary.
summary(model)##
## Call:
## glm(formula = jakosc ~ ., family = binomial(), data = credit2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8041 -0.8943 -0.4744 1.1095 2.2149
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.853e-01 1.456e-01 -2.646 0.008135 **
## konto_czekowe>200 -1.120e+00 3.280e-01 -3.415 0.000637 ***
## konto_czekowe0-200 -5.024e-01 1.777e-01 -2.828 0.004687 **
## konto_czekowebrak -2.027e+00 2.003e-01 -10.120 < 2e-16 ***
## kwota 1.131e-04 2.566e-05 4.410 1.04e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1221.7 on 999 degrees of freedom
## Residual deviance: 1070.5 on 995 degrees of freedom
## AIC: 1080.5
##
## Number of Fisher Scoring iterations: 4Jednym z elementów podsumowania są oszacowane parametry strukturalne modelu, które mogą być wywołane za pomocą:
wsp<-model$coefficients
wsp## (Intercept) konto_czekowe>200 konto_czekowe0-200 konto_czekowebrak
## -0.3852792592 -1.1201088317 -0.5023754383 -2.0267377026
## kwota
## 0.0001131359Wśród oszacowanych wartości parametrów strukturalnych modelu brak jest parametru dla zmiennej stan konta czekowego <0. Wynika to z faktu, iż ten poziom zmiennej został przyjęty jako tzw. poziom odniesenia, w stosunku do którego obliczany jest iloraz szans wystąpienia “złego” kredytu przy pozostałych poziomach zmiennej stan konta czekowego
Aby zinterpretować oszacowane parametry strukturalne modelu należy
obliczyć wartość wyrażenia \(e^{\beta_{i}}\) dla \(i \in \{0, \dots, 4\}\) Można to wykonać za
pomocą funkcji exp.
exp(wsp)## (Intercept) konto_czekowe>200 konto_czekowe0-200 konto_czekowebrak
## 0.6802606 0.3262443 0.6050916 0.1317647
## kwota
## 1.0001131Interpretacja otrzymanych wartości jest następująca:
- dla parametru \(\beta_{0}\) -
parametr ten to tzw. wyraz wolny - obliczona wartość 0.6802606 nie może
być w tym przypadku interpretowana - teoretycznie oznaczałaby szansę
złego kredytu w sytuacji kiedy stan konta czekowego wnioskodawcy jest
<0 (poziom odniesienia), a wnioskowana kwota kredytu wynosi
0;
- dla parametru \(\beta_{1}\) przy
zmiennej konto_czekowe>200 - obliczona wartość wyniosła
0.6802606. Wartość ta jest mniejsza od 1, co oznacza, że szansa
złego kredytu przy koncie czekowym >200 jest mniejsza od
szansy przy stanie konta czekowego <0 (poziom odniesienia).
Obliczając wartość wyrażenia 1- 0.6802606 otrzymujemy 0.3197394, co
oznacza, że w przypadku konta czekowego >200 szansa złego
kredytu jest o tyle mniejsze w porównaniu z sytuacją konta czekowego
<0 przy tej samej wnioskowanej kwocie kredytu;
- w podobny sposób można zinterpretować parametry dla kolejnych
poziomów stanu konta czekowego;
- dla parametru \(\beta_{4}\) przy
zmiennej kwota - obliczona wartość 1.0001131 jest dodatnia, co
oznacza, że wraz ze wzrostem kwoty kredytu szansa złego kredytu
rośnie. Wartość ta jest jednak niewiele większa od 1 co oznacza, że przy
wzroście kwoty kredytu o 1 Euro szansa złego kredytu rośnie
bardzo niewiele. Aby lepiej zrozumieć zależność między jakością kredytu
a kwotą kredytu można założyć wzrost kwoty kredytu o 1000 Euro. Wówczas
szansa złego kredytu przy takim wzroście kwoty rośnie o
0.119784 razy - przy tym samym poziomie zmiennej
konto_czekowe.
Oszacowane wartości parametrów modelu w zależności od poziomu zmiennej konto_czekowe zostały dodane do tabeli z danymicredit2
credit2$wsp_konto<-ifelse(credit2$konto_czekowe == levels(credit2$konto_czekowe)[1], 0,
ifelse(credit2$konto_czekowe == levels(credit2$konto_czekowe)[2], wsp[2],
ifelse(credit2$konto_czekowe == levels(credit2$konto_czekowe)[3], wsp[3],wsp[4])))
head(credit2,10)## konto_czekowe kwota jakosc wsp_konto
## 1 < 0 1169 dobry 0.0000000
## 2 0-200 5951 zly -0.5023754
## 3 brak 2096 dobry -2.0267377
## 4 < 0 7882 dobry 0.0000000
## 5 < 0 4870 zly 0.0000000
## 6 brak 9055 dobry -2.0267377
## 7 brak 2835 dobry -2.0267377
## 8 0-200 6948 dobry -0.5023754
## 9 brak 3059 dobry -2.0267377
## 10 0-200 5234 zly -0.5023754Na tej podstawie można obliczyć wartość funkcji - z zgodnie z zależnością przedstawioną we wprowadzeniu:
credit2$'-z'<-(-(wsp[1] + credit2$wsp_konto + credit2$kwota*wsp[5]))
head(credit2,10)## konto_czekowe kwota jakosc wsp_konto -z
## 1 < 0 1169 dobry 0.0000000 0.2530234
## 2 0-200 5951 zly -0.5023754 0.2143832
## 3 brak 2096 dobry -2.0267377 2.1748842
## 4 < 0 7882 dobry 0.0000000 -0.5064575
## 5 < 0 4870 zly 0.0000000 -0.1656923
## 6 brak 9055 dobry -2.0267377 1.3875718
## 7 brak 2835 dobry -2.0267377 2.0912768
## 8 0-200 6948 dobry -0.5023754 0.1015868
## 9 brak 3059 dobry -2.0267377 2.0659344
## 10 0-200 5234 zly -0.5023754 0.2955017Podstawiając obliczone wartości - z do wzoru funkcji logitowej otrzymujemy następujące wartości prawdopodobieństwa złego kredytu:
credit2$praw_zly<-1/(1+exp(credit2$`-z`))
head(credit2,10)## konto_czekowe kwota jakosc wsp_konto -z praw_zly
## 1 < 0 1169 dobry 0.0000000 0.2530234 0.4370795
## 2 0-200 5951 zly -0.5023754 0.2143832 0.4466085
## 3 brak 2096 dobry -2.0267377 2.1748842 0.1020287
## 4 < 0 7882 dobry 0.0000000 -0.5064575 0.6239757
## 5 < 0 4870 zly 0.0000000 -0.1656923 0.5413286
## 6 brak 9055 dobry -2.0267377 1.3875718 0.1997957
## 7 brak 2835 dobry -2.0267377 2.0912768 0.1099476
## 8 0-200 6948 dobry -0.5023754 0.1015868 0.4746251
## 9 brak 3059 dobry -2.0267377 2.0659344 0.1124522
## 10 0-200 5234 zly -0.5023754 0.2955017 0.4266575Te same wartości można otrzymać bezpośrednio z wyników obliczeń
modelu poprzez wywołanie model$fitted.values. Poniżej 10
pierwszych wartości:
head(model$fitted.values, 10)## 1 2 3 4 5 6 7 8
## 0.4370795 0.4466085 0.1020287 0.6239757 0.5413286 0.1997957 0.1099476 0.4746251
## 9 10
## 0.1124522 0.4266575Obliczona wartość prawdopodobieństwa może zostać wykorzystana do określenia jakości kredytów. Wymaga to jednak przyjęcia granicznej wartości prawdopodobieństwa oddzielające kredyty dobre od złych. Na wstępnym etapie obliczeń można przyjąć p=0,5 co oznacza, że jeśli obliczona wartość złego kredytu jest mniejsza równa 50% to kredyt jest dobry, w przeciwnym razie zły.
p<-0.5
credit2$oszacowana_jakosc<-ifelse(credit2$praw_zly > p, "zly", "dobry")
head(model$fitted.values, 10)## 1 2 3 4 5 6 7 8
## 0.4370795 0.4466085 0.1020287 0.6239757 0.5413286 0.1997957 0.1099476 0.4746251
## 9 10
## 0.1124522 0.4266575Ocenę jakości modelu na podstawie macierzy błędów oraz procentu błędów przedstawiono poniżej.
table(credit2$oszacowana_jakosc, credit2$jakosc, dnn=c("Predykcja","Rzeczywiste"))## Rzeczywiste
## Predykcja dobry zly
## dobry 649 230
## zly 51 70sum(credit2$jakosc != credit2$oszacowana_jakosc) / nrow(credit2)## [1] 0.281Model dla wszystkich zmiennych
Poniżej przedstawiono procedurę budowy modelu logitowego dla
wszystkich zmiennych ze zbioru credit.
Podział danych na zbiór uczący i testowy
set.seed(123) # for reproducibility
random <- sample(1:nrow(credit), 0.75 * nrow(credit))
data_train<-credit[random,]
data_test <- credit[-random, ]Budowa modelu ze wszystkimi zmiennymi
Za pomocą funkcji glm z wykorzystaniem argumentu
family=binomial() zbudowano model logitowy dla wszystkich
zmiennych ze zbioru treningowego.
model<-glm(jakosc~.,data=data_train,family=binomial())
summary(model)##
## Call:
## glm(formula = jakosc ~ ., family = binomial(), data = data_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0381 -0.6965 -0.3608 0.6571 2.6938
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.157e-02 1.643e+00 -0.044 0.965263
## konto_czekowe>200 -8.553e-01 4.266e-01 -2.005 0.044971 *
## konto_czekowe0-200 -2.704e-01 2.546e-01 -1.062 0.288159
## konto_czekowebrak -1.628e+00 2.752e-01 -5.917 3.28e-09 ***
## czas 3.087e-02 1.050e-02 2.939 0.003287 **
## historiaistniejace_spł -5.309e-01 5.386e-01 -0.986 0.324296
## historiakrytyczne -1.235e+00 5.414e-01 -2.281 0.022561 *
## historiaopoznienia -6.739e-01 5.937e-01 -1.135 0.256344
## historiawszystkie_spł 3.703e-01 6.671e-01 0.555 0.578838
## celbiznes -1.128e+00 9.325e-01 -1.210 0.226397
## celedukacja -7.046e-01 9.561e-01 -0.737 0.461164
## celinne -1.496e+00 1.244e+00 -1.202 0.229237
## celmeble -1.028e+00 8.838e-01 -1.163 0.244692
## celprzekwalifikowanie -1.478e+00 1.550e+00 -0.953 0.340372
## celremont 3.154e-01 1.085e+00 0.291 0.771341
## celRTV -1.120e+00 8.777e-01 -1.276 0.201992
## celsam_nowy -3.092e-01 8.806e-01 -0.351 0.725471
## celsam_uzyw -2.054e+00 9.522e-01 -2.157 0.031030 *
## kwota 1.549e-04 5.006e-05 3.095 0.001970 **
## oszczednosci>1000 -1.461e+00 5.889e-01 -2.482 0.013080 *
## oszczednosci100-500 -4.294e-01 3.424e-01 -1.254 0.209721
## oszczednosci500-1000 -8.898e-01 5.222e-01 -1.704 0.088414 .
## oszczednoscibrak -1.162e+00 3.192e-01 -3.642 0.000271 ***
## staz_pracy>7 -6.856e-01 3.437e-01 -1.995 0.046047 *
## staz_pracy1-4 -3.356e-01 2.802e-01 -1.198 0.230929
## staz_pracy4-7 -1.223e+00 3.585e-01 -3.411 0.000648 ***
## staz_pracybezrobotny -5.705e-01 5.122e-01 -1.114 0.265339
## `rata_%doch` 3.493e-01 1.016e-01 3.438 0.000585 ***
## plecM -4.650e-01 2.218e-01 -2.097 0.035986 *
## poreczycieltak -2.401e-01 3.470e-01 -0.692 0.488971
## zamieszkanie -1.291e-03 1.011e-01 -0.013 0.989804
## zabezpieczenienieruchomosc -1.060e+00 5.690e-01 -1.864 0.062370 .
## zabezpieczeniesamochod -5.353e-01 5.391e-01 -0.993 0.320775
## zabezpieczenieubezpieczenie -7.665e-01 5.479e-01 -1.399 0.161803
## wiek -4.076e-03 1.036e-02 -0.393 0.694007
## inne_zobowbrak -3.943e-01 2.779e-01 -1.419 0.155845
## inne_zobowsklep -3.443e-01 4.845e-01 -0.711 0.477331
## rodzaj_mieszwlasne -5.482e-01 2.795e-01 -1.961 0.049845 *
## rodzaj_mieszza_darmo -9.899e-01 6.111e-01 -1.620 0.105250
## l_kredytow 2.102e-01 2.292e-01 0.917 0.359080
## kwalifikacjeniewykwal_rez -8.665e-02 7.846e-01 -0.110 0.912067
## kwalifikacjewykwalifikowany 3.494e-02 7.571e-01 0.046 0.963194
## kwalifikacjewysokie_kwal -2.056e-01 7.542e-01 -0.273 0.785131
## l_osob 3.521e-01 2.853e-01 1.234 0.217188
## telefontak -4.038e-01 2.351e-01 -1.718 0.085846 .
## obcokrajowiectak 1.429e+00 7.694e-01 1.857 0.063352 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 905.90 on 749 degrees of freedom
## Residual deviance: 662.55 on 704 degrees of freedom
## AIC: 754.55
##
## Number of Fisher Scoring iterations: 5Dobór zmiennych do modelu
Dobór zmiennych do modelu został przeprowadzony metodą krokową wstecz za pomocą funkcji `step’. Wynikiem tej funkcji są kolejne modele począwszy od modelu zbudowanego dla wszystkich zmiennych. W kolejnych etapach eliminowana jest jedna zmienna do momentu uzyskania modelu najbardziej dopasowanego. Ocena dopasowania modelu przeprowadzana jest na podstawie kryterium AIC (Akaike Information Criterion) - im mniejsza wartość AIC tym dopasowanie modelu jest lepsze. Eliminacja zmiennych przeprowadzana jest do momentu, kiedy kolejna eliminacja nie powoduje dalszego spadku AIC.
model_opt<-step(model)## Start: AIC=754.55
## jakosc ~ konto_czekowe + czas + historia + cel + kwota + oszczednosci +
## staz_pracy + `rata_%doch` + plec + poreczyciel + zamieszkanie +
## zabezpieczenie + wiek + inne_zobow + rodzaj_miesz + l_kredytow +
## kwalifikacje + l_osob + telefon + obcokrajowiec
##
## Df Deviance AIC
## - kwalifikacje 3 663.24 749.24
## - inne_zobow 2 664.55 752.55
## - zamieszkanie 1 662.55 752.55
## - wiek 1 662.71 752.71
## - poreczyciel 1 663.04 753.04
## - l_kredytow 1 663.39 753.39
## - zabezpieczenie 3 667.93 753.93
## - l_osob 1 664.06 754.06
## <none> 662.55 754.55
## - rodzaj_miesz 2 667.42 755.42
## - telefon 1 665.54 755.54
## - obcokrajowiec 1 666.74 756.74
## - plec 1 666.94 756.94
## - staz_pracy 4 675.87 759.87
## - historia 4 676.42 760.42
## - czas 1 671.31 761.31
## - kwota 1 672.33 762.33
## - cel 9 688.51 762.51
## - `rata_%doch` 1 674.86 764.86
## - oszczednosci 4 683.82 767.82
## - konto_czekowe 3 707.11 793.11
##
## Step: AIC=749.24
## jakosc ~ konto_czekowe + czas + historia + cel + kwota + oszczednosci +
## staz_pracy + `rata_%doch` + plec + poreczyciel + zamieszkanie +
## zabezpieczenie + wiek + inne_zobow + rodzaj_miesz + l_kredytow +
## l_osob + telefon + obcokrajowiec
##
## Df Deviance AIC
## - zamieszkanie 1 663.24 747.24
## - inne_zobow 2 665.28 747.28
## - wiek 1 663.46 747.46
## - poreczyciel 1 663.64 747.64
## - l_kredytow 1 664.16 748.16
## - l_osob 1 664.67 748.67
## - zabezpieczenie 3 668.86 748.86
## <none> 663.24 749.24
## - rodzaj_miesz 2 667.96 749.96
## - telefon 1 667.12 751.12
## - obcokrajowiec 1 667.56 751.56
## - plec 1 667.67 751.67
## - staz_pracy 4 676.76 754.76
## - historia 4 677.21 755.21
## - kwota 1 672.50 756.50
## - czas 1 672.77 756.77
## - cel 9 689.74 757.74
## - `rata_%doch` 1 675.27 759.27
## - oszczednosci 4 684.11 762.11
## - konto_czekowe 3 707.94 787.94
##
## Step: AIC=747.24
## jakosc ~ konto_czekowe + czas + historia + cel + kwota + oszczednosci +
## staz_pracy + `rata_%doch` + plec + poreczyciel + zabezpieczenie +
## wiek + inne_zobow + rodzaj_miesz + l_kredytow + l_osob +
## telefon + obcokrajowiec
##
## Df Deviance AIC
## - inne_zobow 2 665.28 745.28
## - wiek 1 663.46 745.46
## - poreczyciel 1 663.65 745.65
## - l_kredytow 1 664.16 746.16
## - l_osob 1 664.68 746.68
## - zabezpieczenie 3 668.86 746.86
## <none> 663.24 747.24
## - rodzaj_miesz 2 668.10 748.10
## - telefon 1 667.14 749.14
## - obcokrajowiec 1 667.58 749.58
## - plec 1 667.70 749.70
## - staz_pracy 4 676.88 752.88
## - historia 4 677.22 753.22
## - kwota 1 672.50 754.50
## - czas 1 672.78 754.78
## - cel 9 689.74 755.74
## - `rata_%doch` 1 675.31 757.31
## - oszczednosci 4 684.16 760.16
## - konto_czekowe 3 708.17 786.17
##
## Step: AIC=745.28
## jakosc ~ konto_czekowe + czas + historia + cel + kwota + oszczednosci +
## staz_pracy + `rata_%doch` + plec + poreczyciel + zabezpieczenie +
## wiek + rodzaj_miesz + l_kredytow + l_osob + telefon + obcokrajowiec
##
## Df Deviance AIC
## - wiek 1 665.41 743.41
## - poreczyciel 1 665.61 743.61
## - l_kredytow 1 666.14 744.14
## - l_osob 1 666.85 744.85
## <none> 665.28 745.28
## - zabezpieczenie 3 671.71 745.71
## - rodzaj_miesz 2 670.28 746.28
## - obcokrajowiec 1 669.16 747.16
## - telefon 1 669.31 747.31
## - plec 1 669.44 747.44
## - staz_pracy 4 678.59 750.59
## - kwota 1 674.11 752.11
## - czas 1 675.11 753.11
## - cel 9 691.18 753.18
## - historia 4 681.79 753.79
## - `rata_%doch` 1 676.99 754.99
## - oszczednosci 4 685.99 757.99
## - konto_czekowe 3 709.59 783.59
##
## Step: AIC=743.41
## jakosc ~ konto_czekowe + czas + historia + cel + kwota + oszczednosci +
## staz_pracy + `rata_%doch` + plec + poreczyciel + zabezpieczenie +
## rodzaj_miesz + l_kredytow + l_osob + telefon + obcokrajowiec
##
## Df Deviance AIC
## - poreczyciel 1 665.75 741.75
## - l_kredytow 1 666.23 742.23
## - l_osob 1 666.94 742.94
## <none> 665.41 743.41
## - zabezpieczenie 3 671.81 743.81
## - rodzaj_miesz 2 670.80 744.80
## - obcokrajowiec 1 669.31 745.31
## - telefon 1 669.63 745.63
## - plec 1 669.67 745.67
## - staz_pracy 4 679.31 749.31
## - kwota 1 674.24 750.24
## - cel 9 691.18 751.18
## - czas 1 675.47 751.47
## - historia 4 682.06 752.06
## - `rata_%doch` 1 677.19 753.19
## - oszczednosci 4 686.21 756.21
## - konto_czekowe 3 709.69 781.69
##
## Step: AIC=741.75
## jakosc ~ konto_czekowe + czas + historia + cel + kwota + oszczednosci +
## staz_pracy + `rata_%doch` + plec + zabezpieczenie + rodzaj_miesz +
## l_kredytow + l_osob + telefon + obcokrajowiec
##
## Df Deviance AIC
## - l_kredytow 1 666.53 740.53
## - l_osob 1 667.27 741.27
## <none> 665.75 741.75
## - zabezpieczenie 3 672.44 742.44
## - rodzaj_miesz 2 671.01 743.01
## - telefon 1 669.87 743.87
## - obcokrajowiec 1 669.96 743.96
## - plec 1 670.07 744.07
## - staz_pracy 4 679.93 747.93
## - kwota 1 674.59 748.59
## - czas 1 675.78 749.78
## - cel 9 692.26 750.26
## - historia 4 682.26 750.26
## - `rata_%doch` 1 677.85 751.85
## - oszczednosci 4 686.22 754.22
## - konto_czekowe 3 709.75 779.75
##
## Step: AIC=740.53
## jakosc ~ konto_czekowe + czas + historia + cel + kwota + oszczednosci +
## staz_pracy + `rata_%doch` + plec + zabezpieczenie + rodzaj_miesz +
## l_osob + telefon + obcokrajowiec
##
## Df Deviance AIC
## - l_osob 1 668.29 740.29
## <none> 666.53 740.53
## - zabezpieczenie 3 672.93 740.93
## - rodzaj_miesz 2 671.60 741.60
## - telefon 1 670.35 742.35
## - obcokrajowiec 1 670.73 742.73
## - plec 1 670.76 742.76
## - staz_pracy 4 680.22 746.22
## - kwota 1 675.32 747.32
## - czas 1 676.30 748.30
## - historia 4 682.62 748.62
## - cel 9 693.19 749.19
## - `rata_%doch` 1 678.40 750.40
## - oszczednosci 4 687.10 753.10
## - konto_czekowe 3 710.55 778.55
##
## Step: AIC=740.29
## jakosc ~ konto_czekowe + czas + historia + cel + kwota + oszczednosci +
## staz_pracy + `rata_%doch` + plec + zabezpieczenie + rodzaj_miesz +
## telefon + obcokrajowiec
##
## Df Deviance AIC
## <none> 668.29 740.29
## - zabezpieczenie 3 674.77 740.77
## - rodzaj_miesz 2 673.05 741.05
## - plec 1 671.57 741.57
## - telefon 1 672.08 742.08
## - obcokrajowiec 1 672.37 742.37
## - staz_pracy 4 681.62 745.62
## - kwota 1 676.53 746.53
## - czas 1 678.06 748.06
## - historia 4 684.75 748.75
## - cel 9 695.25 749.25
## - `rata_%doch` 1 679.38 749.38
## - oszczednosci 4 688.23 752.23
## - konto_czekowe 3 712.35 778.35Podsumowanie modelu optymalnego:
summary(model_opt)##
## Call:
## glm(formula = jakosc ~ konto_czekowe + czas + historia + cel +
## kwota + oszczednosci + staz_pracy + `rata_%doch` + plec +
## zabezpieczenie + rodzaj_miesz + telefon + obcokrajowiec,
## family = binomial(), data = data_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2323 -0.7021 -0.3759 0.6649 2.7049
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.212e-01 1.338e+00 0.240 0.810315
## konto_czekowe>200 -9.290e-01 4.190e-01 -2.217 0.026611 *
## konto_czekowe0-200 -3.177e-01 2.505e-01 -1.268 0.204667
## konto_czekowebrak -1.617e+00 2.726e-01 -5.932 3e-09 ***
## czas 3.210e-02 1.033e-02 3.106 0.001894 **
## historiaistniejace_spł -7.540e-01 5.158e-01 -1.462 0.143801
## historiakrytyczne -1.299e+00 5.381e-01 -2.415 0.015746 *
## historiaopoznienia -7.582e-01 5.897e-01 -1.286 0.198572
## historiawszystkie_spł 2.784e-01 6.299e-01 0.442 0.658501
## celbiznes -1.002e+00 9.250e-01 -1.083 0.278595
## celedukacja -6.173e-01 9.529e-01 -0.648 0.517128
## celinne -1.497e+00 1.220e+00 -1.227 0.219927
## celmeble -1.025e+00 8.801e-01 -1.165 0.244151
## celprzekwalifikowanie -1.518e+00 1.537e+00 -0.988 0.323267
## celremont 3.454e-01 1.079e+00 0.320 0.748831
## celRTV -1.097e+00 8.721e-01 -1.258 0.208265
## celsam_nowy -2.562e-01 8.752e-01 -0.293 0.769677
## celsam_uzyw -1.964e+00 9.434e-01 -2.082 0.037356 *
## kwota 1.374e-04 4.828e-05 2.847 0.004418 **
## oszczednosci>1000 -1.288e+00 5.639e-01 -2.284 0.022365 *
## oszczednosci100-500 -3.919e-01 3.371e-01 -1.163 0.244979
## oszczednosci500-1000 -8.982e-01 5.073e-01 -1.770 0.076659 .
## oszczednoscibrak -1.110e+00 3.139e-01 -3.536 0.000406 ***
## staz_pracy>7 -6.941e-01 3.184e-01 -2.180 0.029272 *
## staz_pracy1-4 -3.798e-01 2.757e-01 -1.377 0.168453
## staz_pracy4-7 -1.194e+00 3.510e-01 -3.401 0.000671 ***
## staz_pracybezrobotny -6.533e-01 4.471e-01 -1.461 0.144031
## `rata_%doch` 3.225e-01 9.861e-02 3.270 0.001074 **
## plecM -3.871e-01 2.134e-01 -1.814 0.069754 .
## zabezpieczenienieruchomosc -1.068e+00 5.545e-01 -1.926 0.054061 .
## zabezpieczeniesamochod -4.929e-01 5.271e-01 -0.935 0.349705
## zabezpieczenieubezpieczenie -7.716e-01 5.394e-01 -1.431 0.152562
## rodzaj_mieszwlasne -5.234e-01 2.644e-01 -1.980 0.047735 *
## rodzaj_mieszza_darmo -9.284e-01 5.971e-01 -1.555 0.120008
## telefontak -4.198e-01 2.172e-01 -1.932 0.053319 .
## obcokrajowiectak 1.375e+00 7.484e-01 1.837 0.066150 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 905.90 on 749 degrees of freedom
## Residual deviance: 668.29 on 714 degrees of freedom
## AIC: 740.29
##
## Number of Fisher Scoring iterations: 5Przedziały ufności dla oszacowanych wartości parametrów modelu:
confint.default(model_opt)## 2.5 % 97.5 %
## (Intercept) -2.301949e+00 2.9444265965
## konto_czekowe>200 -1.750286e+00 -0.1077777037
## konto_czekowe0-200 -8.085939e-01 0.1732233881
## konto_czekowebrak -2.151647e+00 -1.0828966295
## czas 1.184567e-02 0.0523491876
## historiaistniejace_spł -1.764941e+00 0.2569640202
## historiakrytyczne -2.354116e+00 -0.2447296475
## historiaopoznienia -1.914059e+00 0.3976821042
## historiawszystkie_spł -9.561626e-01 1.5129665159
## celbiznes -2.815310e+00 0.8107817156
## celedukacja -2.484885e+00 1.2503647073
## celinne -3.889248e+00 0.8948976578
## celmeble -2.749979e+00 0.6999297979
## celprzekwalifikowanie -4.531057e+00 1.4943982269
## celremont -1.768965e+00 2.4597750555
## celRTV -2.806739e+00 0.6118858268
## celsam_nowy -1.971549e+00 1.4590541154
## celsam_uzyw -3.813065e+00 -0.1149854402
## kwota 4.280515e-05 0.0002320459
## oszczednosci>1000 -2.393229e+00 -0.1827833649
## oszczednosci100-500 -1.052650e+00 0.2687863135
## oszczednosci500-1000 -1.892506e+00 0.0961650320
## oszczednoscibrak -1.725192e+00 -0.4948235870
## staz_pracy>7 -1.318104e+00 -0.0699972700
## staz_pracy1-4 -9.202074e-01 0.1606974273
## staz_pracy4-7 -1.881941e+00 -0.5059005911
## staz_pracybezrobotny -1.529638e+00 0.2231348076
## `rata_%doch` 1.292271e-01 0.5157813157
## plecM -8.054366e-01 0.0312618831
## zabezpieczenienieruchomosc -2.154882e+00 0.0186396274
## zabezpieczeniesamochod -1.526035e+00 0.5401781727
## zabezpieczenieubezpieczenie -1.828686e+00 0.2855481210
## rodzaj_mieszwlasne -1.041600e+00 -0.0052234275
## rodzaj_mieszza_darmo -2.098745e+00 0.2419717715
## telefontak -8.455919e-01 0.0060036281
## obcokrajowiectak -9.172384e-02 2.8419184093Na podstawie podsumowania modelu można również określić poziom istotności poszczególnych zmiennych na podstawie wartości prawdopodobieństwa Pr(>|z|). Obliczoną wartość prawdopodobieństwa należy porównać z przyjętym poziomem istotności \(\alpha = 0,05\). Jeżeli prawdopodobieństwo Pr(>|z|) jest mniejsze od przyjętego poziomu \(\alpha\) wówczas zmienna jest istotna statystycznie i wpływa na jakość kredytu. W przedstawionym podsumowaniu zmienne te oznaczone są symbolem jednej, dwóch lub trzech *. Aby wskazać tylko te zmienne można wykorzystać poniższy kod:
significant.variables <-summary(model_opt)$coeff[-1,4] < 0.05
names(significant.variables)[significant.variables == TRUE]## [1] "konto_czekowe>200" "konto_czekowebrak" "czas"
## [4] "historiakrytyczne" "celsam_uzyw" "kwota"
## [7] "oszczednosci>1000" "oszczednoscibrak" "staz_pracy>7"
## [10] "staz_pracy4-7" "`rata_%doch`" "rodzaj_mieszwlasne"Obliczenie prawdopodobieństwa złego kredytu na podstawie optymalnego modelu
Poniżej przestawiono sposób obliczenia wartości prawdopodobieństwa złego kredytu na podstawie najlepiej dopasowanego modelu. Obliczona wartość prawdopodobieństwa została dołączono do danych ze zbioru uczącego w postaci zmiennej model.
data_train$model<-predict(model_opt,type="response")
head(data_train$model, 10)## [1] 0.53780186 0.57415193 0.08938812 0.45797020 0.67870592 0.49838413
## [7] 0.02422657 0.05105906 0.06996504 0.06993486Ocena modelu dla p = 50%
Podobnie jak w przypadku zadaniu przykładowym ocenę jakości modelu można przeprowadzić przyjmując p=0,5 co oznacza, że jeśli obliczona wartość złego kredytu jest mniejsza równa 50% to kredyt jest dobry, w przeciwnym razie zły.
p<-0.5
data_train$oszacowana_jakosc<-ifelse(data_train$model > p, "zly", "dobry")
table(data_train$oszacowana_jakosc, data_train$jakosc, dnn=c("Predykcja","Rzeczywiste"))## Rzeczywiste
## Predykcja dobry zly
## dobry 479 102
## zly 52 117sum(data_train$jakosc != data_train$oszacowana_jakosc) / nrow(data_train)## [1] 0.2053333Ocena modelu dla zbioru testowego
Obliczenie prawdopodobieństwa złego kredytu
Ocenę jakości modelu przeprowadzono na podstawie danych znajdujących się w zbiorze danych testowych. W tym celu obliczono wartość prawdopodobieństwa dla wszystkich 250 obserwacji w następujący sposób:
data_test$model <- predict(model_opt,type='response',data_test)Ocena jakości modelu
Krzywa ROC
W poszukiwaniu optymalnego modelu predykcyjnego pomocne mogą być wykorzystane techniki graficzne – do takich należy krzywa ROC (z ang. Receiver Operating Characteristic). Krzywa ROC jest wizualizacją zależności pomiędzy skutecznością klasyfikacji pozytywnych (Czułość) a nieskutecznością klasyfikacji przypadków negatywnych (1-specyficzność) na każdym poziomie prawdopodobieństwa. Krzywa ROC obrazuje więc, jak duży będzie odsetek błędnych klasyfikacji (pozytywnych i negatywnych) dla danego punktu odcięcia. W jaki sposób założone graniczne prawdopodobieństwo przypisania do kategorii pozytywnej (wystąpienie zjawiska) wpływa równocześnie na poprawną i błędną klasyfikację tychże przypadków.
m1_pred <- prediction(data_test$model, data_test$jakosc)
m1_perf <- performance(m1_pred,"tpr","fpr")
#ROC
plot(m1_perf, lwd=2, colorize=TRUE, main="Krzywa ROC", xlab = "1-specyficzność", ylab = "Czułość")
lines(x=c(0, 1), y=c(0, 1), col="red", lwd=1, lty=3);
lines(x=c(1, 0), y=c(0, 1), col="green", lwd=1, lty=4)
Krzywa precyzja - czułość
Krzywa precyzja-czułość przedstawia zależność między precyzją a czułością dla różnych progów. Wysoki obszar pod krzywą oznacza zarówno wysoką czułość, jak i wysoką precyzję, przy czym wysoka precyzja odnosi się do niskiego wskaźnika wyników fałszywie dodatnich, a wysoka czułość odnosi się do niskiego wskaźnika wyników wyników fałszywie ujemnych. Wysokie wyniki dla obu pokazują, że klasyfikator zwraca dokładne wyniki (wysoka precyzja), jak również zwraca większość wszystkich pozytywnych wyników (wysoka czułość).
m1_perf_precision <- performance(m1_pred, measure = "prec", x.measure = "rec")
plot(m1_perf_precision, colorize=TRUE, main="Krzywa precyzja - czułość", xlab = "czułość", ylab = "precyzja" )
Krzywa dokładności klasyfikacji
Krzywa dokładności klasyfikacji przedstawia zależność między dokładnością klasyfikacji (odsetkiem obiektów prawidłowo sklasyfikowanych) a przyjętym prawdopodobieństwem odcięcia (prawdopodobieństwo powyżej którego przypisywana jest klasa złych kredytów).
m1_perf_acc <- performance(m1_pred, measure = "acc")
plot(m1_perf_acc, main="Dokładność klasyfikacji", xlab = "poziom p", ylab = "dokładnośc")
Ocena jakości modelu przy różnym poziomie prawdopodobieństwa
Poniżej przedstawiono obliczone wskaźniki jakości modelu przy różnym poziomie p.
- dokładność - procent kredytów poprawnie sklasyfikowanych;
- precyzja - jaka część kredytów sklasyfikowanych jako dobre są rzeczywiście dobre;
- czułość - jaka część rzeczywiście dobrych kredytów została sklasyfikowana jako dobre;
- miaraF - średnia harmoniczna z precyzji i czułości;
- specyficzność - jaka część rzeczywiście złych kredytów została sklasyfikowana jako złe;
- index Youdena - obliczany jako: \(czułość + specyficzność -1\) - a maksymalna wartość wskaźnika może służyć jako kryterium wyboru optymalnego punktu odcięcia.
ocena<-data.frame(matrix(ncol = 7, nrow = 0))
colnames(ocena)<-c("prawd", "dokladnosc", "precyzja", "czulosc", "miaraF", "specyficznosc", "Youden_Index")
for (p in seq(0.3, 0.7, 0.05)) {
prawd<-p
data_test$predict<-ifelse(data_test$model>p, "zly", "dobry")
dokladnosc<-1-(sum(data_test$predict != data_test$jakosc) / nrow(data_test))
macierz_bledy <- table(data_test$predict, data_test$jakosc, dnn=c("Predykcja","Rzeczywiste"))
precyzja <- round(
macierz_bledy[1,1]/(macierz_bledy[1,1]+macierz_bledy[1,2]),3)
czulosc <- round(
macierz_bledy[1,1]/(macierz_bledy[1,1]+macierz_bledy[2,1]), 3)
miaraF <-round(
2/(1/precyzja+1/czulosc), 3)
specyficznosc<- round(
macierz_bledy[2,2]/(macierz_bledy[2,2]+macierz_bledy[1,2]),3)
Youden_Index <-round(
czulosc+specyficznosc-1,3)
ocena<-rbind(ocena, data.frame(prawd, dokladnosc, precyzja, czulosc, miaraF, specyficznosc, Youden_Index))
}
ocena## prawd dokladnosc precyzja czulosc miaraF specyficznosc Youden_Index
## 1 0.30 0.720 0.832 0.734 0.780 0.691 0.425
## 2 0.35 0.720 0.811 0.763 0.786 0.630 0.393
## 3 0.40 0.736 0.808 0.799 0.803 0.605 0.404
## 4 0.45 0.744 0.803 0.822 0.812 0.580 0.402
## 5 0.50 0.744 0.787 0.852 0.818 0.519 0.371
## 6 0.55 0.760 0.779 0.899 0.835 0.469 0.368
## 7 0.60 0.760 0.766 0.929 0.840 0.407 0.336
## 8 0.65 0.744 0.751 0.929 0.831 0.358 0.287
## 9 0.70 0.716 0.723 0.941 0.818 0.247 0.188