file.name<-'http://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data-numeric' credit.approval<-read.table(file.name) names(credit.approval)[25]<-'target' ######Podział na zbiory dziel<-function(dane, frakcje) { frakcja_train<-frakcje[1] frakcja_valid<-frakcje[2] frakcja_test<-frakcje[3] losuj<-runif(dim(dane)[1]) assign<-ifelse(losuj<frakcja_train, 'train', ifelse(losuj<frakcja_train+frakcja_valid, 'valid', 'test')) assign<-factor(assign, levels=c('train', 'valid', 'test')) return(assign) } dziel2<-function(dane, frakcje) { losuj<-runif(dim(dane)[1]) frakcja_train<-quantile(losuj, frakcje[1]) frakcja_valid<-quantile(losuj, frakcje[2]+frakcje[1]) # frakcja_test<-quantile(losuj, frakcje[3]) assign<-ifelse(losuj<frakcja_train, 'train', ifelse(losuj<frakcja_valid, 'valid', 'test')) assign<-factor(assign, levels=c('train', 'valid', 'test')) return(assign) } table(credit.approval$target) credit.approval$target<-credit.approval$target-1 przydzial<-dziel2(credit.approval, c(.6,.2,.2)) dane$train<-credit.approval[which(przydzial=='train'),] dane$valid<-credit.approval[which(przydzial=='valid'),] dane$test<-credit.approval[which(przydzial=='test'),] dane<-split(credit.approval, przydzial) ols2<-glm(target~., data=dane$train) logit2<-glm(target~., data=dane$train, family=binomial("logit")) probit2<-glm(target~., data=dane$train, family=binomial("probit")) ols2_step<-step(ols2) logit2_step<-step(logit2) probit2_step<-step(probit2) predict(ols2, newdata=dane$valid, type='response') predict(logit2, newdata=dane$valid, type='response') predict(probit2, newdata=dane$valid, type='response') koszt2 <- function(prog, model, data) { decyzja <- ifelse(predict(model, newdata=data, type='response')>prog,1,0) decyzja<-factor(decyzja, levels=c(0,1)) pt <- prop.table(table(decyzja,data$target)) pt[1,2] + pt[2,1] } optymalne_odciecie<-c(optimize(koszt2, c(0,1), model=ols2, data=dane$train)$minimum, optimize(koszt2, c(0,1), model=ols2_step, data=dane$train)$minimum, optimize(koszt2, c(0,1), model=logit2, data=dane$train)$minimum, optimize(koszt2, c(0,1), model=logit2_step, data=dane$train)$minimum, optimize(koszt2, c(0,1), model=probit2, data=dane$train)$minimum, optimize(koszt2, c(0,1), model=probit2_step, data=dane$train)$minimum) koszt_train<-c( koszt2(optymalne_odciecie[1], ols2, data=dane$train), koszt2(optymalne_odciecie[2], ols2_step, data=dane$train), koszt2(optymalne_odciecie[3], logit2, data=dane$train), koszt2(optymalne_odciecie[4], logit2_step, data=dane$train), koszt2(optymalne_odciecie[5], probit2, data=dane$train), koszt2(optymalne_odciecie[6], probit2_step, data=dane$train)) koszt_valid<-c( koszt2(optymalne_odciecie[1], ols2, data=dane$valid), koszt2(optymalne_odciecie[2], ols2_step, data=dane$valid), koszt2(optymalne_odciecie[3], logit2, data=dane$valid), koszt2(optymalne_odciecie[4], logit2_step, data=dane$valid), koszt2(optymalne_odciecie[5], probit2, data=dane$valid), koszt2(optymalne_odciecie[6], probit2_step, data=dane$valid)) credit.approval podsumowanie<-data.frame(model=c("OLS","OLS Step", "Logit", "Logit Step", "Probit", "Probit Step" ),cutoff=optymalne_odciecie, KosztTrain=koszt_train, KosztValid=koszt_valid) set.seed(38166) zmienne_od_czapy<-replicate(200,rnorm(dim(dane$train)[1])) dane_od_czapy<-cbind(dane$train, zmienne_od_czapy) klej4<-function(x) {paste('V', x)} names(dane_od_czapy)[26:45]<-sapply(26:45, klej4) ols_od_czapy<-glm(target~., data=dane_od_czapy) opt_od_czapy<-optimize(koszt2, c(0,1), model=ols_od_czapy, data=dane_od_czapy) dane_od_czapy_valid<-cbind(dane$valid, zmienne_od_czapy) names(dane_od_czapy_valid)[26:45]<-sapply(26:45, klej4) Z<-koszt2(opt_od_czapy$minimum, ols_od_czapy, data=dane_od_czapy_valid) podsumowanie<-rbind(podsumowanie, data.frame(model="od czapy", cutoff=opt_od_czapy$minimum, KosztTrain=opt_od_czapy$objective, KosztValid=Z)) podsumowanie ##########KRZYWE LIFT, GAIN, ROC library('ROCR') prognoza<-predict(probit2_step, newdata=dane$valid, type="response") klasa<-ifelse(prognoza>podsumowanie$cutoff[6],1,0) table(prognoza=klasa, realizacja=dane$valid$target) pred <- prediction(prognoza, dane$valid$target) pred_wizard<- prediction(dane$valid$target,dane$valid$target) pred_random<-prediction(runif(length(dane$valid$target),0,1), dane$valid$target) #########TRAIN SET prognoza_train<-predict(probit2_step, newdata=dane$train, type="response") klasa_train<-ifelse(prognoza_train>podsumowanie$cutoff[6],1,0) table(klasa_train, dane$train$target) pred_train <- prediction(prognoza_train, dane$train$target) ###############GAIN krzywa_gain <- performance(pred,"tpr","rpp") krzywa_gain_wizard <- performance(pred_wizard,"tpr","rpp") krzywa_gain_random <- performance(pred_random,"tpr","rpp") krzywa_gain_train <- performance(pred_train,"tpr","rpp") plot(krzywa_gain_train, col="green", lty=2) plot(krzywa_gain, col="green", add=TRUE) plot(krzywa_gain_wizard, add=TRUE, col="blue") plot(krzywa_gain_random, add=TRUE, col="red") klasyfikacja <- function(prog, model, data) { decyzja <- ifelse(predict(model, newdata=data, type='response')>prog,1,0) decyzja<-factor(decyzja, levels=c(0,1)) pt <- (table(prognoza=decyzja,realizacja=data$target)) return(pt) } rpp <- function(prog, model, data) { decyzja <- ifelse(predict(model, newdata=data, type='response')>prog,1,0) decyzja<-factor(decyzja, levels=c(0,1)) pt <- prop.table(table(decyzja,data$target)) pt[2,1] + pt[2,2] } tpr <- function(prog, model, data) { decyzja <- ifelse(predict(model, newdata=data, type='response')>prog,1,0) decyzja<-factor(decyzja, levels=c(0,1)) pt <- prop.table(table(decyzja,data$target)) pt[2,2] / (pt[2,2]+pt[1,2]) } ##############LIFT krzywa_lift <- performance(pred,"lift","rpp") krzywa_lift_wizard <- performance(pred_wizard,"lift","rpp") krzywa_lift_random <- performance(pred_random,"lift","rpp") krzywa_lift_train <- performance(pred_train,"lift","rpp") plot(krzywa_lift, col="green", lwd=3) plot(krzywa_lift_wizard, add=TRUE, col="blue") plot(krzywa_lift_random, add=TRUE, col="red") plot(krzywa_lift_train, add=TRUE, col="orange", lty=3, lwd=1.5) grid() lift <- function(prog, model, data) { decyzja <- ifelse(predict(model, newdata=data, type='response')>prog,1,0) decyzja<-factor(decyzja, levels=c(0,1)) pt <- prop.table(table(decyzja,data$target)) (pt[2,2]/(pt[2,2]+pt[1,2])) / (pt[2,2]+pt[2,1]) } #############ROC krzywa_roc <- performance(pred,"tpr","fpr") krzywa_roc_wizard <- performance(pred_wizard,"tpr","fpr") krzywa_roc_random <- performance(pred_random,"tpr","fpr") krzywa_roc_train <- performance(pred_train,"tpr","fpr") plot(krzywa_roc, col="green") plot(krzywa_roc_train, col="green", add=TRUE, lty=2) plot(krzywa_roc_wizard, add=TRUE, col="blue") plot(krzywa_roc_random, add=TRUE, col="red") grid() fpr <- function(prog, model, data) { decyzja <- ifelse(predict(model, newdata=data, type='response')>prog,1,0) decyzja<-factor(decyzja, levels=c(0,1)) pt <- prop.table(table(decyzja,data$target)) pt[2,1] / (pt[2,1]+pt[1,1]) } ##############PORÓWNYWANIE MODELI prognoza_probit2<-predict(probit2, newdata=dane$valid, type="response") klasa_probit2<-ifelse(prognoza>podsumowanie$cutoff[5],1,0) pred_probit2 <- prediction(klasa_probit2, dane$valid$target) lift_probit2<-performance(pred_probit2,"lift","rpp") prognoza_ols2<-predict(ols2, newdata=dane$valid, type="response") klasa_ols2<-ifelse(prognoza>podsumowanie$cutoff[1],1,0) pred_ols2 <- prediction(klasa_ols2, dane$valid$target) lift_ols2<-performance(pred_ols2,"lift","rpp") plot(krzywa_lift, col="green") plot(lift_ols2, col="blue", add=TRUE) plot(lift_probit2, col="mediumpurple", add=TRUE) krzywa_gain <- performance(pred,"tpr","rpp") krzywa_probit2 <- performance(pred_probit2,"tpr","rpp") krzywa_ols2 <- performance(pred_ols2,"tpr","rpp") plot(krzywa_gain, col="green") plot(krzywa_probit2, col="blue", add=TRUE) plot(krzywa_ols2, col="mediumpurple", add=TRUE) grid() krzywa_roc <- performance(pred,"tpr","fpr") krzywa_roc_probit2 <- performance(pred_probit2,"tpr","fpr") krzywa_roc_ols2 <- performance(pred_ols2,"tpr","fpr") plot(krzywa_roc, col="green") plot(krzywa_roc_probit2, col="blue", add=TRUE) plot(krzywa_roc_ols2, col="mediumpurple", add=TRUE) grid()