#Biblioteki
library(pROC)
## Type 'citation("pROC")' for a citation.
##
## Dołączanie pakietu: 'pROC'
## Następujące obiekty zostały zakryte z 'package:stats':
##
## cov, smooth, var
#Dodanie danych
dane <- read.csv("C:/Users/grzeg/OneDrive/Pulpit/Projekt_27/german.data-numeric.csv", header = FALSE, sep = " ")
#Dodanie nazw zmiennych
colnames(dane) <- c(
"Stan_konta", "Czas_kredytu", "Historia_kredytowa", "Cel_kredytu", "Kwota_kredytu",
"Oszczednosci", "Zatrudnienie", "Rata_procent", "Status_osobisty", "Inni_gwaranci",
"Czas_mieszkania", "Wlasnosc", "Wiek", "Inne_plany_ratalne", "Mieszkanie",
"Liczba_kredytow", "Praca", "Liczba_osob", "Telefon", "Pracownik_zagraniczny",
"ID1", "ID2", "ID3","ID4", "Class")
#Usuwam zmienne, które nie były opisane
dane <- dane[, !(names(dane) %in% c("ID1", "ID2", "ID3","ID4"))]
summary(dane)
## Stan_konta Czas_kredytu Historia_kredytowa Cel_kredytu
## Min. :1.000 Min. : 4.0 Min. :0.000 Min. : 2.00
## 1st Qu.:1.000 1st Qu.:12.0 1st Qu.:2.000 1st Qu.: 14.00
## Median :2.000 Median :18.0 Median :2.000 Median : 23.00
## Mean :2.577 Mean :20.9 Mean :2.545 Mean : 32.71
## 3rd Qu.:4.000 3rd Qu.:24.0 3rd Qu.:4.000 3rd Qu.: 40.00
## Max. :4.000 Max. :72.0 Max. :4.000 Max. :184.00
## Kwota_kredytu Oszczednosci Zatrudnienie Rata_procent
## Min. :1.000 Min. :1.000 Min. :1.000 Min. :1.000
## 1st Qu.:1.000 1st Qu.:3.000 1st Qu.:2.000 1st Qu.:2.000
## Median :1.000 Median :3.000 Median :3.000 Median :3.000
## Mean :2.105 Mean :3.384 Mean :2.682 Mean :2.845
## 3rd Qu.:3.000 3rd Qu.:5.000 3rd Qu.:3.000 3rd Qu.:4.000
## Max. :5.000 Max. :5.000 Max. :4.000 Max. :4.000
## Status_osobisty Inni_gwaranci Czas_mieszkania Wlasnosc
## Min. :1.000 Min. :19.00 Min. :1.000 Min. :1.000
## 1st Qu.:1.000 1st Qu.:27.00 1st Qu.:3.000 1st Qu.:1.000
## Median :2.000 Median :33.00 Median :3.000 Median :1.000
## Mean :2.358 Mean :35.55 Mean :2.675 Mean :1.407
## 3rd Qu.:3.000 3rd Qu.:42.00 3rd Qu.:3.000 3rd Qu.:2.000
## Max. :4.000 Max. :75.00 Max. :3.000 Max. :4.000
## Wiek Inne_plany_ratalne Mieszkanie Liczba_kredytow
## Min. :1.000 Min. :1.000 Min. :1.000 Min. :0.000
## 1st Qu.:1.000 1st Qu.:1.000 1st Qu.:1.000 1st Qu.:0.000
## Median :1.000 Median :1.000 Median :1.000 Median :0.000
## Mean :1.155 Mean :1.404 Mean :1.037 Mean :0.234
## 3rd Qu.:1.000 3rd Qu.:2.000 3rd Qu.:1.000 3rd Qu.:0.000
## Max. :2.000 Max. :2.000 Max. :2.000 Max. :1.000
## Praca Liczba_osob Telefon Pracownik_zagraniczny
## Min. :0.000 Min. :0.000 Min. :0.000 Min. :0.000
## 1st Qu.:0.000 1st Qu.:1.000 1st Qu.:0.000 1st Qu.:0.000
## Median :0.000 Median :1.000 Median :0.000 Median :0.000
## Mean :0.103 Mean :0.907 Mean :0.041 Mean :0.179
## 3rd Qu.:0.000 3rd Qu.:1.000 3rd Qu.:0.000 3rd Qu.:0.000
## Max. :1.000 Max. :1.000 Max. :1.000 Max. :1.000
## Class
## Min. :1.0
## 1st Qu.:1.0
## Median :1.0
## Mean :1.3
## 3rd Qu.:2.0
## Max. :2.0
#Sprawdzam czy są braki:
sum(is.na(dane))
## [1] 0
#Nie ma
#Przekształcenia wartości ostatniej zmiennej
dane$Class <- ifelse(dane$Class == 1, 0, 1)
#Ustawiam ziarno z nr indeksu
set.seed(145290)
#Losowy podział zbioru na zbiory: uczący, walidacyjny i testowy w proporcjach odpowiednio 40%, 30% i 30%
n <- nrow(dane)
indeksy <- sample(1:n)
#Określenie liczby wierszy w poszczególnych zbiorach
n_uczacy <- floor(0.4 * n)
n_walidacyjny <- floor(0.3 * n)
n_testowy <- n - n_uczacy - n_walidacyjny
#Wyznaczenie indeksów dla każdego zbioru
uczacy_idx <- indeksy[1:n_uczacy]
walidacyjny_idx <- indeksy[(n_uczacy + 1):(n_uczacy + n_walidacyjny)]
testowy_idx <- indeksy[(n_uczacy + n_walidacyjny + 1):n]
#Tworzenie zbiorów
uczacy <- dane[uczacy_idx, ]
walidacyjny <- dane[walidacyjny_idx, ]
testowy <- dane[testowy_idx, ]
dim(uczacy)
## [1] 400 21
dim(walidacyjny)
## [1] 300 21
dim(testowy)
## [1] 300 21
#Budowa modeli:
# Model 1
model1 <- glm(Class ~ Czas_kredytu + Kwota_kredytu + Rata_procent + Wiek + Liczba_kredytow, data = uczacy, family = binomial)
summary(model1)
##
## Call:
## glm(formula = Class ~ Czas_kredytu + Kwota_kredytu + Rata_procent +
## Wiek + Liczba_kredytow, family = binomial, data = uczacy)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.103279 0.526525 -2.095 0.036136 *
## Czas_kredytu 0.033451 0.009748 3.432 0.000600 ***
## Kwota_kredytu -0.287775 0.080708 -3.566 0.000363 ***
## Rata_procent -0.093258 0.102224 -0.912 0.361613
## Wiek 0.197074 0.314150 0.627 0.530446
## Liczba_kredytow 0.651339 0.251936 2.585 0.009729 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 493.67 on 399 degrees of freedom
## Residual deviance: 461.40 on 394 degrees of freedom
## AIC: 473.4
##
## Number of Fisher Scoring iterations: 4
# Model 2
model2 <- glm(Class ~ Stan_konta + Historia_kredytowa + Cel_kredytu + Zatrudnienie + Kwota_kredytu, data = uczacy, family = binomial)
summary(model2)
##
## Call:
## glm(formula = Class ~ Stan_konta + Historia_kredytowa + Cel_kredytu +
## Zatrudnienie + Kwota_kredytu, family = binomial, data = uczacy)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.658058 0.587015 2.825 0.00473 **
## Stan_konta -0.575667 0.102573 -5.612 2e-08 ***
## Historia_kredytowa -0.380274 0.121177 -3.138 0.00170 **
## Cel_kredytu 0.010937 0.004307 2.540 0.01110 *
## Zatrudnienie -0.038831 0.173459 -0.224 0.82286
## Kwota_kredytu -0.218692 0.083110 -2.631 0.00850 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 493.67 on 399 degrees of freedom
## Residual deviance: 416.72 on 394 degrees of freedom
## AIC: 428.72
##
## Number of Fisher Scoring iterations: 4
#W modulu 1 Rata_procent i Wiek są niesitotne statystycznie, a w modelu 2 Zatrudnienie
#Budowa modeli - usunięcie z nich zmiennych nieistotnych statystycznie:
# Model 1
model1 <- glm(Class ~ Czas_kredytu + Kwota_kredytu + Liczba_kredytow, data = uczacy, family = binomial)
summary(model1)
##
## Call:
## glm(formula = Class ~ Czas_kredytu + Kwota_kredytu + Liczba_kredytow,
## family = binomial, data = uczacy)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.12106 0.29338 -3.821 0.000133 ***
## Czas_kredytu 0.03267 0.00972 3.361 0.000776 ***
## Kwota_kredytu -0.29130 0.08056 -3.616 0.000299 ***
## Liczba_kredytow 0.65131 0.25069 2.598 0.009374 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 493.67 on 399 degrees of freedom
## Residual deviance: 462.55 on 396 degrees of freedom
## AIC: 470.55
##
## Number of Fisher Scoring iterations: 4
# Model 2
model2 <- glm(Class ~ Stan_konta + Historia_kredytowa + Cel_kredytu + Kwota_kredytu, data = uczacy, family = binomial)
summary(model2)
##
## Call:
## glm(formula = Class ~ Stan_konta + Historia_kredytowa + Cel_kredytu +
## Kwota_kredytu, family = binomial, data = uczacy)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.563069 0.404839 3.861 0.000113 ***
## Stan_konta -0.575430 0.102515 -5.613 1.99e-08 ***
## Historia_kredytowa -0.383632 0.120295 -3.189 0.001427 **
## Cel_kredytu 0.010966 0.004306 2.547 0.010874 *
## Kwota_kredytu -0.219355 0.083074 -2.640 0.008279 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 493.67 on 399 degrees of freedom
## Residual deviance: 416.77 on 395 degrees of freedom
## AIC: 426.77
##
## Number of Fisher Scoring iterations: 4
#Prezentacja zdolności predykcyjnej modeli na zbiorze uczącym i walidacyjnym
#Model 1
prob_uczacy_m1 <- predict(model1, newdata = uczacy, type = "response")
prob_walidacyjny_m1 <- predict(model1, newdata = walidacyjny, type = "response")
#Model 2
prob_uczacy_m2 <- predict(model2, newdata = uczacy, type = "response")
prob_walidacyjny_m2 <- predict(model2, newdata = walidacyjny, type = "response")
#ROC dla modelu 1
roc_uczacy_m1 <- roc(uczacy$Class, prob_uczacy_m1)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
roc_walidacyjny_m1 <- roc(walidacyjny$Class, prob_walidacyjny_m1)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
#ROC dla modelu 2
roc_uczacy_m2 <- roc(uczacy$Class, prob_uczacy_m2)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
roc_walidacyjny_m2 <- roc(walidacyjny$Class, prob_walidacyjny_m2)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
#AUC
auc(roc_uczacy_m1)
## Area under the curve: 0.6689
auc(roc_walidacyjny_m1)
## Area under the curve: 0.6701
auc(roc_uczacy_m2)
## Area under the curve: 0.7634
auc(roc_walidacyjny_m2)
## Area under the curve: 0.7565
#ROC-wykres porównawczy dla zbioru uczącego model1
plot(roc_uczacy_m1, col = "blue", lwd = 2, main = "ROC – Zbiór uczący")
lines(roc_uczacy_m2, col = "red", lwd = 2)
legend("bottomright", legend = c("Model 1", "Model 2"), col = c("blue", "red"), lwd = 2)
#ROC-wykres porównawczy dla zbioru walidacyjnego model2
plot(roc_walidacyjny_m1, col = "blue", lwd = 2, main = "ROC – Zbiór walidacyjny")
lines(roc_walidacyjny_m2, col = "red", lwd = 2)
legend("bottomright", legend = c("Model 1", "Model 2"), col = c("blue", "red"), lwd = 2)
#KOM: Model 2 ma wyraźnie lepszą zdolność predykcyjną!
#Porównanie zbioru uczącego i walidacyjnego dla obu modeli:
#ROC dla Modelu 1 uczący/walidacyjny
plot(roc_uczacy_m1, col = "blue", lwd = 2, main = "ROC Model 1 – Uczący vs Walidacyjny")
lines(roc_walidacyjny_m1, col = "green", lwd = 2)
legend("bottomright",
legend = c("Zbiór uczący", "Zbiór walidacyjny"),
col = c("blue", "green"),
lwd = 2)
#ROC dla Modelu 2 uczący/walidacyjny
plot(roc_uczacy_m2, col = "red", lwd = 2, main = "ROC Model 2 – Uczący vs Walidacyjny")
lines(roc_walidacyjny_m2, col = "orange", lwd = 2)
legend("bottomright",
legend = c("Zbiór uczący", "Zbiór walidacyjny"),
col = c("red", "orange"),
lwd = 2)
#Prezentacja zdolności predykcyjnej najlepszego moidelu (modelu2) na zbiorze walidacyjnym i testowym
prob_testowy_m2 <- predict(model2, newdata = testowy, type = "response")
#ROC dla testowego
roc_testowy_m2 <- roc(testowy$Class, prob_testowy_m2)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
#AUC dla testowego i walidacyjnego
auc_testowy_m2 <- auc(roc_testowy_m2)
auc_walidacyjny_m2 <- auc(roc_walidacyjny_m2)
auc(roc_testowy_m2)
## Area under the curve: 0.7768
auc(roc_walidacyjny_m2)
## Area under the curve: 0.7565
#Wykres ROC dla Modelu 2 testowy/walidacyjny
plot(roc_walidacyjny_m2, col = "blue", lwd = 2, main = "ROC – Model 2 (Walidacyjny vs Testowy)")
lines(roc_testowy_m2, col = "red", lwd = 2)
legend("bottomright",
legend = c(paste("Walidacyjny (AUC =", round(auc_walidacyjny_m2,4), ")"),paste("Testowy (AUC =", round(auc_testowy_m2,4), ")")),
col = c("blue","red"),
lwd = 2)
#Krzywa kosztów złych decyzji
#Parametry
koszt_FP <- 150
koszt_FN <- 1000
#Progi odcięcia
progi <- seq(0.5, 0.9, by = 0.01)
#Predykcje prawdopodobieństwa na zbiorze walidacyjnym
prawd_walidacyjny <- predict(model2, newdata = walidacyjny, type = "response")
prawdziwa_klasa <- walidacyjny$Class
#Funkcja do obliczania kosztu
koszt_calkowity <- sapply(progi, function(prog){
pred_klasa <- ifelse(prawd_walidacyjny >= prog, 1, 0)
FP <- sum(pred_klasa == 1 & prawdziwa_klasa == 0)
FN <- sum(pred_klasa == 0 & prawdziwa_klasa == 1)
koszt <- FP * koszt_FP + FN * koszt_FN
return(koszt)
})
#Wykresu
plot(progi, koszt_calkowity, type = "l", col = "darkred", lwd = 2,
xlab = "Próg odcięcia", ylab = "Koszt złych decyzji [PLN]",
main = "Krzywa kosztów złych decyzji – Model 2 (Walidacyjny)")
grid()
#Optymalny próg odcięcia, dla którego koszty złych decyzji będą najniższe.
#Znalezienie
min_koszt <- min(koszt_calkowity)
optymalny_prog <- progi[which.min(koszt_calkowity)]
#Wyniku
cat("Minimalny koszt złych decyzji:", min_koszt, "PLN\n")
## Minimalny koszt złych decyzji: 54600 PLN
cat("Optymalny próg odcięcia:", optymalny_prog, "\n")
## Optymalny próg odcięcia: 0.52
#Wykres
plot(progi, koszt_calkowity, type = "l", col = "darkred", lwd = 2,
xlab = "Próg odcięcia", ylab = "Koszt złych decyzji [PLN]",
main = "Krzywa kosztów złych decyzji – Model 2 (Walidacyjny)")
grid()
abline(v = optymalny_prog, col = "blue", lwd = 2, lty = 2)
text(optymalny_prog, min_koszt + 1000, paste("Optymalny próg =", round(optymalny_prog,2)), col="blue")
#Predykcje prawdopodobieństwa na zbiorze testowym dla najlepszego modelu (model2)
prawd_testowy <- predict(model2, newdata = testowy, type = "response")
prawdziwa_klasa_test <- testowy$Class
#Przypisanie klas według progu
pred_klasa_test <- ifelse(prawd_testowy >= optymalny_prog, 1, 0)
#Obliczenie liczby FP i FN
FP <- sum(pred_klasa_test == 1 & prawdziwa_klasa_test == 0)
FN <- sum(pred_klasa_test == 0 & prawdziwa_klasa_test == 1)
#Całkowity koszt
koszt_calkowity_test <- FP * koszt_FP + FN * koszt_FN
#Wyniki
cat("Koszt złych decyzji na zbiorze testowym przy progu", optymalny_prog, ":", koszt_calkowity_test, "PLN\n")
## Koszt złych decyzji na zbiorze testowym przy progu 0.52 : 57150 PLN
cat("Liczba fałszywych pozytywów (FP):", FP, "\n")
## Liczba fałszywych pozytywów (FP): 21
cat("Liczba fałszywych negatywów (FN):", FN, "\n")
## Liczba fałszywych negatywów (FN): 54