#  READING CREDSCO.TXT

dd <- read.table("F:/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))

# 
# LETS TAKE THE FACTORIAL COORDINATES OF THE MCA ANALYSIS TO PERFORM THE CLUSTERING  
# 

# CREATION OF THE DATA FRAME OF CATEGORICAL VARIABLES

dcat <- data.frame(edatR,antigR,plaçR,despesesR,ingressosR,patrimoniR,carrecsR,importR,preuR,ratfinR,estalviR,Vivenda,Estat_civil,Registres,Tipus_feina)


# LOADING ACM FUNTION

source("C:/Users/aluja/Documents/txt/docent/MD/R/acm.r") # LA TROBAREU A LA PAGINA WEB DE L'ASSIGNATURA

ac1 <- acm(dcat)

attributes(ac1)

# HISTOGRAM OF EIGENVALUES OF THE MCA (WHY ARE SO DIFFERENT?)

barplot(ac1$vaps)

# WHAT IS THE SIGNIFICANT DIMENSION IN MCA?

i <- 1

while (ac1$vaps[i] > mean(ac1$vaps)) i = i+1

nd = i-1


#
# CLUSTERING
#

#
# LETS TAKE THE INDIVIDUALS DEFINED BY FACTORIAL COORDINATES ISSUED FROM THE MCA
# 

Psi = ac1$rs[,1:nd]

# KMEANS RUN, BUT HOW MANY CLASSES?

k1 <- kmeans(Psi,5)

attributes(k1)

k1$size

k1$withinss

k1$centers

# LETS COMPUTE THE DECOMPOSITION OF INERTIA

Bss <- sum(rowSums(k1$centers^2)*k1$size)
Bss
Wss <- sum(k1$withinss)
Wss
Tss <- sum(rowSums(Psi^2))
Tss

Bss+Wss

Ib1 <- 100*Bss/(Bss+Wss)
Ib1

# LETS REPEAT THE KMEANS RUN WITH K=5

k2 <- kmeans(Psi,5)
k2$size

Bss <- sum(rowSums(k2$centers^2)*k2$size)
Bss
Wss <- sum(k2$withinss)
Wss

Ib2 <- 100*Bss/(Bss+Wss)
Ib2

# WHY WE HAVE OBTAINED DIFFERENT RESULTS?, AND WHICH RUN IS BETTER?

# NOW TRY K=8

k3 <- kmeans(Psi,100)
k3$size

Bss <- sum(rowSums(k3$centers^2)*k3$size)
Wss <- sum(k3$withinss)

Ib3 <- 100*Bss/(Bss+Wss)
Ib3


# HIERARCHICAL CLUSTERING

d  <- dist(Psi)
h1 <- hclust(d,method="ward")  # NOTICE THE COST
plot(h1)

# BUT WE ONLY NEED WHERE THERE ARE THE LEAPS OF THE HEIGHT

# WHERE ARE THER THE LEAPS? WHERE WILL YOU CUT THE DENDREOGRAM?, HOW MANY CLASSES WILL YOU OBTAIN?

nc = 8

c1 <- cutree(h1,nc)

c1[1:20]

# LETS SEE THE PARTITION VISUALLY

plot(Psi[,1],Psi[,2],col=c1,main="Clustering of credit data in 8 classes")
legend("topleft",c("c1","c2","c3","c4","c5","c6","c7","c8"),pch=1,col=c(1:8))

# LETS SEE THE QUALITY OF THE HIERARCHICAL PARTITION

cdg <- aggregate(as.data.frame(Psi),list(c1),mean)[,2:(nd+1)]
cdg


Bss <- sum(rowSums(cdg^2)*as.numeric(table(c1)))

Ib4 <- 100*Bss/Tss
Ib4

# LETS CONSOLIDATE THE PARTITION

k5 <- kmeans(Psi,centers=cdg)
k5$size

Bss <- sum(rowSums(k5$centers^2)*k5$size)
Wss <- sum(k5$withinss)

Ib5 <- 100*Bss/(Bss+Wss)
Ib5

#
# CLUSTERING OF LARGE DATA SETS
#

# FIRST 2 KMEANS WITH K=14

n1 = 14

k1 <- kmeans(Psi,n1)
k2 <- kmeans(Psi,n1)

table(k2$cluster,k1$cluster)

clas <- (k2$cluster-1)*n1+k1$cluster

freq <- table(clas)
freq[1:10]

# WHAT DO WE HAVE IN VECTOR freq?

cdclas <- aggregate(as.data.frame(Psi),list(clas),mean)[,2:(nd+1)]

# SECOND HIERARCHICAL CLUSTERING UPON THE CENTROIDS OF CROSSING THE 2 KMEANS PARTITIONS

d2 <- dist(cdclas)
h2 <- hclust(d2,method="ward",members=freq)  # COMPARE THE COST

plot(h2)
barplot(h2$height[(dim(cdclas)[1]-40):(dim(cdclas)[1]-1)])

c2 <- cutree(h2,nc)

# WARNING c2 NOT ALLOW TO CLASSIFY DIRECTLY THE 4446 INDIVIDUALS BUT THE ELEMENTS OF freq
# WARNING cdclas CONTAINS THE COORDINATES OF THE CROSSINGS BUT THEY NEED TO BE WEIGHED BY  freq

#cdg <- aggregate(cdclas,list(c2),mean)[,2:(nd+1)]  ALL ELEMENTS COUNT EQUAL

cdg <- aggregate((diag(freq/sum(freq)) %*% as.matrix(cdclas)),list(c2),sum)[,2:(nd+1)]
cdg


# CONSOLIDATION

k6 <- kmeans(Psi,centers=cdg)
k6$size

Bss <- sum(rowSums(k6$centers^2)*k6$size)
Wss <- sum(k6$withinss)

Ib6 <- 100*Bss/(Bss+Wss)
Ib6


# PROBALISTIC CLUSTERING 

library(mclust)

emc <- Mclust(Psi,G=7:9) # AGAIN THE COST

print(emc)

attributes(emc)

# NOW FOR EACH INDIVIDUAL WE HAVE A PROBABILITY TO BELONG TO EACH CLASS

emc$z[1:10,]

# LETS SEE TO WHICH CLASS IS ASSGINED EVERY INDIVIDUAL (BY MAX PROB)

emc$classification[1:10]

# LETS COMPUTE THE QUALITY OF THE PARTITION

cdg <- aggregate(as.data.frame(Psi),list(emc$classification),mean)[,2:(nd+1)]
cdg

Bss <- sum(rowSums(cdg^2)*as.numeric(table(c1)))

Ib7 <- 100*Bss/Tss
Ib7

# LETS SEE THE PARTITION VISUALLY

plot(Psi[,1],Psi[,2],col=emc$classification,main="Clustering of credit data in 8 classes")
legend("topleft",c("c1","c2","c3","c4","c5","c6","c7","c8"),pch=1,col=c(1:8))