############################################################
#
# MBD - Clustering supervisado (I)
#
############################################################

############################################################
#
# KNN
#
############################################################

############################################################
# Variables
############################################################
# Datos de clientes de una companyia telefonica
# age: edad
# annualincome: ingresos anuales
# calldroprate: tasa de llamadas perdidas
# callfailurerate: tasa de llamadas cortadas
# callingnum: numero de telefono
# customerid: identificador del cliente
# customersuspended:
# education: formacion educativa
# gender: sexo
# homeowner: propietario de vivienda
# maritalstatus: estado civil
# monthlybilledamount: factura mensual
# noadditionallines: ?
# numberofcomplaints: numero de quejas
# numberofmonthunpaid: numero de meses sin pagar
# numdayscontractequipmentplanexpiring: numero de dias para expirar el contrato
# occupation: ocupacion
# penaltytoswitch: penalizacion por cambio de linea
# state: estado de USA
# totalminsusedinlastmonth: minutos en ultimo mes
# unpaidbalance: balance de impagos
# usesinternetservice: servicio de internet contratado
# usesvoiceservice: servicio de voz contratado
# percentagecalloutsidenetwork
# totalcallduration: duracion total de las llamadas
# avgcallduration: duracion media de las llamadas
# churn: baja del cliente
# year: anyo
# month: mes

rm(list=ls())

############################################################
# Cargar paquetes
############################################################
# install.packages('class')
# install.packages('deldir')
# install.packages('kknn')
library(deldir)
library(kknn)
library(class)

############################################################
# Leer e inspeccionar los datos
############################################################

d0 <- read.table('edw.csv',header=TRUE,sep=',')
dim(d0)
#20468 clientes y 29 variables
summary(d0)
#Mirar categoricas y numericas. La educacion hay 2 opciones, crear 3 variables dummy que valdria 0, 1, 2 o 3
#en funcion de si son de bachelor or equivalent, etc. Haciendo esto estamos valorando de que de un estudio a otro
#hay la misma diferencia de dificultad.

############################################################
# Eliminar variables irrelevantes y transformar algunas
############################################################
var.rm <- c('customerid','noadditionallines','state','year','month','callingnum')
d <- d0[,-which(names(d0) %in% var.rm)]
#eliminamos variables como año (siempre 2015), numero de telefono ya que no da info, etc.

d$education <- relevel(d$education,'High School or below')
#pone la categoria de referencia la mas bajas de todas, de esta forma al transformarla en niveles te cogera
#high school como la mas pequeña. En este caso solo estaba high school desordenada. Si no fuese asi habria que
#indicarle el orden con factor(d$education, levels c('High school,...,))

d$occupationTech <- ifelse(d$occupation=='Technology Related Job',1,0)
#la tranformamos en 2 dummies. 1 o 0

d$occupationNonTech <- ifelse(d$occupation=='Non-technology Related Job',1,0)
d$occupation <- NULL
d <- d[,names(d)[c(1:21,23,24,22)]] # Poner la variable respuesta la ultima (Cuestion totalmente estetica!!!)

############################################################
# Tranformar todas las variables en numericas
############################################################
for (i in 1:ncol(d)) d[,i] <- as.numeric(d[,i])
#De esta forma las categoricas que se han transformado en dummies se transforman en numericas

############################################################
# Dividir la muestra
############################################################
##-- Dividir en muestra de entrenamiento y muestra test
p <- 0.7                 # Proporcion en muestra de entrenamiento
n <- dim(d)[1]           # numero de observaciones
set.seed(12345)
train.sel <- sample(c(FALSE,TRUE),n,rep=TRUE,prob=c(1-p,p))
#cramos un vector donde true cuando te envia a la muestra train i false test
train <- d[train.sel,]
test <- d[!train.sel,]

############################################################
# Comparar capacidad predictiva
############################################################
##-- 1-NN Train + Test
knn1 <- knn(train[,-ncol(d)], test=test[,-ncol(d)], cl=train$churn, k = 1)
#le quitamos la variable respuesta y más cosas (nose!!! :( )
t <- table(knn1,test$churn)
#0 y 1 lo que ha predicho el modelo. No darse de baja =1. nos interesa la Proporcion de aciertos la diagonal decreciente
sum(diag(t))/sum(t)
#prob de acierto de 84%. Si el objetivo fuera predecir los que se dan de baja (71) no los esta cazando bien.

##-- 1-NN Cross-validation
knn2 <- knn.cv(d[,-ncol(d)], cl=d$churn, k = 1)
t <- table(knn2,d$churn)
t
#proporcion de aciertos
sum(diag(t))/sum(t)
#se ha mejorado un poco pero no mucho

##-- Naive (Asignar a la categoria mayoritaria)
max(prop.table(table(test$churn)))
#asignandolos todos al 0 tenemos 90% de acierto pero lo que pasa es que nos interesa los del 1

############################################################
# Numero de grupos
############################################################
p <- c()
#prob de acierto
K <- seq(1,21,2)
#vector que nos diga que vayamos probando los vecinos proximos impares.

for(k in K){
  cat('Iteration:',k,'of',max(K),'\n')
  knn <- knn(train[,-ncol(d)], test=test[,-ncol(d)], cl=train$churn, k = k)
  t <- table(knn,test$churn)
  p[(k+1)/2] <- sum(diag(t))/sum(t)
}
plot(K,p,pch=19,type='b')
#prob de acierto en funcion de k. que no nos engañe la escala!!! de 0.84 a 0.9!
#No aplicamos la regla del codo!! lo que tenemos que coger es el valor mas alto
#a partir de 9 siempre cogemos el valor mas alto
cbind(K,p)
#ver el plot con numeros. en el valor 21 si hacemos t <- table(knn,test$churn) vemos que de los 21 vecinos acabmos
#cogiendolos todos y los asignan al 0.
#si nos quedamos con el 5 ganamos un poco respecto al 0.84 que teniamos antes. Viendo el t no se xd


############################################################
# Usar ACP --> Reduccion dimensionalidad
############################################################
#intentar proyector todos los puntos en aquelespacio que recoja la mayor variabilidad posible
res.acp0 <- princomp(d[,-ncol(d)])
screeplot(res.acp0,type='lines')
#Nos dice el % de variabilidad que recoge cada una de las dimensiones. Nos dice que podemos proyector todos los datos
#en una sola dimension. Por tanto reducimos todos los datos en una dimension.
res.acp <- res.acp0$scores[,1]
#nos quedamos con una componente ( si hubieran salido 3 pos habriamos cogido 3)
#recogemos las variables y las llevamos a una unica variable
train2 <- data.frame(c1=res.acp[train.sel])
test2 <- data.frame(c1=res.acp[!train.sel])


K <- seq(1,21,2)
p <- c()
for(k in K){
  cat('Iteration:',k,'of',max(K),'\n')
  knn <- knn(train2, test2, cl=train$churn, k = k)
  t <- table(knn,test$churn)
  p <- c(p,sum(diag(t))/sum(t))
}
#ahora va mas rapido y hemos conseguido una capacidad predictiva de 0.95. puede parecer poco predictivo coger menos datos

plot(K,p,pch=19,type='b')
cbind(K,p)

##-- 1-NN
knn1 <- knn(train2, test2, cl=train$churn, k = 1)
t <- table(knn1,test$churn)
t
#si vemos la tabla la proporcion de los que estan asignados a que realmente se vayan a dar de baja es mucho mayor.
# si llegase un cliente futuro esperariamos ver estas predicciones.

############################################################
# Anyadir un minimo de votos (parametro l)
############################################################
knn2 <- knn(train2, test2, cl=train$churn, l=2, k = 2)
#cogemos numero par (2 vecinos mas proximos) exijo que o los dos vecinos son iguales o no se asigna a ningun sitio
t <- table(knn2,test$churn)
t
#funciona un poco mejor que la anterior pero hacemos la pequeña trampa de que no esta toda la muestra test por
#los que han empatado.
#los NAs son los que habian empate.
sum(diag(t))/sum(t)
#el % no es real porque es sobre los datos que hemos asignado y no sobre el total
tmiss <- table(knn2,test$churn,useNA = 'always')
tmiss
#dice que 444 no han sido asignados.


############################################################
# Kernel
############################################################
##-- Con datos originales
kknn1 <- kknn(factor(churn)~., train, test,k=1)
#para un vecino proximo
fit <- fitted(kknn1)
t <- table(fit,test$churn)
sum(diag(t))/sum(t)
#con todas estas mejoras hemos llegado del 84 al 95%

p <- c()
K <- seq(1,21,2)
for(k in K){
  cat('Iteration:',k,'of',max(K),'\n')
  kknn <- kknn(factor(churn)~., train, test,k=k)
  t <- table(fitted(kknn),test$churn)
  p <- c(p,sum(diag(t))/sum(t))
}
plot(K,p,pch=19,type='b')
cbind(K,p)

##-- Con ACP --> NO mejora
kknn3 <- kknn(factor(train$churn)~., train2, test2,k=2)
fit <- fitted(kknn3)
t <- table(fit,test$churn)
sum(diag(t))/sum(t)

K <- seq(1,21,2)
for(k in K){
  cat('Iteration:',k,'of',max(K),'\n')
  kknn <- kknn(factor(train$churn)~., train2, test2,k=k)
  t <- table(fitted(kknn),test$churn)
  p[(k+1)/2] <- sum(diag(t))/sum(t)
}
plot(K,p,pch=19,type='b')
cbind(K,p)

############################################################
# Ejercicio 6-1
############################################################

##-- 1. Lee los datos "bank-additional-full.csv"

datos <- read.table('bank-additional-full.csv',header=TRUE,sep=';')

##-- 2. Realiza una depuracion
# a. Ordena las categoricas ordinales
# b. Convierte en dummies las categoricas nominales
# c. Todas las variables deben ser de clase "numeric"
# d. Elimina varible pdays
# Puedes usar la siguiente funcion manual para pasar una variable a dummy

datos$pdays <- NULL
to.dummy <- function(x,ordinal,ord=NULL){
  if(ordinal){
    y <- as.numeric(factor(x,levels=ord))
  }else{
    nl <- nlevels(x)
    y <- matrix(ncol=nl-1,nrow=length(x))
    for(i in 1:(nl-1)){
      y[,i] <- ifelse(x==levels(x)[i],1,0)
    }
    colnames(y) <- paste0(levels(x)[1:(nl-1)],'.dummy')
  }
  return(y)
}
job2 <- to.dummy(datos$job,FALSE)
marital2 <- to.dummy(datos$marital,FALSE)
education2 <- to.dummy(datos$education,TRUE,
                       c('unknown','illiterate','basic.4y','basic.6y','basic.9y',
                         'professional.course','high.school','university.degree'))
default2 <- to.dummy(datos$default,FALSE)
housing2 <- to.dummy(datos$housing,TRUE,c('unknown','no','yes'))
loan2 <- to.dummy(datos$loan,TRUE,c('unknown','no','yes'))
month2 <- to.dummy(datos$month,FALSE)
day_of_week2 <- to.dummy(datos$day_of_week,FALSE)
#poutcome2 <- to.dummy(datos$poutcome,FALSE)

datos$job <- NULL
datos$marital <- NULL
datos$default <- NULL
datos$housing <- NULL
datos$loan <- NULL
datos$month <- NULL
datos$day_of_week <- NULL
datos$poutcome <- NULL

datos <- cbind(datos,job2,marital2,education2,default2,housing2,loan2,month2,day_of_week2)#,poutcome2)
dim(datos)
for(i in 1:length(datos)) datos[,i] <- as.numeric(datos[,i])

##-- 3. Divide la muestra: Train + test en una submuestra de 10000

k <- seq(1,21,2)
p <- c()
for (k in k) {
  cat('Iteration', k, 'of', max(k), '\n')
  knn1 <- knn(train, test, cl=y.train, k = k)
  t <- table(knn1, y.test)
  p[(k+1)/2] <- sum(diag(t))/sum(t)

}

cbind(k,p)

##-- 4. Estudia el numero de grupos idoneo en una submuestra de 10000




############################################################
#
# Naive Bayes
#
############################################################
rm(list=ls())

############################################################
# Variables
############################################################
# Objetivo: clasificar los productos
# id: identificador del producto
# feat_XX: caracter�stica XX
# target: tipo de producto

############################################################
# Cargar paquetes
############################################################
# install.packages('e1071')
# install.packages('ineq')
library(e1071)
library(ineq)

############################################################
# Leer datos
############################################################
setwd('C:/Users/jcortes/Google Drive/Docencia/Master/17_18/Datasets')
datos <- read.csv2('products.csv',header=TRUE)

############################################################
# Inspeccionar datos
############################################################
dim(datos)                                 # Dimension
summary(datos)                             # Descriptiva
#mucha variedad
table(apply(apply(datos,2,is.na),2,sum))   # Tabla con numero de missings
#construye un dataframe con true o flase dependiendo de NAs. vemos que todas las columnas tienen 0 por lo tanto no hay missings NAs

############################################################
# Premisa de independencia
############################################################
cor.matrix <- cor(datos[,-c(1,length(datos))])              # Matriz de correlaciones
cor.num <- as.numeric(cor.matrix)                           # Todas las correlaciones en un vector
t.cor <- table(cut(cor.num[cor.num!=1],br=seq(-1,1,0.1)))/2 # Categorazi�n de las correlaciones en intervalos de 0.1
#tabla categorizada
t.cor
barplot(t.cor, las=2)
#manera de encontrar si hay correlaciones altas o no. Hay diferentes manera de hacerlos.
#son variables bastante independientes. La mayoria estan por el 0 y por tanto poco correlacionados.

############################################################
# Dividir la muestra
############################################################
d <- datos[,-1]                            # Eliminamos identificador
p <- 0.7                                   # Proporcion en muestra de entrenamiento
n <- nrow(d)                               # Numero de observaciones
set.seed(12345)
train.sel <- sample(c(FALSE,TRUE),n,rep=TRUE,prob=c(1-p,p))
train <- d[train.sel,]
test <- d[!train.sel,]

############################################################
# Aplicar bayes
############################################################
nb <- naiveBayes(target ~ ., train)
nb$levels       # Clases
nb$apriori      # A priori
#pertenencia a cada clase
nb$tables       # A posteriori
#probabilidades a posteriori
#para la carac[,1] la media de la desviacion tipo¿? no se si cierto...
#nos gustaria que las medias[,1] fueran diferentes

############################################################
# Capacidad predictiva
############################################################
##-- Global
preds <- predict(nb, newdata = test)
t <- table(preds, test$target)
t
#predicciones respecto a la realidad. queremos valores en la diagonal decreciente solo!!! no pasa del todo
#la clase 5 la predice bastante bien
p.acierto <- sum(diag(t))/sum(t)
p.acierto
#53% al ser 9 clases no esta mal ya que lo normal seria un 1/9

##-- Se pueden pedir las probabilidades de cada clase para inspecci�n visual
preds2 <- predict(nb, newdata = test,type = "raw")
#dame todas las probabilidades tete
#lo que me da es la probabilidad de que el individuo y vaya a la clase x
head(preds2)
heatmap(scale(preds2[1:50,]),Rowv=NA,Colv=NA,col = cm.colors(256))
#ver si en la asignacion se tiene dudas o no.
#en el indivio 1 vemos que hay una predominancia en la clase 2. vemos que en cada linia hay solo un lila por tanto es bien

##-- Proporci�n de acierto por clase
barplot(diag(prop.table(t,2)))
#la prob por columnas de la tabla de confusion. en la clase 6 un 73% van al...
#la clase 5 es la que mas acierta y la 3 menos.


############################################################
# Importancia de las variables
############################################################
# Medir la importancia de las variables en funcion del indice Gini
graphics.off()                  # Cerrar ventanas
nbt <- nb$tables                # Tablas con medias
var.imp <- c()                  # Gini

##-- Calcular indice y graficar distribucion
for (i in 1:(ncol(d)-1)){
  ##-- Crear ventana grafica cada 9 graficos
  if((i-1) %% 9==0){
    windows()
    par(mfrow=c(3,3))
  }

  x <- nbt[[i]][,1]
  barplot(x,main=names(d)[i])  # grafico
  var.imp[i] <- ineq(x)        # GINI
  #indice gini es la medida de la desigualdad entre las clases
  #las variables que tenian valores muy parecidos gini mu bajo
}
#La feat 93 tiene pinta de ser buena predictora ( es bien!)
#la 65 es mal :(

##-- Los mas predictores
sel.pr <- order(var.imp, decreasing=TRUE)[1:9]
#9 variables que mas discriminan
windows()
par(mfrow=c(3,3))
for (i in sel.pr){
  x <- nbt[[i]][,1]
  barplot(x,main=names(d)[i])
}


##-- Los menos predictores
sel.pr <- order(var.imp, decreasing=FALSE)[1:9]
#9 variables que menos discriminan
windows()
par(mfrow=c(3,3))
for (i in sel.pr){
  x <- nbt[[i]][,1]
  barplot(x,main=names(d)[i])
}

############################################################
# Intento de mejora 1: Seleccionar mejores predictores
############################################################
##-- Selecionar los mejores predictores
hist(var.imp)
sel.pr <- which(var.imp>0.5)
#todos los que tengan un indice mayor que 0.5 (menos predictores)
train1 <- train[,c(paste0('feat_',sel.pr),'target')]

##-- Aplicar bayes nuevamente
nb1 <- naiveBayes(target ~ ., train1, type="class")
preds <- predict(nb1, newdata = test)
t <- table(preds, test$target)
p.acierto1 <- sum(diag(t))/sum(t)
p.acierto1

############################################################
# Intento de mejora 2: Quitar variables correlacionadas
############################################################
##-- Sistema 1: Correlaciones grandes
high.corr <- which(cor.matrix>0.8,arr.ind = TRUE)
#quitamos las variables con correlacion mayor de 0.8
t.high.corr <- sort(table(as.numeric(high.corr))/2,decreasing=TRUE)
t.high.corr
sel.rm <- names(t.high.corr)[t.high.corr>=2]
train2 <- train[,-which(names(train) %in% paste0('feat_',sel.rm))]
#escoger tu el criterio. Un criterio posible para cada variable sumar todas las correlaciones con todas las otras
#como todas son numericas puedes hacer una regresion lineal y las que tienen un R2 mas grande ir eliminandolas.


##-- Aplicar bayes nuevamente
nb2 <- naiveBayes(target ~ ., train2, type="class")
preds <- predict(nb2, newdata = test)
t <- table(preds, test$target)
p.acierto2 <- sum(diag(t))/sum(t)
p.acierto2

##-- Sistema 2: Suma de correlaciones
cor.sum <- sort(apply(cor.matrix,2,sum),decreasing=TRUE)
cor.sum
sel.rm <- names(cor.sum)[1:30]
train3 <- train[,-which(names(train) %in% sel.rm)]

##-- Aplicar bayes nuevamente
nb3 <- naiveBayes(target ~ ., train3, type="class")
preds <- predict(nb3, newdata = test)
t <- table(preds, test$target)
p.acierto3 <- sum(diag(t))/sum(t)
p.acierto3
#se ha subido un 0.8 % alomejor porque ya no se puede mas.

############################################################
# Ejercicio 6-2
############################################################
##-- 1. Reduce el numero de variables predictoras a 20 y aplica el KNN con diferente n�mero de vecinos
## �Mejora la capacidad predictiva medida por la proporci�n de aciertos?

##-- 2. Obten una medida probab�listica del KNN y promediala con la obtenida con el NB. Asigna a un grupo segun esta nueva probabilidad
## �Mejora la capacidad predictiva medida por la proporci�n de aciertos? (15')