#  READING CREDSCO.TXT

dd <- read.table("C:/Users/aluja/Documents/txt/docent/Curs MD english/R/credsco.txt",header=T)

dim(dd)

names(dd)

# FIRST CONTROL: MAXIMS, MINIMUMS, ... (FINDING ERRORS AND OUTLIERS)
# WHICH RESPONSE?
# WHICH ARE CATEGORICAL AND WHICH ARE CONTINUOUS
# ARE THERE MISSINGS?

summary(dd)

# DEALING WITH MISSINGS: ELIMINATE FROM DICTAMENT MANCANT, VIVENDA MANCANT, ESTAT CIVIL MANCANT I TIPUS DE FEINA MANCANT (VALORS MANCANTS DE LES VARAIBLES CATEGORIQUES, son pocs)

table(dd[,1]==0)
table(dd[,3]==0)
table(dd[,6]==0)
table(dd[,8]==0)

dd <- dd[dd[,1] != 0 & dd[,3] != 0 & dd[,6] != 0 & dd[,8] != 0,]
dim(dd)


# LOOK FOR MISSING VALUES AMONG THE CONTINUOUS VARIABLES (99999999 MISSING CODE)

# LETS WORK WITH DIRECTLY WITH THE VARIABLES OF THE DATA FRAME

objects()

attach(dd)

objects()

table(Ingressos == 99999999)

table(Patrimoni == 99999999)

table(Carrecs_pat == 99999999)

table(Antig_feina == 0)

table(Ingressos == 0)

# LES VAS. CONTINUES TENEN MASSES VALORS MANCANTS, NO ELS PODEM ELIMINAR COM SI RES
# HEM DE FER UN TRACTAMENT DE LES DADES MANCANTS
# PASSEM ELS VALORS MANCANTS A NA, INCLOSOS ELS 0 DELS INGRESSOS

Ingressos[Ingressos == 99999999 | Ingressos == 0] <- NA
Patrimoni[Patrimoni == 99999999] <- NA
Carrecs_pat[Carrecs_pat == 99999999] <- NA

# WARNING: NOW dd[,10] Ingressos DOESNT HAVE THE SAME CONTENT

summary(dd[,10])
summary(Ingressos)


# IMPUTATION By THE 1NN

library(class)

# FOR EVERY INDIVIDUAL WITH INGRESSOS MISSING WE LOOK FOR THE MOST SIMILAR INDIVIDUAL ACCORDING THE REMAINING VARIABLES AND WE COPY THE VALUE OF INGRESSOS ON THE FIRST 

aux = dd[,-10]
dim(aux)
aux1 = aux[!is.na(Ingressos),]
dim(aux1)
aux2 = aux[is.na(Ingressos),]
dim(aux2)
knn.ing = knn(aux1,aux2,Ingressos[!is.na(Ingressos)])   # NEITHER AUX1, AUX2 CAN CONTAIN NAs)
Ingressos[is.na(Ingressos)] = knn.ing

# IMPUTATION OF PATRIMONI
aux = dd[,-11]
aux1 = aux[!is.na(Patrimoni),]
aux2 = aux[is.na(Patrimoni),]
knn.pat = knn(aux1,aux2,Patrimoni[!is.na(Patrimoni)])
Patrimoni[is.na(Patrimoni)] = knn.pat

# IMPUTATION OF CARRECS PATRIMONIALS
aux = dd[,-12]
aux1 = aux[!is.na(Carrecs_pat),]
aux2 = aux[is.na(Carrecs_pat),]
knn.car = knn(aux1,aux2, Carrecs_pat[!is.na(Carrecs_pat)])
Carrecs_pat[is.na(Carrecs_pat)] = knn.car
Carrecs_pat[Patrimoni==0] <- 0

dd[,10] <- Ingressos
dd[,11] <- Patrimoni
dd[,12] <- Carrecs_pat

# VERIFY AGAIN MAX AND MIN

dim(dd)
summary(dd)


# DECLARE CATEGORICAL 

Dictamen    <- as.factor(Dictamen)
Vivenda     <- as.factor(Vivenda)
Estat_civil <- as.factor(Estat_civil)
Registres   <- as.factor(Registres)
Tipus_feina <- as.factor(Tipus_feina)

levels(Dictamen)
levels(Vivenda)
levels(Estat_civil)
levels(Registres)
levels(Tipus_feina)


levels(Dictamen) <- c("positiu","negatiu")
levels(Vivenda) <- c("lloguer","escriptura","contr_privat","ignora_cont","pares","altres viv")
levels(Estat_civil) <- c("solter","casat","vidu","separat","divorciat")
levels(Registres) <- c("reg_no","reg_si")
levels(Tipus_feina) <- c("fixe","temporal","autonom","altres sit")

# WARNING

is.factor(Dictamen)
is.factor(dd[,1])

# DERIVATION OF NEW VARIABLES: “FEATURE EXTRACTION”

# RATIO OF FINANCEMENT 

Rati_fin = 100*Import_sol/Preu_finan

hist(Rati_fin)

# CAPACITY TO SAVE

Estalvi <- (Ingressos-Despeses-(Carrecs_pat/100))/(Import_sol/Plaç)

hist(Estalvi)

# RECODING THE CONTINUOUS VARIABLES TO CATEGORICAL (TO ASSESS NON LINEARITIES)

antigR <- cut(Antig_feina, breaks=c(-1,1,3,8,14,99))
plaçR <- cut(Plaç, breaks=c(0,12,24,36,48,99))
edatR <- cut(Edat, breaks=c(0,25,30,40,50,99))
despesesR <- cut(Despeses, breaks=c(0,40,50,60,80,9999))
ingressosR <- cut(Ingressos, breaks=c(0,80,110,140,190,9999))
patrimoniR <- cut(Patrimoni, breaks=c(-1,0,3000,5000,8000,999999))
carrecsR <- cut(Carrecs_pat, breaks=c(-1,0,500,1500,2500,999999))
importR <- cut(Import_sol, breaks=c(0,600,900,1100,1400,99999))
preuR <- cut(Preu_finan, breaks=c(0,1000,1300,1500,1800,99999))
ratfinR <- cut(Rati_fin, breaks=c(0,50,70,80,90,100))
estalviR <- cut(Estalvi, breaks=c(-99,0,2,4,6,99))

levels(antigR) <- paste("Antig",levels(antigR))
levels(plaçR) <- paste("Plaç",levels(plaçR))
levels(edatR) <- paste("Edat",levels(edatR))
levels(despesesR) <- paste("Desp",levels(despesesR))
levels(ingressosR) <- paste("Ingr",levels(ingressosR))
levels(patrimoniR) <- paste("Patr",levels(patrimoniR))
levels(carrecsR) <- paste("Carr",levels(carrecsR))
levels(importR) <- paste("Import",levels(importR))
levels(preuR) <- paste("Preu",levels(preuR))
levels(ratfinR) <- paste("Ratfin",levels(ratfinR))
levels(estalviR) <- paste("Estalv",levels(estalviR))


# DECISION TREES   CART

library(rpart)

dtot <- data.frame(Dictamen,Antig_feina,Plaç,Edat,Despeses,Ingressos,Patrimoni,Carrecs_pat,
Import_sol,Preu_finan,Rati_fin,Estalvi,edatR,antigR,plaçR,despesesR,ingressosR,patrimoniR,
carrecsR,importR,preuR,ratfinR,estalviR,Vivenda,Estat_civil,Registres,Tipus_feina)

n <- nrow(dtot)

# HOLDOUT OF A 1/3 OF DATA TO ESTIMATE THE MISCLASSIFICATION PROB. OF THE CHOSEN TREE 

learn <- sample(1:n, round(0.67*n))

nlearn <- length(learn)
ntest <- n - nlearn

# SELECTION OF MODEL BY CROSSVALIDATION

# FIRST WE BUILD A MAXIMAL TREE (WITH LOW VALUE OF cp) AND EQUIPROBABILITY OF CLASSES
 
ct1 = rpart(Dictamen~Antig_feina+Plaç+Edat+Despeses+Ingressos+Patrimoni+Carrecs_pat+Import_sol+
Preu_finan+Rati_fin+Estalvi+edatR+antigR+plaçR+despesesR+ingressosR+patrimoniR+carrecsR+
importR+preuR+ratfinR+estalviR+Vivenda+Estat_civil+Registres+Tipus_feina, data=dtot[learn,], parms=list(prior=c(.50,.50), split='gini'), control=rpart.control(cp=0.001, xval=10, maxdepth=15))


# FIRST LETS TAKE A LOOK ON THIS OBJECT ct1

print(ct1)  # COULD YOU GUESS WHICH INFORMATION DELIVERS? COULD YOU FOLLOW THE SPLITTING PROCESS?

# LETS TAKE A LOOK VISUALLY

plot(ct1)

# NOW WE WILL SELECT THE OPTIMAL TREE (WITH LOWER CROSSVALIDATION COST)
# WE LOOK AT THE COST OF THE TREES WITH INCREASING cp
# WE SELECT THE TREE WITH LOWEST CROSS-VALIDATED COST (Min XERROR)

ct1$cptable

# THE OPTIMAL TREE CORRESPONDS TO A cp =  0.003

# NOW LETS OBTAIN THIS TREE

ct1 = rpart(Dictamen~Antig_feina+Plaç+Edat+Despeses+Ingressos+Patrimoni+Carrecs_pat+Import_sol+
Preu_finan+Rati_fin+Estalvi+edatR+antigR+plaçR+despesesR+ingressosR+patrimoniR+carrecsR+
importR+preuR+ratfinR+estalviR+Vivenda+Estat_civil+Registres+Tipus_feina, data=dtot[learn,], parms=list(prior=c(.50,.50), split='gini'), control=rpart.control(cp= 0.00285, xval=0, maxdepth=15))

# MIREM L'ARBRE QUE HEM FET

ct1
plot(ct1)
text(ct1, use.n=T)

summary(ct1)

# CACULATE THE ERROR RATE IN THE LEARNING SAMPLE

ct1l=predict(ct1,data=dtot[learn,])

# LETS HAVE A LOOK ON THE PREDICTIONS

ct1l[1:10,]

# ERROR RATE IN THE LEARNING SAMPLE

ct1lp=NULL
ct1lp[ct1l[,1]<0.5]="pred_neg.learn"
ct1lp[ct1l[,1]>=0.5]="pred_pos.learn"
table(Dictamen[learn],ct1lp)
aa <- table(Dictamen[learn],ct1lp)
error_rate.learn <- 100*(aa[1,1]+aa[2,2])/nlearn 
error_rate.learn


# ERROR RATE IN THE TEST SAMPLE

ct1t=predict(ct1,newdata=dtot[-learn,])
ct1tp=NULL
ct1tp[ct1t[,1]<0.5]="pred_neg.test"
ct1tp[ct1t[,1]>=0.5]="pred_pos.test"
table(Dictamen[-learn],ct1tp)
aa <- table(Dictamen[-learn],ct1tp)
error_rate.test <- 100*(aa[1,1]+aa[2,2])/ntest 
error_rate.test


# NOW WE WILL REPEAT THE SAME BUT CHANGING THE cp (=0.002)

ct2 = rpart(Dictamen~Antig_feina+Plaç+Edat+Despeses+Ingressos+Patrimoni+Carrecs_pat+Import_sol+
Preu_finan+Rati_fin+Estalvi+edatR+antigR+plaçR+despesesR+ingressosR+patrimoniR+carrecsR+
importR+preuR+ratfinR+estalviR+Vivenda+Estat_civil+Registres+Tipus_feina, data=dtot[learn,], parms=list(prior=c(.50,.50), split='gini'), control=rpart.control(cp= 0.002, xval=0, maxdepth=15))


# CACULATE THE ERROR RATE IN THE LEARNING SAMPLE

ct2l=predict(ct2,data=dtot[learn,])


# ERROR RATE IN THE LEARNING SAMPLE

ct2lp=NULL
ct2lp[ct2l[,1]<0.5]="pred_neg.learn"
ct2lp[ct2l[,1]>=0.5]="pred_pos.learn"
table(Dictamen[learn],ct2lp)
aa <- table(Dictamen[learn],ct2lp)
error_rate2.learn <- 100*(aa[1,1]+aa[2,2])/nlearn 
error_rate2.learn


# ERROR RATE IN THE TEST SAMPLE

ct2t=predict(ct2,newdata=dtot[-learn,])
ct2tp=NULL
ct2tp[ct2t[,1]<0.5]="pred_neg.test"
ct2tp[ct2t[,1]>=0.5]="pred_pos.test"
table(Dictamen[-learn],ct2tp)
aa <- table(Dictamen[-learn],ct2tp)
error_rate2.test <- 100*(aa[1,1]+aa[2,2])/ntest 
error_rate2.test

# CONCENTRATION CURVE

# THE POSITIVE PREDICTIONS ON THE TEST SAMPLE

pred.test <- ct1t[,1]

# THE TOTAL NUMBER OF TEST INDIVIDUALS IN EACH LEAVE 

totn <- table(-pred.test)/ntest
ac_totn <- 100*cumsum(as.numeric(totn))

# RANKING THE PREDICTIONS

rank_pred.test <- rank(pred.test)

# COMPUTING HOW MANY POSITIVE ARE IN EACH LEAVE (LEAVE = EQUAL RANK INDIVIDUALS)


Dict.test <- Dictamen[-learn]
table(Dict.test)
npos <- table(Dict.test)[1]

ac_true.pos <- 100*cumsum(rev(as.numeric(tapply(Dict.test=="positiu",rank_pred.test,sum))))/npos

# PLOT

plot(ac_totn,ac_true.pos,type="l", main="Concentration Curve")
lines(ac_totn, ac_totn, col="red")


# ANOTHER WAY OF DOING THE CONCENTRATION CURVE

ac_tot <- 100*(1:ntest)/ntest

ord_pred.test <- order(pred.test,decreasing=T)

ac_pos_test <- 100*cumsum(Dict.test[ord_pred.test]=="positiu")/npos

plot(ac_tot,ac_pos_test,type="l", main="Concentration Curve")
lines(ac_tot, ac_tot, col="red")


# ROC CURVE

nneg <- ntest - npos

ac_fals.pos <- 100*cumsum(rev(as.numeric(tapply(Dict.test=="negatiu",rank_pred.test,sum))))/nneg

# PLOT

plot(ac_fals.pos,ac_true.pos,type="l", main="ROC curve")
lines(ac_fals.pos, ac_fals.pos, col="red")

# OR

ac_neg_test <- 100*cumsum(Dict.test[ord_pred.test]=="negatiu")/nneg

plot(ac_neg_test,ac_pos_test,type="l", main="Concentration Curve")
lines(ac_neg_test, ac_neg_test, col="red")

# NOW TRY USING THE INFORMATION AS SPLIT CRITERION