MODELISATION PREDICTIVE ET APPRENTISSAGE STATISTIQUE AVEC R

/* ==================================================================== */
/* STEPHANE TUFFERY
/* MODELISATION PREDICTIVE ET APPRENTISSAGE STATISTIQUE  AVEC R
/* ==================================================================== */



# ---------------------------------------------------------------------------------------------------------
# Installation des packages
# ---------------------------------------------------------------------------------------------------------

install.packages('ada')
install.packages('ade4')
install.packages('arules')
install.packages('arulesViz')
install.packages('boot')
install.packages('C50')
install.packages('car')
install.packages('caret')
install.packages('CHAID')
install.packages('combinat')
install.packages('corrplot')
install.packages('doSNOW')
install.packages('e1071')
install.packages('extraTrees')
install.packages('FactoMineR')
install.packages('foreach')
install.packages('foreign')
install.packages('gbm')
install.packages('ggplot2')
install.packages('glmnet')
install.packages('gmodels')
install.packages('grplasso')
install.packages('ipred')
install.packages('kernlab')
install.packages('lattice')
install.packages('leaps')
install.packages('LiblineaR')
install.packages('MASS')
install.packages('missForest')
install.packages('nnet')
install.packages('plsRglm')
install.packages('prim')
install.packages('pROC')
install.packages('questionr')
install.packages('randomForest')
install.packages('randtoolbox')
install.packages('rgl')
install.packages('rgrs')
install.packages('ROCR')
install.packages('rpart')
install.packages('rpart.plot')
install.packages('sas7bdat')
install.packages('snow')
install.packages('speedglm')
install.packages('tree')



# ---------------------------------------------------------------------------------------------------------
# CHAPITRE 2
# Lecture du fichier texte
# ---------------------------------------------------------------------------------------------------------

# noms des variables
myVariableNames<-c("Comptes","Duree_credit","Historique_credit","Objet_credit",
"Montant_credit","Epargne","Anciennete_emploi","Taux_effort",
"Situation_familiale","Garanties","Anciennete_domicile","Biens","Age",
"Autres_credits","Statut_domicile","Nb_credits","Type_emploi",
"Nb_pers_charge","Telephone","Etranger","Cible")

# nom des variables en anglais
myVariableNamesE<-c("checking_status","duration","credit_history",
"purpose","credit_amount","savings","employment","installment_rate",
"personal_status","other_parties","residence_since","property_magnitude",
"age","other_payment_plans","housing","existing_credits","job",
"num_dependents","telephone","foreign_worker","class")

# lecture du fichier texte sur Internet
credit = read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data",h=FALSE,col.names=myVariableNames)

# ajout d'un identifiant (le no de ligne)
credit$Cle <- seq(1,nrow(credit))

# recodage de certaines variables
credit$Comptes <- as.factor(substr(credit$Comptes,3,3))
credit$Anciennete_emploi <- substr(credit$Anciennete_emploi,3,3)
credit$Epargne <- substr(credit$Epargne,3,3)
credit$Epargne[credit$Epargne == "5" ] <- "0"
credit$Etranger <- NULL



# ---------------------------------------------------------------------------------------------------------
# Passage de SAS à R avec le package sas7bdat qui sait lire une table SAS non compressée
# ---------------------------------------------------------------------------------------------------------

install.packages("sas7bdat")
library(sas7bdat)
credit = read.sas7bdat("C:\\Users\\TUFFERY\\Documents\\...\\Etude de cas\\credit.sas7bdat")



# ---------------------------------------------------------------------------------------------------------
# Passage de SPSS à R avec le package foreign
# ---------------------------------------------------------------------------------------------------------

install.packages("foreign")
library(foreign)
credit = read.spss("C:\\Users\\TUFFERY\\Documents\\...\\Etude de cas\\credit.sav",use.value.labels=TRUE, max.value.labels=Inf, to.data.frame=TRUE)



# ---------------------------------------------------------------------------------------------------------
# Préparation des données
# ---------------------------------------------------------------------------------------------------------

credit$Cible[credit$Cible == 1] <- 0 # crédits OK
credit$Cible[credit$Cible == 2] <- 1 # crédits KO
#credit$Cible <- credit$Cible - 1
credit$Cible <- factor(credit$Cible)

# transformation de l'ensemble des variables qualitatives en facteurs
varquali <- c("Comptes", "Epargne", "Historique_credit", "Objet_credit",
 "Situation_familiale", "Garanties", "Biens", "Autres_credits",
 "Statut_domicile", "Type_emploi", "Anciennete_emploi", "Telephone", "Nb_pers_charge")
indices <- names(credit) %in% varquali
for (i in (1:ncol(credit))) { if (indices[i] == 1) { credit[,i] <- factor(credit[,i]) } }

# liste des variables quantitatives
varquanti <- c("Duree_credit", "Montant_credit", "Taux_effort",
"Anciennete_domicile", "Age", "Nb_credits")

# exclusion de la clé de la liste des variables analysées
names(credit)
vars <- -grep('Cle', names(credit))



# ---------------------------------------------------------------------------------------------------------
# Constitution d'un échantillon de test
# ---------------------------------------------------------------------------------------------------------

# récupération des clés de l'échantillon de test
cles = read.table("http://blogperso.univ-rennes1.fr/stephane.tuffery/public/cles.txt",h=FALSE,col.names="Cle")
id <- cles$Cle

# variante : tirage comme RANUNI de SAS avec graine = 123
library(randtoolbox)
a <- 397204094
b <- 0
m <- 2^(31)-1
set.generator(name="congruRand", mod=m, mult=a, incr=b, seed=123)
id <- (1:1000)[which(runif(1000) < 0.66)]
head(id,20)
all.equal(id,cles$Cle)
set.generator("default")

# échantillon de validation
test  <- credit[-id,]



# ---------------------------------------------------------------------------------------------------------
# CHAPITRE 3
# Premières analyses descriptives
# ---------------------------------------------------------------------------------------------------------

# statistiques de base
summary(credit)

# histogrammes de toutes les variables
library(Hmisc)
hist.data.frame(credit[,1:16])

# histogramme
library(lattice)
histogram(~Duree_credit | Cible , data = credit, type="percent", col="grey", breaks=10)
histogram(~Montant_credit | Cible , data = credit, type="percent", col="grey", breaks=10)
histogram(~Age | Cible , data = credit, type="percent", col="grey", breaks=10)

# statistiques descriptives par groupes
aggregate(credit[,c("Age","Duree_credit")],by=list(Cible=credit$Cible),mean)
aggregate(credit[,c("Age","Duree_credit")],by=list(Cible=credit$Cible),summary)
by(credit[,c("Age","Duree_credit","Montant_credit")],list(Cible=credit$Cible),summary)
by(credit[,1:20],list(Cible=credit$Cible),summary)
tapply(credit$Age,credit$Cible,summary)
#chi² de Kruskal-Wallis
kruskal.test(credit$Age~credit$Cible)$statistic
kruskal.test(credit$Duree_credit~credit$Cible)$statistic
kruskal.test(credit$Montant_credit~credit$Cible)$statistic
#chi² de Kruskal-Wallis normalisé
sqrt(kruskal.test(credit$Age~credit$Cible)$statistic/sum(!is.na(credit$Age)))
sqrt(kruskal.test(credit$Duree_credit~credit$Cible)$statistic/sum(!is.na(credit$Duree_credit)))
sqrt(kruskal.test(credit$Montant_credit~credit$Cible)$statistic/sum(!is.na(credit$Montant_credit)))

# superposition des fonctions de densité des bons et mauvais dossiers
plot(density(credit$Age[credit$Cible==0]),main="Fonction de densité",col="blue",
xlim = c(18,80), ylim = c(0,0.1),lwd=2)
par(new = TRUE)
plot(density(credit$Age[credit$Cible==1]),col="red",lty=3,lwd=2,
xlim = c(18,80), ylim = c(0,0.1),xlab = '', ylab = '',main=' ')
legend("topright",c("Cible=0","Cible=1"),lty=c(1,3),col=c("blue","red"),lwd=2)

# taux d'impayés par décile de l'âge
q <- quantile(credit$Age, seq(0, 1, by=0.1))
q[1] <- q[1]-1 # comme les intervalles sont en (a,b] et que q[1]=19
# et que la base contient deux clients de 19 ans
# il faut que la 1ère borne vaille 18 pour ne pas exclure ces deux clients
qage <- cut(credit$Age,q)
tab <- table(qage,credit$Cible)
prop.table(tab,1) # affichage % en ligne
barplot(prop.table(tab,1)[,2],las=3,main="Age",ylab="Taux impayés",density=0)
abline(h=.3,lty=2)

# fonction d'affichage du taux d'impayés par quantile
ti = function(x,pas=0.1)
{
q <- unique(quantile(x, seq(0, 1, by=pas)))
unique(q)
q[1] <- q[1]-1 # comme les intervalles sont en (a,b]
qx <- cut(x,q)
tab <- table(qx,credit$Cible)
print(prop.table(tab,1)) # affichage % en ligne
barplot(prop.table(tab,1)[,2],las=3,main=deparse(substitute(x)),ylab="Taux impayés",density=0,horiz=F)
abline(h=prop.table(table(credit$Cible))[2],lty=2)
}
	 
ti = function(x,pas=0.1)
{
q <- unique(quantile(x, seq(0, 1, by=pas)))
unique(q)
q[1] <- q[1]-1 # comme les intervalles sont en (a,b]
qx <- cut(x,q)
tab <- table(qx,credit$Cible)
print(prop.table(tab,1)) # affichage % en ligne
par(mar = c(4, 7, 4, 2))
barplot(prop.table(tab,1)[,2],las=1,main=deparse(substitute(x)),ylab="Taux impayés \n\n",density=0,horiz=T)
abline(v=prop.table(table(credit$Cible))[2],lty=2)
}

ti(credit$Age)
ti(credit$Duree_credit,pas=0.05)
ti(credit$Montant_credit,pas=0.05)

# discrétisation des variables continues
Age <- cut(credit$Age,c(0,25,Inf))
tab <- table(Age,credit$Cible)
prop.table(tab,1) # affichage % en ligne
Duree_credit <- cut(credit$Duree_credit,c(0,15,36,Inf))
tab <- table(Duree_credit,credit$Cible)
prop.table(tab,1) # affichage % en ligne
Montant_credit <- cut(credit$Montant_credit,c(0,4000,Inf))
tab <- table(Montant_credit,credit$Cible)
prop.table(tab,1) # affichage % en ligne

# tableaux croisés - version 1
ct <- function(x)
{ print(names(credit)[x])
prop.table(table(credit[,x],credit$Cible),1) }
#for (i in (1:ncol(credit))) print(ct(i))
for (i in (1:ncol(credit)))
{ if ((names(credit[i])) %in% c("Cle","Cible","Duree_credit","Montant_credit","Age"))
	{} else { print(ct(i)) }
}

# tableaux croisés - version 2 un peu plus jolie
ct <- function(x)
{ cat("\n",names(credit)[x])
prop.table(table(credit[,x],credit$Cible),1) }
for (i in (1:ncol(credit)))
{ if (!(names(credit[i])) %in% c("Cle","Cible","Duree_credit","Montant_credit","Age"))
	{ print(ct(i)) }
}

# tableaux croisés - version 3 avec effectifs
ct <- function(x)
{ cat("\n",names(credit)[x],"\n")
cbind(prop.table(table(credit[,x],credit$Cible),1),
addmargins(table(credit[,x]))[-(nlevels(credit[,x])+1)]) }
for (i in (1:ncol(credit)))
{ if (!(names(credit[i])) %in% c("Cle","Cible","Duree_credit","Montant_credit","Age"))
	{ print(ct(i)) }
}

# tableaux croisés comme SAS et SPSS
install.packages('gmodels')
library(gmodels)
CrossTable(credit$Comptes,credit$Cible,prop.chisq=F,chisq = T,format="SAS")

# chi² et V de Cramer des variables explicatives
# les 3 variables continues étant sous leur forme discrétisées
npred <- -grep('(Cle|Cible|Duree_credit|Montant_credit|Age)', names(credit))
credit2 <- credit[,npred]
credit2$Age <- Age
credit2$Duree_credit <- Duree_credit
credit2$Montant_credit <- Montant_credit

# 1e solution pour le V de Cramer :
# utilisation du package rgrs ou questionr
library(rgrs)
#install.packages("questionr")
library(questionr)
cramer  <- matrix(NA,ncol(credit2),3)
for (i in (1:ncol(credit2)))
{     cramer[i,1] <- names(credit2[i])
	cramer[i,2] <- cramer.v(table(credit2[,i],credit$Cible))
	cramer[i,3] <- chisq.test(table(credit2[,i],credit$Cible))$p.value
}
colnames(cramer) <- c("variable","V de Cramer","p-value chi2")

# 2e solution pour le V de Cramer :
# utilisation de la formule : ici, racine carrée de (chi2 / effectif)
cramer  <- matrix(NA,ncol(credit2),3)
for (i in (1:ncol(credit2)))
{     cramer[i,1] <- names(credit2[i])
	cramer[i,2] <- sqrt(chisq.test(table(credit2[,i],credit$Cible))$statistic/
(length(credit2[,i])))
	cramer[i,3] <- chisq.test(table(credit2[,i],credit$Cible))$p.value
}
colnames(cramer) <- c("variable","V de Cramer","p-value chi2")

# affichage des variables par V de Cramer décroissants
vcramer <- cramer [order(cramer[,2], decreasing=T),]
vcramer
# graphique
par(mar = c(8, 4, 4, 0))
barplot(as.numeric(vcramer[,2]),col=gray(0:nrow(vcramer)/nrow(vcramer)),
names.arg=vcramer[,1], ylab='V de Cramer', ylim=c(0,0.35),cex.names = 0.8, las=3)



# ---------------------------------------------------------------------------------------------------------
# Liaisons entre variables explicatives
# ---------------------------------------------------------------------------------------------------------

# V de Cramer des paires de variables explicatives
install.packages('questionr')
library(questionr)
cramer  <- matrix(NA,ncol(credit2),ncol(credit2))
# variante 1
for (i in (1:ncol(credit2)))
{     for (j in (1:ncol(credit2)))
	{
	cramer[i,j] <- cramer.v(table(credit2[,i],credit2[,j]))
	}
}
# variante 2
vcram <- function(i,j) { cramer.v(table(credit2[,i],credit2[,j])) }
i <- (1:ncol(credit2))
j <- (1:ncol(credit2))
cramer <- outer(i,j,Vectorize(vcram))
# fin variantes
colnames(cramer) <- colnames(credit2)
rownames(cramer) <- colnames(credit2)
cramer
#install.packages("corrplot")
library(corrplot)
corrplot(cramer)
corrplot(cramer, method="shade", shade.col=NA, tl.col="black", tl.srt=45)
par(omi=c(0.4,0.4,0.4,0.4))
corrplot(cramer,type="upper",tl.srt=45,tl.col="black",tl.cex=1,diag=F,addCoef.col="black",addCoefasPercent=T)

# croisement des variables explicatives les plus fortement liées

tab <- table(credit$Biens,credit$Statut_domicile)
prop.table(tab,1) # affichage % en ligne

tab <- table(Duree_credit,Montant_credit)
prop.table(tab,1) # affichage % en ligne
tab <- table(Duree_credit,credit$Cible)
cbind(prop.table(tab,1),addmargins(tab,2)) # affichage % en ligne
tab <- table(Montant_credit,credit$Cible)
cbind(prop.table(tab,1),addmargins(tab,2)) # affichage % en ligne

tab <- table(credit$Type_emploi,credit$Telephone)
prop.table(tab,1) # affichage % en ligne

tab <- table(credit$Biens,credit$Statut_domicile)
prop.table(tab,1) # affichage % en ligne

tab <- table(credit$Nb_credits,credit$Historique_credit)
prop.table(tab,1) # affichage % en ligne



# ---------------------------------------------------------------------------------------------------------
# CHAPITRE 4
# Discrétisation automatique des variables continues
# ---------------------------------------------------------------------------------------------------------

# package pour calcul AUC
library(pROC)

decoup = function(base,x,y,h=0,k=0,pAUC=0,nbmod=3,calcul=1,algo="Nelder-Mead",graphe=0)
{

# renommage des variables
attach(base)
Xt <- x
Yt <- y
detach(base)

# traitement spécifique des valeurs manquantes
X  <- Xt[is.na(Xt)==0]
Y  <- Yt[is.na(Xt)==0]

seuils <- as.vector(quantile(X, seq(0, 1, length=(nbmod+1))))[2:nbmod]

fitauc = function(s)
{
s2 <- c(-Inf,unique(s),Inf)
qX <- cut(X,s2)
tab <- table(qX,Y)
logit <- glm(Y~qX, family=binomial(link = "logit"))
qXn <- predict(logit, newdata=base[is.na(Xt)==0,], type="response")
resultat <- auc(Y,qXn, partial.auc=c(1,pAUC), partial.auc.correct=FALSE) *
(1-sum((table(qX)/length(X))^2))/(1-(1-h)*(sum((table(qX)/length(X))^2))) *
((1-(1-k)*(sum((prop.table(tab,1)[,2])^2)))/(1-sum((prop.table(tab,1)[,2])^2)))
return(-resultat)
}

# application du découpage
Applical = function() {
sf <- unique(c(-Inf,est$par,Inf))
qX <- cut(Xt,sf)
tab <- table(addNA(qX),Yt)
cat("\n","Résultat du découpage :","\n")
cat("\n","Seuils     % Négat. % Posit. # +  #     % #","\n")
print(cbind(prop.table(tab,1)*100,tab[,2],table(addNA(qX)),table(addNA(qX))*100/length(Xt)))
cat("\n","Indicateur de convergence (0 = convergence optimisation)","\n")
cat(est$convergence) # vérifier qu'on obtient 0
cat("\n","AUC (partielle) maximisée :","\n")
cat(-est$value)
cat("\n","Homogénéité des classes (0 <- faible ... forte -> 1) :","\n")
cat(1-sum((prop.table(table(addNA(qX),Yt),1)[,2])^2))
return(qX)
}

# calcul aire sous la courbe ROC
Gini = function(t) {
cat("\n","AUC avant découpage :","\n")
logit <- glm(Y~X, family=binomial(link = "logit"))
g1 <- auc(Y,predict(logit, newdata=base[is.na(Xt)==0,], type="response"))
cat(g1)
cat("\n","AUC après découpage :","\n")
logit <- glm(Yt~t,family=binomial(link = "logit"))
g2 <- auc(Yt,predict(logit, newdata=base, type="response"))
cat(g2)
cat("\n","% Evolution AUC avant/après découpage :","\n")
cat(100*(g2-g1)/g1)
cat("\n")
}

# recherche des seuils optimaux à partir des valeurs initiales précédentes (calcul = 1)
# ou utilisation de bornes prescrites (calcul = 0)
if (calcul == 1) { est = optim(seuils,fitauc,method=algo) }
else { est <- list() ; est$par <- bornes ; est$convergence <- 0 ; est$value <- 0}

cat("\n","---------------------------------------------------------------------------","\n")
cat("\n","Discrétisation de", substitute(x), "en", nbmod, "classes ( algorithme ",algo,") \n")
cat("\n","---------------------------------------------------------------------------","\n")

qX1 <- Applical()
Gini(qX1)

# courbe de densité
if (graphe == 1) {
x11()
layout(matrix(c(1, 2), 2, 1, byrow = TRUE),heights=c(3, 1))
par(mar = c(2, 4, 2, 1))
base0 <- X[Y==0]
base1 <- X[Y==1]
xlim1 <- range(c(base0,base1))
ylim1 <- c(0,max(max(density(base0)$y),max(density(base1)$y)))
plot(density(base0),main=" ",col="blue",ylab=paste("Densité de ",deparse(substitute(x))),
xlim = xlim1, ylim = ylim1 ,lwd=2)
par(new = TRUE)
plot(density(base1),col="red",lty=3,lwd=2,
xlim = xlim1, ylim = ylim1,xlab = '', ylab = '',main=' ')
legend("topright",c(paste(deparse(substitute(y))," = 0"),paste(deparse(substitute(y))," = 1")),
lty=c(1,3),col=c("blue","red"),lwd=2)
abline(v=est$par,lty=2)
texte <- c("Chi² de Kruskal-Wallis = \n\n",round(kruskal.test(X~Y)$statistic,digits=3))
text(xlim1[2]*0.8, ylim1[2]*0.5, texte,cex=0.75)
plot(X~Y,horizontal = TRUE,xlab= deparse(substitute(y)),col=c("blue","red"))
}

# fin de la fonction de discrétisation automatique
}

# suppression des objets existants de l'environnement global (GlobalEnv) pour éviter un conflit dans la fonction "attach"
# quand on lui passe en paramètre l'un de ces objets
rm(Age)
rm(Duree_credit)
rm(Montant_credit)
for (i in (2:5)) decoup(credit,Age,Cible,nbmod=i,h=0,k=0,pAUC=0.8,graphe=1)
for (i in (2:5)) decoup(credit,Duree_credit,Cible,nbmod=i,h=0,k=0,pAUC=0.8,graphe=1)
for (i in (2:5)) decoup(credit,Montant_credit,Cible,nbmod=i,h=0,k=0,pAUC=0.8,graphe=1)

# seuils fixés
bornes <- c(26,47)
decoup(credit,age,cible,nbmod=3,calcul=0,h=0,k=0,pAUC=0.8,graphe=1)



# ---------------------------------------------------------------------------------------------------------
# Discrétisation des variables
# ---------------------------------------------------------------------------------------------------------

# création d'un data frame avec variables continues discrétisées
# au vu des taux d'impayés et effectifs de chaque modalité

# discrétisation des variables continues
npred <- -grep('(Cle|Duree_credit|Montant_credit|Age)', names(credit))
credit2 <- credit[,npred]
credit2$Age <- cut(credit$Age,c(0,25,Inf),right=TRUE)
credit2$Duree_credit <- cut(credit$Duree_credit,c(0,15,36,Inf),right=TRUE)
credit2$Montant_credit <- cut(credit$Montant_credit,c(0,4000,Inf),right=TRUE)

# recodage et fusion des modalités
library(car)

# attention aux simples et doubles quotes dans le "recode"
credit2$Comptes <- recode(credit2$Comptes,"4='Pas de compte';1='CC < 0 euros';2='CC [0-200 euros[';3='CC > 200 euros' ")
table(credit2$Comptes)

credit2$Historique_credit <- recode(credit2$Historique_credit,"'A30'='Impayés passés';
'A31'='Impayé en cours dans autre banque';
c('A32','A33')='Pas de crédits ou en cours sans retard';'A34'='Crédits passés sans retard'")
table(credit2$Historique_credit)

credit2$Objet_credit <- recode(credit2$Objet_credit,"'A40'='Voiture neuve';
'A41'='Voiture occasion';c('A42','A44','A45')='Intérieur';
'A43'='Vidéo - HIFI';c('A46','A48')='Etudes';
'A47'='Vacances';'A49'='Business';else='Autres'")
table(credit2$Objet_credit)

credit2$Epargne <- recode(credit2$Epargne,
"0='Sans épargne';1:2='< 500 euros';3:4='> 500 euros'")
table(credit2$Epargne)

credit2$Anciennete_emploi <- recode(credit2$Anciennete_emploi,
"1:2='Sans emploi ou < 1 an';3='entre 1 et 4 ans';4:5='depuis au moins 4 ans'")
table(credit2$Anciennete_emploi)

credit2$Situation_familiale<- recode(credit2$Situation_familiale,
"'A91'='Homme divorcé/séparé';'A92'='Femme divorcée/séparée/mariée';
c('A93','A94')='Homme célibataire/marié/veuf';'A95'='Femme célibataire'")
table(credit2$Situation_familiale)

credit2$Garanties <- recode(credit2$Garanties,"'A103'='Avec garant';else='Sans garant'")
table(credit2$Garanties)

credit2$Biens <- recode(credit2$Biens,"'A121'='Immobilier';'A124'='Aucun bien';
else='Non immobilier'")
table(credit2$Biens)

credit2$Autres_credits <- recode(credit2$Autres_credits,"'A143'='Aucun crédit extérieur';else='Crédits extérieurs'")
table(credit2$Autres_credits)

credit2$Statut_domicile <- recode(credit2$Statut_domicile,"'A152'='Propriétaire';else='Non Propriétaire'")
table(credit2$Statut_domicile)

# structure du fichier de crédit après transformations
summary(credit2)

# échantillons d'apprentissage et de validation
train  <- credit2[id,]
valid  <- credit2[-id,]



# ---------------------------------------------------------------------------------------------------------
# CHAPITRE 5
# Régression logistique avec sélection pas à pas
# ---------------------------------------------------------------------------------------------------------

# formule avec l'ensemble des prédicteurs
predicteurs <- -grep('(Cle|Cible)', names(credit2))
# formule sans interaction
formule <- as.formula(paste("y ~ ",paste(names(credit2[,predicteurs]),collapse="+")))
# formule avec interactions
formule <- as.formula(paste("y ~ ( ",paste(names(credit2[,predicteurs]),collapse="+"),")^2"))
formule

# sélection de modèle ascendante
logit <- glm(Cible~1,data=train,family=binomial(link = "logit"))
summary(logit)
# recherche maximale
selection <- step(logit,direction="forward",trace=TRUE,k = log(nrow(train)),
scope=list(upper=formule))
selection
summary(selection)
# recherche limitée aux principaux prédicteurs
selection <- step(logit,direction="forward",trace=TRUE,k = 2,
scope=list(upper=~Comptes + Duree_credit + Historique_credit + Epargne + Garanties + Age + Autres_credits 
+ Taux_effort + Montant_credit + Statut_domicile + Anciennete_emploi))

# application du modèle à un jeu de données
train.ascbic <- predict(selection, newdata=train, type="response")
valid.ascbic <- predict(selection, newdata=valid, type="response")

# aire sous la courbe ROC
library(ROCR)
pred <- prediction(train.ascbic,train$Cible,label.ordering=c(0,1))
pred <- prediction(valid.ascbic,valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# sélection sur l'AIC sans l'objet de crédit - avec interactions
selection <- step(logit,direction="forward",trace=TRUE,k = 2,
scope=list(upper=~(Comptes + Duree_credit + Historique_credit + 
    Epargne + Garanties + Age + Autres_credits + Taux_effort + 
    Montant_credit + Statut_domicile + Anciennete_emploi)^2))

# sélection de modèle descendante
logit <- glm(Cible~.,data=train,family=binomial(link = "logit"))
selection <- step(logit,direction="backward",trace=TRUE,k = log(nrow(train)))
selection <- step(logit,direction="backward",trace=TRUE,
scope=list(lower=~Comptes+Historique_credit+Duree_credit+Age+Epargne+Garanties+Autres_credits))
selection



# ---------------------------------------------------------------------------------------------------------
# Régression logistique avec sélection globale
# ---------------------------------------------------------------------------------------------------------

# sélection globale
library(leaps)
train  <- credit2[id,]
valid  <- credit2[-id,]

# méthode de Lawless et Singhal
logit <- glm(Cible~.,data=train,family=binomial(link = "logit"))
summary(logit)
pi <- predict(logit, newdata=train, type="response")
poids <- pi * (1 - pi)
y <- as.numeric(train[,"Cible"])
y[y == 1] <- 0 # crédits OK
y[y == 2] <- 1 # crédits KO
train$Cible_lin = log (pi/(1 - pi)) + ( (y - pi) / (pi*(1 - pi)) )
reg <- lm(Cible_lin~.,data=train[,-17],weights=poids)
summary(reg)
fw <- regsubsets(Cible_lin~.,wt=poids,data=train[,-17],nbest=1,nvmax=40,really.big=T)
fw <- regsubsets(Cible_lin~.,wt=poids,data=train[,-17],nbest=1000,nvmax=16,really.big=T)


# fonction leaps

# variables explicatives mises sous forme matricielle
x <- data.frame(model.matrix( ~ .,data=train[,-which(names(train)=="Cible")]))
y <- as.numeric(train[,"Cible"])
y[y == 1] <- 0 # crédits OK
y[y == 2] <- 1 # crédits KO

selec <- leaps(x,y,method="Cp",nbest=1) # impossible si plus de 31 variables
selec <- leaps(x,y,method="Cp",nbest=1,strictly.compatible=F)
plot(selec$size-1,selec$Cp,xlab="# predicteurs",ylab="Cp") # on soustrait 1 à size à cause de l'intercept
selec$size
selec$Cp
selec$which
# modèle correspondant au Cp minimum
best.model <- selec$which[which((selec$Cp == min(selec$Cp))),]
print(best.model)
colnames(x)[best.model] # variables du meilleur modèle

# avec R² ou R² ajusté
#selec <- leaps(x,train[,"Cible"],method="r2",nbest=1,strictly.compatible=F)
selec <- leaps(x,train[,"Cible"],method="adjr2",nbest=1,strictly.compatible=F)
plot(selec$size-1,selec$adjr2,xlab="# predicteurs",ylab="R2 ajusté")
which((selec$adjr2 == max(selec$adjr2)))

# ajustement du modèle sélectionné
formule <- as.formula(paste("y ~ ",paste(colnames(x)[best.model],collapse="+")))
z <- cbind(x,y)
logit <- glm(formule,data=z,family=binomial(link = "logit"))
summary(logit)

# application du modèle à un autre jeu de données
xt <- data.frame(model.matrix( ~ .,data=valid[,-which(names(valid)=="Cible")]))
prob <- predict(logit, newdata=xt, type="response")

# aire sous la courbe ROC
library(ROCR)
pred <- prediction(prob,valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# recherche de l'aire sous la courbe ROC maximale
selec <- leaps(x,y,method="Cp",nbest=10,strictly.compatible=F)
nmodel <- nrow(selec$which)
aucfw <- matrix(NA,nmodel,3)
models <- matrix(NA,nmodel,ncol(selec$which)+1)
for (i in 1:nmodel)
{
best.model <- selec$which[i,]
formule <- as.formula(paste("y ~ ",paste(colnames(x)[best.model],collapse="+")))
logit <- glm(formule,data=z,family=binomial(link = "logit"))
prob <- predict(logit, newdata=xt, type="response")
pred <- prediction(prob,valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]
aucfw[i,1] <- selec$size[i]-1
aucfw[i,2] <- performance(pred,"auc")@y.values[[1]]
models[i,1:ncol(selec$which)] <- selec$which[i,1:ncol(selec$which)]
models[i,ncol(selec$which)+1] <- performance(pred,"auc")@y.values[[1]]
aucfw[i,3] <- selec$Cp[i]
}
colnames(aucfw) <- c("taille","AUC","Cp")
selglob <- aucfw [order(aucfw[,2], decreasing=T),]
selglob

# recherche des meilleures variables par CART appliqué à l'ensemble des modèles précédents
# avec leur AUC comme variable à expliquer
library(rpart)
colnames(models) <- c(names(x),"AUC")
ensmodels <- as.data.frame(models)
cart <- rpart(AUC ~ . ,data = ensmodels,method="anova",parms=list(split="gini"),control=list(maxdepth=3))
cart
plot(cart,branch=.2, uniform=T, compress=T, margin=.1)
text(cart, fancy=T, use.n=T, pretty=0, all=T, cex=1.6, fwidth = 2)

# recherche des meilleures variables par RF appliquées à l'ensemble des modèles précédents
# avec leur AUC comme variable à expliquer
# et en utilisant la fonction "importance"
library(randomForest)
set.seed(235)
rf <- randomForest(AUC ~ . ,data = ensmodels, importance=TRUE, ntree=500, mtry=6, replace=T, keep.forest=T,nodesize=5)
varImpPlot(rf,cex=1.4)
importance(rf,type=2,scale=F)[order(importance(rf,type=2,scale=F), decreasing=TRUE),]

# fonction regsubsets

#fw <- regsubsets(model.matrix( ~ 0+Comptes+Duree_credit+Historique_credit+Objet_credit+Epargne+Anciennete_emploi+Situation_familiale+Garanties+Biens+Age,data=train),train$Cible)
fw <- regsubsets(Cible~.,data=train,nbest=1,nvmax=36)
fw <- regsubsets(Cible~.,data=train,nbest=1,nvmax=16,really.big=T)

# augmentation de la marge du bas pour éviter la légende tronquée pour les montants de crédit
(newOmi <- par()$omi)    
newOmi[1] <- 3     
par(omi=newOmi)
plot(fw,scale="bic")
plot(fw,scale="r2")

selec <- summary(fw)
apply(selec$which,1,sum) # nb de prédicteurs du modèle (y compris la constante)
plot(apply(selec$which,1,sum)-1, selec$bic, ,xlab="# predicteurs",ylab="BIC")

# meilleur modèle selon critère BIC
best.model <- selec$which[which((selec$bic == min(selec$bic))),]
min(selec$bic)
print(best.model)
x <- data.frame(model.matrix( ~ .,data=train[,-which(names(train)=="Cible")]))
colnames(x)[best.model] # variables du meilleur modèle

# recherche de l'aire sous la courbe ROC maximale
#library(ROCR)
fw <- regsubsets(Cible~.,data=train,nbest=1000,nvmax=16,really.big=T)
selec <- summary(fw)
x <- data.frame(model.matrix( ~ .,data=train[,-17]))
y <- as.numeric(train[,"Cible"])
y[y == 1] <- 0 # crédits OK
y[y == 2] <- 1 # crédits KO
z <- cbind(x,y)
xt <- data.frame(model.matrix( ~ .,data=valid[,-17]))
nmodel <- nrow(selec$which)
aucfw <- matrix(NA,nmodel,3)
models <- matrix(NA,nmodel,ncol(selec$which)+1)
for (i in 1:nmodel)
{ if (apply(selec$which,1,sum)[i] > 0)
{
best.model <- selec$which[i,]
formule <- as.formula(paste("y ~ ",paste(colnames(x)[best.model],collapse="+")))
logit <- glm(formule,data=z,family=binomial(link = "logit"))
prob <- predict(logit, newdata=xt, type="response")
pred <- prediction(prob,valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]
aucfw[i,1] <- apply(selec$which,1,sum)[i]-1
aucfw[i,2] <- performance(pred,"auc")@y.values[[1]]
models[i,1:ncol(selec$which)] <- selec$which[i,1:ncol(selec$which)]
models[i,ncol(selec$which)+1] <- performance(pred,"auc")@y.values[[1]]
aucfw[i,3] <- selec$bic[i]
}
}
colnames(aucfw) <- c("taille","AUC","BIC")
selglob <- aucfw [order(aucfw[,2], na.last = NA, decreasing=T),]
head(selglob,20)
hist(selglob[,2],breaks=100,xlim=c(0.6,0.8))



# ---------------------------------------------------------------------------------------------------------
# Régression logistique avec combinaisons de variables
# ---------------------------------------------------------------------------------------------------------

library(combinat)
library(ROCR)

# version for-loop
# -------------------------------------------------------

CombiReg = function(apprent, validat, varY, varX, p)
{
y <- apprent[,varY] # variable à expliquer
cible <- validat[,varY] # variable à prédire
size <- length(varX) # nb total de variables
combi <- combn(size, p) # combinaison de variables
perf  <- matrix(0,dim(combi)[2],1) # dim(combi)[2] = nb de combinaisons
predi <- matrix(" ",dim(combi)[2],p)

for(i in (1:dim(combi)[2]))
{

# sélection des prédicteurs
s <- combi[,i] # ie combinaison de variables
predicteurs <- varX[s]

# écriture de la formule avec les prédicteurs sélectionnés
if (p > 1) {
formule <- as.formula(paste("y ~ ",paste(names(apprent[,predicteurs]),collapse="+")))
} else {
formule <- as.formula(paste("y ~ ",predicteurs))
}
# ajustment du modèle logit
logit <- glm(formule, data=apprent, family=binomial(link = "logit"))
# application du modèle logit à l'échantillon de validation
scores <- predict(logit, newdata=validat, type="response")
pred <- prediction(scores,cible,label.ordering=c(0,1))
#perf[i] <- (performance(pred,"auc")@y.values[[1]]-0.5)*2
perf[i] <- performance(pred,"auc")@y.values[[1]]
predi[i,] <- predicteurs

} # fin de la boucle

cr <- list(perf,predi)

} # fin de la fonction

varx <- c(varquali, varquanti)
cr <- CombiReg(apprent=train, validat=valid, varY="Cible", varX=varx, p=7)
# affichage des combinaisons par AUC décroissantes
resultat <- cbind(cr[[1]],cr[[2]])[order(-cr[[1]]),]
# affichage des 100 meilleures combinaisons
resultat[1:min(100,dim(resultat)[1]),]


# version vectorisée
# -------------------------------------------------------

CombiRegR = function(apprent, validat, varY, varX, p)
{
y <- apprent[,varY] # variable à expliquer
cible <- validat[,varY] # variable à prédire
size <- length(varX) # nb total de variables
combi <- combn(size, p) # combinaison de variables
predi <- matrix(" ",dim(combi)[2],p)

f = function(i)
{
# sélection des prédicteurs
s <- combi[,i] # ie combinaison de variables
predicteurs <- varX[s]

# écriture de la formule avec les prédicteurs sélectionnés
if (p > 1) {
formule <- as.formula(paste("y ~ ",paste(names(apprent[,predicteurs]),collapse="+")))
} else {
formule <- as.formula(paste("y ~ ",predicteurs))
}
# ajustment du modèle logit
logit <- glm(formule, data=apprent, family=binomial(link = "logit"))
# application du modèle logit à l'échantillon de validation
scores <- predict(logit, newdata=validat, type="response")
pred <- prediction(scores,cible,label.ordering=c(0,1))
return(list(auc = performance(pred,"auc")@y.values[[1]],predi = predicteurs))

} # fin de la fonction à vectoriser

cr <- matrix(unlist(Vectorize(f)(seq(1,dim(combi)[2]))),ncol=dim(combi)[2],byrow=F)

} # fin de la fonction

cr <- CombiRegR(apprent=train, validat=valid, varY="Cible", varX=varx, p=7)
# affichage des combinaisons par AUC décroissantes
resultat <- data.frame(t(cr))
colnames(resultat) <- c('AUC',paste("V",1:(ncol(resultat)-1),sep=""))
resultat$AUC <- as.numeric(as.character(resultat$AUC))
resultat <- resultat[order(-resultat$AUC),]
head(resultat,100)



# ---------------------------------------------------------------------------------------------------------
# Régression logistique avec les modalités définitives
# ---------------------------------------------------------------------------------------------------------

# création d'un data frame avec variables continues discrétisées
# pour régression logistique

# discrétisation des variables continues
npred <- -grep('(Cle|Duree_credit|Montant_credit|Age)', names(credit))
credit2 <- credit[,npred]
credit2$Age <- cut(credit$Age,c(0,25,Inf),right=TRUE)
credit2$Duree_credit <- cut(credit$Duree_credit,c(0,15,36,Inf),right=TRUE)
credit2$Montant_credit <- cut(credit$Montant_credit,c(0,4000,Inf),right=TRUE)

# choix de la modalité de référence pour l'âge
credit2$Age <- relevel(credit2$Age,ref="(25,Inf]")

# recodage et fusion des modalités
library(car)

# attention aux simples et doubles quotes dans le "recode"
credit2$Comptes <- recode(credit2$Comptes,"4='Pas de compte';1='CC < 0 euros';2='CC [0-200 euros[';3='CC > 200 euros' ")
#credit2$Comptes  <- factor(credit2$Comptes,levels=c('4','1','2','3'),labels=c("Pas de compte","CC < 0 euros","CC [0-200 euros[","CC  >= 200 euros"))
table(credit2$Comptes)

credit2$Historique_credit <- recode(credit2$Historique_credit,"c('A30','A31')='Crédits en impayé';
c('A32','A33')='Pas de crédits ou en cours sans retard';'A34'='Crédits passés sans retard'")
table(credit2$Historique_credit)

credit2$Objet_credit <- recode(credit$Objet_credit,"'A40'='Voiture neuve';
'A41'='Voiture occasion';c('A42','A43','A44','A45')='Intérieur';c('A46','A48','A49','A410')='Etudes-business-Autres';
'A47'='Vacances'")
table(credit2$Objet_credit)

credit2$Epargne <- recode(credit2$Epargne,
"c(0,3,4)='Pas épargne ou > 500 euros';1:2='< 500 euros'")
table(credit2$Epargne)

credit2$Anciennete_emploi<- recode(credit2$Anciennete_emploi,
"1:2='Sans emploi ou < 1 an';3='E [1-4[ ans';4:5='E GE 4 ans'")
table(credit2$Anciennete_emploi)

credit2$Situation_familiale<- recode(credit2$Situation_familiale,
"'A91'='Homme divorcé/séparé';'A92'='Femme divorcée/séparée/mariée';
c('A93','A94')='Homme célibataire/marié/veuf';'A95'='Femme célibataire'")
table(credit2$Situation_familiale)

credit2$Garanties <- recode(credit2$Garanties,"'A103'='Avec garant';else='Sans garant'")
table(credit2$Garanties)

credit2$Biens <- recode(credit2$Biens,"'A121'='Immobilier';'A124'='Aucun bien';
else='Non immobilier'")
table(credit2$Biens)

credit2$Autres_credits <- recode(credit2$Autres_credits,"'A143'='Aucun crédit extérieur';else='Crédits extérieurs'")
table(credit2$Autres_credits)

credit2$Statut_domicile <- recode(credit2$Statut_domicile,"'A152'='Propriétaire';else='Non Propriétaire'")
table(credit2$Statut_domicile)

# structure du fichier de crédit après transformations
summary(credit2)

# échantillons d'apprentissage et de validation
train  <- credit2[id,]
valid  <- credit2[-id,]

# régression logistique
logit <- glm(Cible~Comptes+Historique_credit+Duree_credit+Age+Epargne+Garanties+Autres_credits,data=train,family=binomial(link = "logit"))
summary(logit)

# grille de score
VARIABLE=c("",gsub("[0-9]", "", names(unlist(logit$xlevels))))
MODALITE=c("",as.character(unlist(logit$xlevels)))
MODALITE=c("",unlist(logit$xlevels))
names=data.frame(VARIABLE,MODALITE,NOMVAR=c("(Intercept)",paste(VARIABLE,MODALITE,sep="")[-1]))
regression=data.frame(NOMVAR=names(coefficients(logit)),COEF=as.numeric(coefficients(logit)))
param = merge(names,regression,all.x=TRUE)[-1]
param$COEF[is.na(param$COEF)] <- 0 
param
# calcul du poids total pour normalisation
mini=aggregate(data.frame(min = param$COEF), by = list(VARIABLE = param$VARIABLE), min)
maxi=aggregate(data.frame(max = param$COEF), by = list(VARIABLE = param$VARIABLE), max)
total=merge(mini,maxi)
total$diff = total$max - total$min
poids_total = sum(total$diff)
# calcul des poids par modalité
grille = merge(param,mini,all.x=TRUE)
grille$delta = grille$COEF - grille$min
grille$POIDS = round((100*grille$delta) / poids_total)
grille[which(VARIABLE!=""),c("VARIABLE","MODALITE","POIDS")]
grille[order(grille$VARIABLE,grille$MODALITE)[which(VARIABLE!="")],c("VARIABLE","MODALITE","POIDS")]

# application du modèle à un jeu de données
valid$logit <- predict(logit, newdata=valid, type="response")

# aire sous la courbe ROC
library(ROCR)
pred <- prediction(valid$logit,valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # AUC = 0,76202
#auc(valid$Cible,valid$logit) # avec package pROC

# courbe ROC
perf <- performance(pred,"tpr","fpr")
plot(perf,main='Courbe ROC')
segments(0,0,1,1,lty=3) # ajout diagonale en pointillés 
# courbe de lift
lift <- performance(pred,"lift","rpp")
plot(lift,main='Courbe de lift')
# courbe de lift usuelle
lift <- performance(pred,"tpr","rpp")
plot(lift,main='Courbe de lift')
segments(0,0,prop.table(table(credit$Cible))[2],1,lty=3) # ajout diagonale en pointillés gris
segments(prop.table(table(credit$Cible))[2],1,1,1,lty=3) # ajout diagonale en pointillés gris

# courbe ROC avec intervalles de confiance
library(pROC)
roc <- plot.roc(valid$Cible,valid$logit,main="", percent=TRUE, ci=TRUE)
# calcul des intervalles de confiance par simulation de Monte-Carlo
roc.se <- ci.se(roc,specificities=seq(0, 100, 5))
plot(roc.se, type="shape", col="grey")

# superposition des fonctions de répartition des bons et mauvais dossiers
plot.ecdf((valid$logit[valid$Cible==0]),main="Fonction de répartition du score",col="blue",pch=16)
plot.ecdf((valid$logit[valid$Cible==1]),col="red",pch=17,add=T)
legend("bottomright",c("Score=0","Score=1"),pch=c(16,17),col=c("blue","red"),lwd=1)

# Kolmogorov-Smirnov
perf <- performance(pred,"tpr", "fpr")
max(perf@y.values[[1]]-perf@x.values[[1]])
ks <- perf@y.values[[1]]-perf@x.values[[1]]
(seuil <- pred@cutoffs[[1]][which.max(ks)])
segments(seuil,1-perf@y.values[[1]][which.max(ks)],seuil,1-perf@x.values[[1]][which.max(ks)],col='black',lty=3)

# superposition des fonctions de densité des bons et mauvais dossiers
plot(density(valid$logit[valid$Cible==0]),main="Fonction de densité du score",col="blue",
xlim = c(-0.2,1.1), ylim = c(0,3),lwd=2)
par(new = TRUE)
plot(density(valid$logit[valid$Cible==1]),col="red",lty=3,lwd=2,
xlim = c(-0.2,1.1), ylim = c(0,3),xlab = '', ylab = '',main=' ')
legend("topright",c("Score=0","Score=1"),
lty=c(1,3),col=c("blue","red"),lwd=2)
abline(v=seuil,col="grey")

# box-plot
plot(logit~Cible,data=valid,horizontal = TRUE)

# application du modèle à l'ensemble des données
credit2$logit <- predict(logit, newdata=credit2, type="response")
pred <- prediction(credit2$logit,credit2$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # AUC = 0,7899667

# taux d'impayés par décile du score
q <- quantile(credit2$logit, seq(0, 1, by=0.05))
qscore <- cut(credit2$logit,q)
tab <- table(qscore,credit2$Cible)
ti <- prop.table(tab,1)[,2] # affichage % en ligne
par(mar = c(7, 4, 2, 0))
barplot(as.numeric(ti),col=gray(0:length(ti)/length(ti)),
names.arg=names(ti), ylab='Taux impayés', ylim=c(0,1),cex.names = 0.8, las=3)
abline(v=c(7.3,18.1),col="red")

library(car)
zscore <- recode(credit2$logit,as.factor.result=T,"lo:0.117='Faible';
0.117:0.486='Moyen';
0.486:hi='Fort'")
table(zscore)
tab <- table(zscore,credit2$Cible)
prop.table(tab,1) # affichage % en ligne

# discrétisation automatique du score
for (i in (3:4)) decoup(credit2,logit,Cible,nbmod=i,h=1,k=0,pAUC=0.8)
decoup(credit2,logit,Cible,nbmod=3,h=1,k=0,pAUC=0.8)



# ---------------------------------------------------------------------------------------------------------
# Régression logistique probit
# ---------------------------------------------------------------------------------------------------------

# probit
probit <- glm(Cible~Comptes+Historique_credit+Duree_credit+Age+Epargne+Garanties+Autres_credits,data=train,family=binomial(link = "probit"))
summary(probit)

logit$coefficients/probit$coefficients
pi/sqrt(3)

valid$probit <- predict(probit, newdata=valid, type="response")
pred <- prediction(valid$probit,valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # AUC = 0.7576286



# ---------------------------------------------------------------------------------------------------------
# Régression logistique log-log
# ---------------------------------------------------------------------------------------------------------

# log-log
cloglog <- glm(Cible~Comptes+Historique_credit+Duree_credit+Age+Epargne+Garanties+Autres_credits,data=train,family=binomial(link = "cloglog"))
summary(cloglog)
# modèle sur cible inversée
ccloglog <- glm(ifelse(train$Cible==0,1,0)~Comptes+Historique_credit+Duree_credit+Age+Epargne+Garanties+Autres_credits,data=train,family=binomial(link = "cloglog"))
summary(ccloglog)

valid$cloglog <- predict(cloglog, newdata=valid, type="response")
pred <- prediction(valid$cloglog,valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

cor(valid$logit,valid$probit,valid$cloglog)



# ---------------------------------------------------------------------------------------------------------
# CHAPITRE 6
# Régression logistique pénalisée ridge
# ---------------------------------------------------------------------------------------------------------

# régression elasticnet (méthode de descente des coordonnées)
library(glmnet)

# le vecteur de prédicteurs x de glmnet doit être numérique => toutes
# les variables qui sont des facteurs sont remplacées par les indicatrices
# de leurs modalités (sauf la modalité de référence)
x <- model.matrix( ~ . -1 ,data=train[,-which(names(train)=="Cible")])
y <- train[,"Cible"]

# validation croisée
set.seed(235)
cvfit = cv.glmnet(x,y,alpha=0.001, family = "binomial",type="auc",nlambda=100)
cvfit$lambda # valeurs de pénalisation
cvfit$lambda[1] # plus petit lambda annulant tous les coefficients
# à noter que ce lambda devrait être infini pour alpha = 0 (ridge)
# et que pour avoir un lambda fini, on remplace alpha = 0 par alpha = 0.001
cvfit$lambda[100] # lambda précédent divisé par 10000
cvfit$lambda.min # lambda donnant le plus petit taux d’erreur (cross-validé)
which(cvfit$lambda==cvfit$lambda.min) # valeur de lambda donnant l'AUC maximale
# représentation graphique de l'erreur (= AUC) en fonction de la pénalisation
plot(cvfit)
abline(h=cvfit$cvm[which(cvfit$lambda==cvfit$lambda.min)],col='blue',lty=2)
abline(v=log(cvfit$lambda.min),col='blue',lty=2)

# calcul de la régression pour une plage de valeurs de lambda
fit = glmnet(x,y,alpha=0,family = "binomial",lambda=seq(cvfit$lambda[1],
cvfit$lambda[100],length=10000),standardize = T)

# affichage des coefficients
plot(fit)
plot(fit,xvar="lambda",label="T")

# prédiction sur une plage de lambda
# sans précision de "s", toute la séquence de lambda utilisée en apprentissage est reprise
ypred <- predict(fit,newx=x,type="response")
length(fit$lambda)
library(ROCR)
roc <- function(x) { performance(prediction(ypred[,x],y),"auc")@y.values[[1]] }
vauc <- Vectorize(roc)(1:length(fit$lambda))
which.max(vauc) # rang de la pénalisation donnant la plus forte AUC
fit$lambda[which.max(vauc)] # pénalisation donnant la plus forte AUC
vauc[which.max(vauc)] # plus forte AUC

# prédiction sur une plage de lambda sur la base de test
xt <- model.matrix( ~ . -1,data=valid[,-which(names(valid)=="Cible")])
yt <- as.numeric(valid[,"Cible"])
ytpred <- predict(ridge,newx=xt,type="response")
vauc <- Vectorize(roc)(1:length(ridge$lambda))
# affichage de l'AUC
plot(vauc~log(ridge$lambda),col='blue',lty=2,cex=0.5,pch=16)
#abline(v=log(ridge$lambda[which.max(vauc)]),col='black',lty=2)
abline(h=vauc[which.max(vauc)],col='black',lty=2)
vauc[which.max(vauc)] # 0.7806553
which.max(vauc)
ridge$lambda[which.max(vauc)] # pénalisation donnant la plus forte AUC
log(ridge$lambda[which.max(vauc)]) # log-pénalisation donnant la plus forte AUC
tail(ridge$lambda[which(vauc >= 0.780655)])
vauc[which(ridge$lambda==1.5)] # 0.7806553

# affichage des coefficients
plot(fit,xvar="lambda",label="T")
abline(v=log(1.5),col='black',lty=2)

# coefficients donnant l'AUC maximale en test
coef(ridge,s=ridge$lambda[which.max(vauc)])
coef(ridge,s=1.5)

# équivalent au logit$xlevels précédent 
niveaux <- Vectorize(levels)(train[,1:ncol(train)])

# traitement des variables quantitatives
# auxquelles on ajoute une modalité unique fictive
niveaux$Taux_effort <- ""
niveaux$Anciennete_domicile <- ""
niveaux$Nb_credits <- ""
niveaux

# grille de score
VARIABLE=c("",gsub("[0-9]", "", names(unlist(niveaux))))
MODALITE=c("",as.character(unlist(niveaux)))
names=data.frame(VARIABLE,MODALITE,NOMVAR=c("(Intercept)",paste(VARIABLE,MODALITE,sep="")[-1]))
#coef <- ridge$beta[,which.max(vauc)]
coef <- ridge$beta[,which(ridge$lambda==1.5)]
regression=data.frame(NOMVAR=names(coef),COEF=as.numeric(coef))
param = merge(names,regression,all.x=TRUE)[-1]
param$COEF[is.na(param$COEF)] <- 0
param
# calcul du poids total pour normalisation
mini=aggregate(data.frame(min = param$COEF), by = list(VARIABLE = param$VARIABLE), min)
maxi=aggregate(data.frame(max = param$COEF), by = list(VARIABLE = param$VARIABLE), max)
total=merge(mini,maxi)
total$min[total$min==total$max] <- 0
total$diff = abs(total$max - total$min)
poids_total = sum(total$diff)
# calcul des poids par modalité
grille = merge(param,mini,all.x=TRUE)
grille$delta = grille$COEF - grille$min
grille$delta[grille$VARIABLE=="Anciennete_domicile"] <- 
abs(max(train$Anciennete_domicile)*grille$COEF[grille$VARIABLE=="Anciennete_domicile"])
grille$delta[grille$VARIABLE=="Nb_credits"] <- 
abs(max(train$Nb_credits)*grille$COEF[grille$VARIABLE=="Nb_credits"])
grille$delta[grille$VARIABLE=="Taux_effort"] <- 
abs(max(train$Taux_effort)*grille$COEF[grille$VARIABLE=="Taux_effort"])
grille$POIDS = round((100*grille$delta) / poids_total)
grille[which(grille$VARIABLE!=""&grille$VARIABLE!="Cible"),c("VARIABLE","MODALITE","POIDS")]


# calcul d'un modèle pénalisé après suppression des variables les moins discriminantes
# !!!! on applique auparavant les regroupements définitifs de modalités !!!!
# ceux du meilleur modèle logistique simple
credit2$Anciennete_domicile <- NULL
credit2$Nb_pers_charge <- NULL
credit2$Type_emploi <- NULL
credit2$Telephone <- NULL
credit2$Nb_credits <- NULL
train  <- credit2[id,]
valid  <- credit2[-id,]

x <- model.matrix( ~ . -1 ,data=train[,-which(names(train)=="Cible")])
y <- train[,"Cible"]

set.seed(235)
cvfit = cv.glmnet(x,y,alpha=0, family = "binomial",type="auc",nlambda=100)
fit = glmnet(x,y,alpha=0,family = "binomial",lambda=seq(cvfit$lambda[1],
cvfit$lambda[100],length=10000),standardize = T)
fit = glmnet(x,y,alpha=0,family = "binomial",lambda=seq(cvfit$lambda[1],
cvfit$lambda[100],length=10000),standardize = T)
# prédiction sur une plage de lambda sur la base de test
xt <- model.matrix( ~ . -1,data=valid[,-which(names(valid)=="Cible")])
yt <- valid[,"Cible"]
ytpred <- predict(fit,newx=xt,type="response")
roc <- function(x) { performance(prediction(ytpred[,x],yt),"auc")@y.values[[1]] }
vauc <- Vectorize(roc)(1:length(fit$lambda))
vauc[which.max(vauc)]
fit$lambda[which.max(vauc)] # pénalisation donnant la plus forte AUC
log(fit$lambda[which.max(vauc)])
plot(vauc~log(fit$lambda),col='blue',lty=2,cex=0.5,pch=16)
abline(v=log(fit$lambda[which.max(vauc)]),col='black',lty=2)
abline(h=vauc[which.max(vauc)],col='black',lty=3)
coef(fit,s=fit$lambda[which.max(vauc)])



# ---------------------------------------------------------------------------------------------------------
# CHAPITRE 7
# Régression logistique pénalisée lasso
# ---------------------------------------------------------------------------------------------------------

# régression elasticnet (méthode de descente des coordonnées)
library(glmnet)

# N.B.
# les résultats présentés dans le livre ont été obtenus
# sans regrouper la modalité A410 "Autres" avec "Etudes-Business" (variable Objet)
credit2$Anciennete_domicile <- NULL
credit2$Nb_pers_charge <- NULL
credit2$Type_emploi <- NULL
credit2$Telephone <- NULL
credit2$Nb_credits <- NULL
# les deux variables Biens et Situation_familiale ne sont supprimées
# que dans le cadre du "relaxed lasso"
#credit2$Biens <- NULL
#credit2$Situation_familiale <- NULL
# mais cela ne fait pas monter l'AUC
train  <- credit2[id,]
valid  <- credit2[-id,]

x <- model.matrix( ~ . -1 ,data=train[,-which(names(train)=="Cible")])
y <- as.numeric(train[,"Cible"])

# validation croisée
set.seed(235)
cvfit = cv.glmnet(x,y,alpha=1, family = "binomial",type="auc",nlambda=100)
cvfit$lambda[1] # plus petit lambda annulant tous les coefficients
coef(cvfit,s=0.1571349)
coef(cvfit,s=0.158)
cvfit$lambda[100] # lambda précédent divisé par 10000
cvfit$lambda.min # lambda donnant le plus petit taux d’erreur (cross-validé)
which(cvfit$lambda==cvfit$lambda.min) # 83 c'est la 83e valeur de lambda
# sur 100 qui donne l'AUC maximale
# c'est une valeur de lambda plus proche de son minimum
# (et donc de la régression logistique) que de son maximum (annulant tous les coefficients)
# sachant qu'il n'y a que 84 valeurs de lambda
cvfit$cvm[which(cvfit$lambda==cvfit$lambda.min)]
# plus faible erreur (i.e. plus forte AUC) = 0.8054402
# représentation graphique de l'erreur (= AUC) en fonction de la pénalisation
plot(cvfit)
abline(h=cvfit$cvm[which(cvfit$lambda==cvfit$lambda.min)],col='blue',lty=2)
abline(v=log(cvfit$lambda.min),col='blue',lty=2)

# calcul de la régression pour une plage de valeurs de lambda
fit = glmnet(x,y,alpha=1,family = "binomial",lambda=seq(cvfit$lambda[1],
cvfit$lambda[length(cvfit$lambda)],length=10000),standardize = T)
coef(fit,s=log(cvfit$lambda[length(cvfit$lambda)]))
CrossTable(credit2$Objet_credit,credit2$Cible,prop.chisq=F,chisq = T,format="SAS")

# affichage des coefficients
plot(fit,xvar="lambda",label="T")

# prédiction sur une plage de lambda sur la base de test
xt <- model.matrix( ~ . -1,data=valid[,-which(names(valid)=="Cible")])
yt <- valid[,"Cible"]
ytpred <- predict(fit,newx=xt,type="response")
library(ROCR)
roc <- function(x) { performance(prediction(ytpred[,x],yt),"auc")@y.values[[1]] }
vauc <- Vectorize(roc)(1:length(fit$lambda))

plot(vauc~log(fit$lambda),col='blue',lty=2,cex=0.5,pch=16)
abline(v=log(fit$lambda[which.max(vauc)]),col='black',lty=2)
abline(h=vauc[which.max(vauc)],col='black',lty=3)

# pénalisation optimale
which.max(vauc) # rang de la pénalisation donnant la plus forte AUC
fit$lambda[which.max(vauc)] # pénalisation donnant la plus forte AUC
log(fit$lambda[which.max(vauc)])

# Kolmogorov-Smirnov
performance(prediction(ytpred[,which.max(vauc)],yt),"auc")@y.values[[1]]
pred <- prediction(ytpred[,which.max(vauc)],yt,label.ordering=c(0,1))
perf <- performance(pred,"tpr", "fpr")
max(attr(perf,'y.values')[[1]]-attr(perf,'x.values')[[1]])

# affichage des coefficients
plot(fit,xvar="lambda",label="T")
abline(v=log(fit$lambda[which.max(vauc)]),col='black',lty=2)

# coefficients donnant l'AUC maximale en test
(coef <- coef(fit,s=fit$lambda[which.max(vauc)]))
sum(coef[,1]==0) # 2
sum(coef[,1]!=0) # 24

# grille de score lasso
# équivalent au logit$xlevels précédent 
niveaux <- Vectorize(levels)(train[,1:ncol(train)])
niveaux$Taux_effort <- ""
niveaux

# grille de score
VARIABLE=c("",gsub("[0-9]", "", names(unlist(niveaux))))
MODALITE=c("",as.character(unlist(niveaux)))
names=data.frame(VARIABLE,MODALITE,NOMVAR=c("(Intercept)",paste(VARIABLE,MODALITE,sep="")[-1]))
(coef <- fit$beta[,which.max(vauc)])
regression=data.frame(NOMVAR=names(coef),COEF=as.numeric(coef))
regression
param = merge(names,regression,all.x=TRUE)[-1]
param$COEF[is.na(param$COEF)] <- 0
param
# calcul du poids total pour normalisation
mini=aggregate(data.frame(min = param$COEF), by = list(VARIABLE = param$VARIABLE), min)
maxi=aggregate(data.frame(max = param$COEF), by = list(VARIABLE = param$VARIABLE), max)
total=merge(mini,maxi)
total$min[total$min==total$max] <- 0
total
total$diff = abs(total$max - total$min)
poids_total = sum(total$diff)
poids_total
# calcul des poids par modalité
grille = merge(param,mini,all.x=TRUE)
grille
grille$delta = grille$COEF - grille$min
grille$delta[grille$VARIABLE=="Taux_effort"] <- 
abs(max(train$Taux_effort)*grille$COEF[grille$VARIABLE=="Taux_effort"])
grille$POIDS = round((100*grille$delta) / poids_total)
grille[which(grille$VARIABLE!=""&grille$VARIABLE!="Cible"),c("VARIABLE","MODALITE","POIDS")]

# coefficients des 10000 modèles générés pour les différentes pénalisations
coef <- as.matrix(coef(fit))
#coefnul = function(x){sum(coef[,x]==0)}
coefnnul = function(x){sum(coef[,x]!=0)-1} # on soustrait "1" pour l'intercept
coefn0 <- Vectorize(coefnnul)(1:length(fit$lambda))
lasso <- cbind(coefn0,vauc)
head(lasso)
tail(lasso)

# AUC max pour chaque nb de coef nuls
aggregate(lasso[,"vauc"],by=list(lasso[,"coefn0"]),max)
auccoef0 <- aggregate(lasso[,"vauc"],by=list(lasso[,"coefn0"]),max)
barplot(auccoef0$x,names.arg=auccoef0$Group.1,las=3,ylim=c(0,0.8))

# application sur l'ensemble de la base
xt <- model.matrix( ~ . -1,data=credit2[,-which(names(valid)=="Cible")])
yt <- credit2[,"Cible"]
ytpred <- predict(fit,newx=xt,type="response",s=fit$lambda[which.max(vauc)])
performance(prediction(ytpred,yt),"auc")@y.values[[1]] # 0.8160429



# ---------------------------------------------------------------------------------------------------------
# Adaptive lasso
# ---------------------------------------------------------------------------------------------------------

# régression elasticnet (méthode de descente des coordonnées)
library(glmnet)

# calcul d'un modèle pénalisé après suppression des variables les moins discriminantes
# !!!! on applique auparavant les regroupements définitifs de modalités !!!!
# ceux du meilleur modèle logistique simple
credit2$Anciennete_domicile <- NULL
credit2$Nb_pers_charge <- NULL
credit2$Type_emploi <- NULL
credit2$Telephone <- NULL
credit2$Nb_credits <- NULL
train  <- credit2[id,]
valid  <- credit2[-id,]

# échantillon d'apprentissage
x <- model.matrix( ~ . ,data=train[,-which(names(train)=="Cible")])
y <- train[,"Cible"]
# échantillon de validation
xt <- model.matrix( ~ . ,data=valid[,-which(names(valid)=="Cible")])
yt <- valid[,"Cible"]

# coefficients des moindres carrés ordinaires
beta.ols <- lm(y~x,data=train)$coefficients[-1]
beta.ols
penal <- pmax(1,abs(beta.ols)^(-1),na.rm = TRUE)
penal

# validation croisée
set.seed(235)
cvfit = cv.glmnet(x,y,alpha=1,family = "binomial",type="auc",nlambda=100)
# calcul de la régression pour une plage de valeurs de lambda
fit = glmnet(x,y,alpha=1,family = "binomial",lambda=seq(cvfit$lambda[1],
cvfit$lambda[length(cvfit$lambda)],length=10000),standardize = T,
penalty.factor=penal)
# prédiction sur une plage de lambda sur la base de test
ytpred <- predict(fit,newx=xt,type="response")
roc <- function(x) { performance(prediction(ytpred[,x],yt),"auc")@y.values[[1]] }
vauc <- Vectorize(roc)(1:length(fit$lambda))
vauc[which.max(vauc)]
coef(fit,s=fit$lambda[which.max(vauc)])



# ---------------------------------------------------------------------------------------------------------
# Group lasso
# ---------------------------------------------------------------------------------------------------------

# régression elasticnet (méthode de descente des coordonnées)
library(grplasso)

# group lasso avec détection automatique des groupes formées par les modalités des facteurs
summary(credit2)
length(id) # 644
0.1571349 * 644
lambda <- 0.003673408*length(id)
# pénalisation donnant précédemment (avec toutes les discrétisation) la plus forte AUC en test
lambda #  2.365675
lambda <- 0.003651067*length(id)
# pénalisation donnant précédemment la plus forte AUC en test
lambda #  2.351287
0.1571349*length(id) # plus petit lambda annulant précédemment tous les coefficients
lambda <- seq(100,1,length=1000)
# mise au format numérique de la variable à expliquer
# et transformation du codage 1/2 en 0/1
credit3 <- credit2
credit3$Cible <- as.numeric(credit2$Cible)-1
# group lasso sur une plage de valeur, en commençant par la plus forte pénalisation
gplasso <- grplasso(Cible ~ . ,data = credit3[id,vars], lambda = lambda, model = LogReg(), center=TRUE, standardize=TRUE)
coe <- coefficients(gplasso)
	
# affichage des variables dans l'ordre de leur sélection
lambda[1]
coe[,370]
modannul = function(i){sort(rownames(coe)[(coe[,i]!=0)])}
modannul(1)
for (j in (length(lambda)-1):1)
	{
	if (length(modannul(j)) < length(modannul(j+1)))
		{cat("\n",j,modannul(j+1)[-match(modannul(j),modannul(j+1))])}
	}


## coefficients en fonction de la pénalisation
plot(gplasso)
#plot(gplasso, log = "x")
abline(v=lambda[370],lty=3)
plot(gplasso, log = "x")


# autre essai de group lasso

# 1ère variable (intercept) not penalized => valeur NA
group <- c(NA, 2, 2, 2, 3, 3, 4, 4, 4, 4, 5, 6, 6, 7,
8, 8, 9, 10, 10, 11, 12, 13, 14, 14, 15)
group <- c(NA, 2, 2, 2, 3, 3, 4, 4, 4, 4, 5, 6, 6, 7,
8, 8, 9, 10, 10, 11, 12, 13, 14, 14, 15, 16, 17, 17, 18,
19, 20, 20, 21, 21, 22, 23, 23, 24)
# sans intercept
group <- c(2, 2, 2, 2, 3, 3, 4, 4, 4, 4, 5, 6, 6, 7,
8, 8, 9, 10, 10, 11, 12, 13, 14, 14, 15, 16, 17, 17, 18,
19, 20, 20, 21, 21, 22, 23, 23, 24)

gplasso <- grplasso(x, y , index = group, lambda = lambda, model = LogReg(), center=FALSE, standardize = TRUE)
coefficients(gplasso)
# avec le lambda optimal du lasso ordinaire
gplasso <- grplasso(x, y , index = group, lambda = 0.003651067*length(y), model = LogReg(), center=TRUE, standardize=TRUE)
coefficients(gplasso)



# ---------------------------------------------------------------------------------------------------------
# Régression logistique pénalisée elastic net
# ---------------------------------------------------------------------------------------------------------

# régression elasticnet (méthode de descente des coordonnées)
library(glmnet)

# calcul d'un modèle pénalisé après suppression des variables les moins discriminantes
# !!!! on applique auparavant les regroupements définitifs de modalités !!!!
# ceux du meilleur modèle logistique simple
credit2$Anciennete_domicile <- NULL
credit2$Nb_pers_charge <- NULL
credit2$Type_emploi <- NULL
credit2$Telephone <- NULL
credit2$Nb_credits <- NULL
train  <- credit2[id,]
valid  <- credit2[-id,]

# échantillon d'apprentissage
x <- model.matrix( ~ . -1 ,data=train[,-which(names(train)=="Cible")])
y <- train[,"Cible"]
# échantillon de validation
xt <- model.matrix( ~ . -1,data=valid[,-which(names(valid)=="Cible")])
yt <- valid[,"Cible"]

# calcul de l'AUC maximale en test, pour chaque valeur de alpha
elastic <- function(a)
{
set.seed(235)
cvfit = cv.glmnet(x,y,alpha=a, family = "binomial",type="auc",nlambda=100)
# calcul de la régression pour une plage de valeurs de lambda
fit = glmnet(x,y,alpha=a,family = "binomial",lambda=seq(cvfit$lambda[1],
cvfit$lambda[length(cvfit$lambda)],length=10000),standardize = T)
# prédiction sur une plage de lambda sur la base de test
ytpred <- predict(fit,newx=xt,type="response")
roc <- function(x) { performance(prediction(ytpred[,x],yt),"auc")@y.values[[1]] }
vauc <- Vectorize(roc)(1:length(fit$lambda))
return(vauc[which.max(vauc)])
}
elastic(0)
elastic(0.001)
elastic(0.1)
elastic(1)



# ---------------------------------------------------------------------------------------------------------
# CHAPITRE 8
# Régression logistique PLS
# ---------------------------------------------------------------------------------------------------------

install.packages('plsRglm')
install.packages('fields')
library(fields)
library(plsRglm)
train  <- credit2[id,]
valid  <- credit2[-id,]
summary(train)

# échantillon d'apprentissage
x <- model.matrix( ~ . -1 ,data=train[,-which(names(train)=="Cible")])
y <- as.numeric(train[,"Cible"])
# recodage en 0/1
y[y == 1] <- 0 # crédits OK
y[y == 2] <- 1 # crédits KO
# échantillon de validation
xt <- model.matrix( ~ . -1,data=valid[,-which(names(valid)=="Cible")])
yt <- as.numeric(valid[,"Cible"])
# recodage en 0/1
yt[yt == 1] <- 0 # crédits OK
yt[yt == 2] <- 1 # crédits KO

#pls1 <- plsRglm(y,x,nt=1,modele="pls-glm-logistic")
pls1 <- plsRglm(y,x,nt=1,dataPredictY=xt,modele="pls-glm-logistic")
pls2 <- plsRglm(y,x,nt=2,dataPredictY=xt,modele="pls-glm-logistic")
pls3 <- plsRglm(y,x,nt=3,dataPredictY=xt,modele="pls-glm-logistic")
pls4 <- plsRglm(y,x,nt=4,dataPredictY=xt,modele="pls-glm-logistic")
plsm <- plsRglm(y,x,nt=36,dataPredictY=xt,modele="pls-glm-logistic",pvals.expli=T,sparseStop=T)

# quelques sorties
plsm$computed_nt
pls1$InfCrit
barplot(pls1$InfCrit[,2])
print(pls1)
pls1$Std.Coeffs
pls1$Coeffs
print(plsm)
str(pls1)
pls1$family

# test sur variables centrées-réduites
summary(scale(x))
pls1b <- plsRglm(y,scale(x),nt=1,modele="pls-glm-logistic")
pls1b$Std.Coeffs
pls1b$Coeffs

# comparaison des signes des coefficients
cor(x,y)
cor(train$Taux_effort,as.numeric(train$Cible))
cor(x,y)*pls1$Coeffs[-1]
cor(x,y)*pls2$Coeffs[-1]
pls1$Coeffs
pls2$Coeffs
pls1$Coeffs*pls2$Coeffs
# variables avec des signes différents pour les coefficients de PLS1 et PLS2
rownames(pls1$Coeffs)[which(pls1$Coeffs*pls2$Coeffs<0)]
# variables avec des signes différents pour les coefficients de PLS2
# et la corrélation de la variable avec la variable Y à expliquer
# à noter que "+1" qui précède le "which" à cause de l'intercept manquant
# dans cor(x,y) et dans pls2$Coeffs[-1]
rownames(pls2$Coeffs)[1+which(cor(x,y)*pls2$Coeffs[-1]<0)]
rownames(pls1$Coeffs)[1+which(cor(x,y)*pls1$Coeffs[-1]<0)]
# coefficients opposés pour 1 et 2 composantes
pls1$Coeffs[which(pls1$Coeffs*pls2$Coeffs<0)]
pls2$Coeffs[which(pls1$Coeffs*pls2$Coeffs<0)]
cbind(rownames(pls1$Coeffs),pls1$Coeffs[1],pls2$Coeffs[1])[which(pls1$Coeffs*pls2$Coeffs<0)]
# la modalité Intérieur de Objet_credit est plus risquée que la moyenne
prop.table(table(train$Objet_credit,train$Cible),1)

# application du modèle
valid$pls1 <- pls1$ValsPredictY
valid$pls2 <- pls2$ValsPredictY
valid$pls3 <- pls3$ValsPredictY
valid$pls4 <- pls4$ValsPredictY

# calcul AUC
library(ROCR)
pred <- prediction(pls1$ValsPredictY,valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# enregistrement des coefficients
write.table(pls1$Coeffs, file = "plslogit.csv", quote = FALSE, sep = ";", na = " ", row.names = TRUE, col.names = NA)

# bootstrap pour avoir des intervalles de confiance des coefficients
# c'est un bootstrap dont les coefficients sont estimés par moindres carrés
bootpls1 <- bootplsglm(plsRglm(y,x,nt=1,modele="pls-glm-logistic"), sim="ordinary", stype="i")
bootpls1 <- bootplsglm(pls1, typeboot="fmodel_np", R=1000, statistic=coefs.plsRglmnp, sim="ordinary", stype="i")

# bootstrap pour avoir des intervalles de confiance des coefficients :
# c'est un bootstrap dont les coefficients sont estimés par la loi désignée lors la régression PLS
# et qui dans le cas du logit va donc effectuer les estimations par maximisation de la vraisemblance
ptm <- proc.time()
set.seed(123)
bootpls10000 <- bootplsglm(pls1, typeboot = "plsmodel", R = 10000, sim = "ordinary", stype = "i")
proc.time() - ptm
bootpls1

library(boot)
bootpls1 <- boot(data=cbind(y,x), statistic=coefs.plsRglm, sim="ordinary", stype="i", R=250, nt=3, modele="pls-glm-logistic")
bootpls1 <- bootplsglm(plsm,sim="ordinary", stype="i", R=1000)
bootpls1
bootpls2 <- bootplsglm(pls2,sim="ordinary", stype="i", R=1000)
warnings()
bootpls2
summary(coefs.plsRglmnp)

# intervalles de confiance
library(boot)
boot.ci(bootpls1, conf = c(0.90,0.95,0.99), type = c("norm","basic","perc","bca"), index=1)
boot.ci(bootpls1, conf = c(0.90,0.95,0.99), type = "all", index=37)
# Bootstrap confidence intervals à 95 % (normal, basic, percentile, et en option Bca)
# les résultats sont les mêmes que par la fonction boot.ci, mais sont fournis pour l'ensemble des prédicteurs
(ic.pls <- confints.bootpls(bootpls10000,typeBCa=TRUE))
# calculs des IC
str(bootpls1)
qnorm(0.95)
mean(bootpls1$t[,37]) # moyenne des coefficients bootstrapés d'un prédicteur
sd(bootpls1$t[,37]) # écart-type des coefficients bootstrapés d'un prédicteur
range(quantile(bootpls1$t[,22]))
?quantile
quantile(bootpls1$t[,22],c(0.90,0.95,0.975,0.99)) # quantiles des coefficients bootstrapés
# IC bootstrap normal
# à noter que l'IC normal n'est calculé ni à partir de la moyenne des coefs bootstrappés
# ni à partir du coefficient exact
bootpls1$t0[37] - (qnorm(0.975)*sd(bootpls1$t[,37]))
bootpls1$t0[37] + (qnorm(0.975)*sd(bootpls1$t[,37]))
mean(bootpls1$t[,37]) - (2*sd(bootpls1$t[,37]))
mean(bootpls1$t[,37]) + (2*sd(bootpls1$t[,37]))
# mais l'IC normal est calculé à partir du coefficient exact - le biais
# le biais étant la différence entre la moyenne des coefs bootstrappés et le coefficient exact
# l'IC normal est donc donné par :
(2*bootpls1$t0[37]) - mean(bootpls1$t[,37]) - (qnorm(0.975)*sd(bootpls1$t[,37]))
(2*bootpls1$t0[37]) - mean(bootpls1$t[,37]) + (qnorm(0.975)*sd(bootpls1$t[,37]))
# à noter qu'une variante calculée en remplaçant l'écart-type des coefs bootstrappés
# par la racine carrée de l'estimation bootstrap de la variance
# donne un résultat très légèrement différent
(bootvar <- mean(sapply(1:10000, function(i) (bootpls1$t[i,37] - mean(bootpls1$t[,37]))^2)))
(2*bootpls1$t0[37]) - mean(bootpls1$t[,37]) - (qnorm(0.975)*sqrt(bootvar))
(2*bootpls1$t0[37]) - mean(bootpls1$t[,37]) + (qnorm(0.975)*sqrt(bootvar))

# IC bootstrap basique (bootstrap pivotal confidence interval) 
# N.B. C'est le type 6 de calcul de quantiles qui permet de retrouver
# exactement le calcul du la fonction boot.ci pour les méthodes basique et percentile
(2*bootpls1$t0[37]) - quantile(bootpls1$t[,37],0.975,type=6)
(2*bootpls1$t0[37]) - quantile(bootpls1$t[,37],0.025,type=6)
# IC bootstrap percentile
quantile(bootpls1$t[,37],c(0.025,0.975),type=6)
# Boxplot plots for bootstrap distributions
boxplots.bootpls(bootpls1,prednames=FALSE,articlestyle=FALSE)
boxplots.bootpls(bootpls1,prednames=FALSE)
warnings()
help(boxplots.bootpls)
# Plot bootstrap confidence intervals
# Le choix typeBCa=T allonge fortement le temps de calcul
plots.confints.bootpls(confints.bootpls(bootpls1,typeBCa=T),legendpos = "bottomright")



# ---------------------------------------------------------------------------------------------------------
# CHAPITRE 9
# Arbres de décision
# ---------------------------------------------------------------------------------------------------------

library(tree)
arbre <- tree(Cible ~ Age+Duree_credit,data=credit)
arbre

# arbre de décision CART
install.packages("rpart.plot")
library(rpart)
library(rpart.plot)

# paramètres de rpart
# methode = class ou anova pour arbre de régression
# minsplit = 20 par défaut
# minbucket = 1/3*minsplit par défaut

set.seed(235)
cart <- rpart(Cible ~ . ,data = credit[id,vars],method="class",parms=list(split="gini"),cp=0)
#cart <- rpart(Cible ~ . ,data = credit[id,vars],method="class",parms=list(split="information"),cp=0)
cart <- rpart(Cible ~ . ,data = credit[id,vars],method="class",parms=list(split="gini"),cp=0.05)
cart <- rpart(Cible ~ . ,data = credit[id,vars],method="class",parms=list(split="gini"),control=list(minbucket=30,minsplit=30*2,maxdepth=2))

# scissions au milieu des intervalles des variables continues
cart <- rpart(Cible ~ Montant_credit ,data = credit[id,vars],method="class",parms=list(split="gini"),cp=0.01)
table(credit$Montant_credit)

cart # commande équivalente à "print(cart)"
summary(cart,digits=3) # plus d'informations sur les scissions

# affichage graphique de l'arbre
plot(cart,branch=.2, uniform=T, compress=T, margin=.1)
text(cart, fancy=T,use.n=T,pretty=0,all=T, cex=1)

# affichage graphique un peu amélioré de l'arbre
plot(cart,branch=.2, uniform=T, compress=T, margin=.1) # tracé de l'arbre
text(cart, fancy=T,use.n=T,pretty=0,all=T, cex=.6) # ajout des légendes des noeuds
# cex est un facteur multiplicateur de la taille des caractères
# création d'un fichier graphique au format postscript
post(cart, file = "C:/Users/.../Etude de cas/CART.ps",
title = "CART pour credit scoring")
# affichage amélioré avec package rpart.plot
prp(cart,type=2,extra=1,split.box.col="lightgray")

install.packages("rpart.utils")
library(rpart.utils)
rpart.rules.table(cart)

# calcul erreur par resusbstitution
sum(predict(cart,type="class") != credit[id,"Cible"])/nrow(credit[id,])
0.4766839 * 193/644
0.9637306 * 193/644
# calcul écart-type de l'erreur xerror
x <- 0.9637306*0.29969
sqrt((x*(1-x))/644)/0.29969

# on peut afficher l'erreur xerror calculée par validation croisée
# et le nb nsplit de scissions
# en fonction du coefficient de complexité
set.seed(235)
printcp(cart)
print(cart$cptable)
(28*0.002590674)+0.4766839
(26*0.002590674)+0.4818653
(0.059585492*2)+0.8808290 # = 1

# ajout de l'AUC à côté de l'erreur xerror
library(ROCR)
set.seed(235)
auc <- matrix(NA,nrow(cart$cptable)-1,4)
for(i in 2:nrow(cart$cptable))
{
cartp <- prune(cart,cp=cart$cptable[i,"CP"])
predc <- predict(cartp,type="prob",test)[,2]
pred <- prediction(predc,test$Cible,label.ordering=c(0,1))
auc[i-1,1] <- cart$cptable[i,"CP"]
auc[i-1,2] <- cart$cptable[i,"nsplit"]+1
auc[i-1,3] <- cart$cptable[i,"xerror"]
auc[i-1,4] <- performance(pred,"auc")@y.values[[1]]
} # fin de la boucle
colnames(auc) <- c("CP","nfeuilles","erreur","AUC")
auc

# on peut aussi afficher cela dans un graphique
# montrant la décroissance erreur relative en fonction du coefficient de complexité
plotcp(cart) # calcul par validation croisée
summary(cart,digits=3)

# élaguage
prunedcart = prune(cart,cp=0.0129534)
prunedcart4f = prune(cart,cp=0.0328152)
prunedcart7f = prune(cart,cp=0.0310881)
prunedcart9f = prune(cart,cp=0.0155441)
prunedcart15f = prune(cart,cp=0.0103627)
# choix valeur coefficient de pénalisation de la complexité pour élagage
plot(prunedcart4f,branch=.2, uniform=T, compress=T, margin=.1)
text(prunedcart4f, fancy=T,use.n=T,pretty=0,all=T,cex=.5)
# affichage avec package rpart.plot
prp(prunedcart4f,type=2,extra=1,split.box.col="lightgray")

# élaguage automatique au minimum d'erreur + 1 écart-type
xerr <- cart$cptable[,"xerror"]
xerr
minxerr <- which.min(xerr)
seuilerr <- cart$cptable[minxerr, "xerror"]+cart$cptable[minxerr, "xstd"]
xerr [xerr < seuilerr][1]
mincp <- cart$cptable[names(xerr [xerr < seuilerr][1]), "CP"]
prunedcart <- prune(cart,cp=mincp)
names(xerr [xerr < seuilerr][1])

# élaguage automatique au minimum d'erreur
xerr <- cart$cptable[,"xerror"]
minxerr <- which.min(xerr)
mincp <- cart$cptable[minxerr, "CP"]
prunedcart <- prune(cart,cp=mincp)
# code plus compact :
prunedcart <- prune(cart, cp=cart$cptable[which.min(cart$cptable[,"xerror"]),"CP"])
# affichage de l'arbre élagué
plot(prunedcart,branch=.2, uniform=T, compress=T, margin=.1)
text(prunedcart, fancy=T,use.n=T,pretty=0,all=T,cex=.5)

# élagage au nombre minimum de feuilles
mincart <- prune(cart,cp=cart$cptable[min(2,nrow(cart$cptable)), "CP"])
plot(mincart,branch=.2, uniform=T, compress=T, margin=.1)
text(mincart, fancy=T,use.n=T,pretty=0,all=T,cex=.5)
prp(mincart,type=2,extra=1,split.box.col="lightgray")

# for stumps : rpart.control(maxdepth=1,cp=-1,minsplit=0)
cart <- rpart(Cible ~ . ,data = credit[id,vars],method="class",parms=list(split="gini"),control=list(maxdepth=1,cp=-1,minsplit=0))
plot(cart)

# mesure des performances sur l'échantillon de test
test$CART <- predict(cart,type="prob",test)
test$CART <- predict(prunedcart,type="prob",test)
test$CART4f <- predict(prunedcart4f,type="prob",test)
test$CART7f <- predict(prunedcart7f,type="prob",test)
test$CART9f <- predict(prunedcart9f,type="prob",test)
test$mincart <- predict(mincart,type="prob",test)
head(test["CART4f"],5)

# aire sous la courbe ROC
library(ROCR)
pred <- prediction(test$CART9f[,2],test$Cible,label.ordering=c(0,1))
auc <- performance(pred,"auc")
performance(pred,"auc")@y.values[[1]]
auc
# équivalent : attr(performance(pred,"auc"),'y.values')[[1]]
# récupère l'attribut "y.values" de l'objet "performance"
perf <- performance(pred,"tpr","fpr")
plot(perf,main='Courbe ROC')
segments(0,0,1,1,col='grey',lty=3) # ajout diagonale en pointillés gris
# courbe de lift
perf <- performance(pred,"lift","rpp")
plot(perf,main='Courbe de lift')
# courbe de lift usuelle
pred <- prediction(ytpred,yt,label.ordering=c(0,1))
lift <- performance(pred,"tpr","rpp")
plot(lift,main='Courbe de lift')

# comparaison de courbes ROC
pred7f <- prediction(test$CART7f[,2],test$Cible,label.ordering=c(0,1))
pred9f <- prediction(test$CART9f[,2],test$Cible,label.ordering=c(0,1))
plot(performance(pred9f,"tpr","fpr"),col='red',lty=2,main='AUC de modèles CART')
plot(performance(pred7f,"tpr","fpr"), col='black',add=TRUE,lty=1)
segments(0,0,1,1,lty=3)
#legend(0.6,0.6,c('7 feuilles','9 feuilles'),col=c('red','black'),lwd=3)
legend("bottomright",c('9 feuilles','7 feuilles'),col=c('red','black'),
lty=c(2,1),lwd=3)

# comparaison de courbes ROC avec intervalles de confiance
library(pROC)
roc <- plot.roc(test$Cible,test$CART7f[,2],col='black',lty=1,ci=TRUE)
plot.roc(test$Cible,test$CART9f[,2],add=TRUE,col='red',lty=2,ci=TRUE)
# calcul des intervalles de confiance par simulation de Monte-Carlo
roc.se <- ci.se(roc,specificities=seq(0,1,.01),boot.n=10000)
plot(roc.se, type="shape", col="#0000ff22")
#plot(roc.se, type="shape", col=c(rgb(1, 1, 0, 0.2), rgb(0,1, 1, 0.2), rgb(1, 0, 1, 0.2)))
legend("bottomright",c('9 feuilles','7 feuilles'),col=c('red','black'),lty=c(2,1),lwd=3)

# gestion des valeurs manquantes
cart <- rpart(Cible ~ . ,data = credit[id,vars],method="class",parms=list(split="gini"),cp=0.05,control=list(usesurrogate=0,maxcompete=0))
summary(cart)
test2 <- test
ms <- sample(dim(test)[1],100)
test2$Objet_credit[ms] <- test2$Duree_credit[ms] <- test2$Comptes[ms] <- NA
summary(test2)
testNA <- predict(cart,type="prob",test2,na.action = na.omit)
summary(testNA)
head(test2$Comptes,100)
head(testNA,100)

# erreur
test$cart <- predict(prunedcart9f,type="class",test)
err_cart <- sum(test$cart != test$Cible) / nrow(test)
err_cart # 0.247 pour CART à 9 feuilles



# ---------------------------------------------------------------------------------------------------------
# Arbre de décision CHAID
# ---------------------------------------------------------------------------------------------------------

# install.packages("partykit")
install.packages("CHAID", repos="http://R-Forge.R-project.org")

library("CHAID")

credit2$Anciennete_domicile <- NULL
credit2$Taux_effort <- NULL
credit2$Nb_credits <- NULL

str(credit2)
# paramètres par défaut
chaid <- chaid(Cible ~ . ,data = credit2[id,vars])

# avec des paramètres de contrôle
chaid_control(alpha2 = 0.05, alpha3 = -1, alpha4 = 0.05,
               minsplit = 20, minbucket = 7, minprob = 0.01,
               stump = F, maxheight = -1)

ch.arbre <- chaid_control(alpha2 = 0.05, alpha3 = -1, alpha4 = 0.05,
               minsplit = 30, minbucket = 10, stump = F, maxheight = 3)

chaid <- chaid(Cible ~ . ,data = credit2[id,vars], control=ch.arbre)
chaid

# stump
ch.arbre <- chaid_control(alpha4 = 1, stump = TRUE)
chaid <- chaid(Cible ~ . ,data = credit2[id,vars], control=ch.arbre)
chaid

# affichage
plot(chaid)

# AUC
library(ROCR)
predc <- predict(chaid,type="prob",test)[,2]
pred <- prediction(predc,test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]



# ---------------------------------------------------------------------------------------------------------
# Arbre C5.0
# ---------------------------------------------------------------------------------------------------------

# arbre de décision C5.0
install.packages("C50")
library(C50)

c50 <- C5.0(Cible ~ . ,data = credit[id,vars],rules=T)
C5imp(c50)
summary(c50)

C50 <- predict(c50,type="prob",credit[-id,vars])
head(C50[,2])
pred <- prediction(C50[,2],credit[-id,"Cible"],label.ordering=c(0,1))
(auc <- performance(pred,"auc")@y.values[[1]]) # 0.6469617

c50 <- C5.0(Cible ~ . ,data = credit[id,vars],rules=T,
control = C5.0Control(minCases=50))
# mincases = 30 ou 50 => AUC = 0.6480689



# ---------------------------------------------------------------------------------------------------------
# CHAPITRE 10
# Algorithme PRIM
# ---------------------------------------------------------------------------------------------------------

library(prim)

# sur 2 variables
# ------------------------------------------

x <- credit[,c("Duree_credit","Age")]
y <- as.numeric(credit[,"Cible"])
y[y == 1] <- 0 # crédits OK
y[y == 2] <- 1 # crédits KO

prim <- prim.box(x, y, peel.alpha=0.05, paste.alpha=0.01,mass.min=0.05, pasting=T, verbose=T,threshold.type=1,threshold=0.3)
summary(prim,print.box=T)

#plot(prim,col = "transparent",xlim=range(credit$Duree_credit),ylim=range(credit$Age))
plot(prim,col = "transparent")
points(x[y == 1, ],col="red")
points(x[y == 0, ],col="green")
points(x)


# sur toutes les variables quantis
# ------------------------------------------

x <- credit[,varquanti]
y <- as.numeric(credit[,"Cible"])
y[y == 1] <- 0 # crédits OK
y[y == 2] <- 1 # crédits KO

prim <- prim.box(x, y, peel.alpha=0.05, paste.alpha=0.01,mass.min=0.05, pasting=T, verbose=T,threshold.type=1,threshold=0.3)
summary(prim,print.box=T)
prim$box


# sur toutes les variables quantis et indicatrices des variables qualis
# ------------------------------------------

formule <- as.formula(paste("y = ~ 0 +",paste(names(credit2[,varquali]),collapse="+")))

x <- cbind(credit[,varquanti],data.frame(model.matrix( formule,data=credit2)))
y <- as.numeric(credit[,"Cible"])
y[y == 1] <- 0 # crédits OK
y[y == 2] <- 1 # crédits KO

prim <- prim.box(x, y, peel.alpha=0.5, paste.alpha=0.01,mass.min=0.05, pasting=T, verbose=T,threshold.type=1,threshold=0.3)
summary(prim)
summary(prim,print.box=T)
prim$box

# intersection des boîtes et des données de test
b <- data.frame(row.names(prim$x[[1]]))
b$Cle <- b1$row.names.prim.x..1...
b$row.names.prim.x..1... <- NULL
boxtest = merge(b,test)[,c("Cle","Cible")]
boxtest
table(box1test$Cible)

# intersection des boîtes et des données de test : fonction
maprim=function(i=1)
{b <- data.frame(row.names(prim$x[[i]]));
b$Cle <- b$row.names.prim.x..i...;
b$row.names.prim.x..i... <- NULL;
boxtest = merge(b,test)[,c("Cle","Cible")];
tab <- table(boxtest$Cible);
return(tab)}
#maprim(7)

# aire sous la courbe ROC
library(ROCR)
a <- NULL
b <- NULL
for (i in 1:prim$num.class) {
a <- c(a,c(rep(0,maprim(i)[1]),rep(1,maprim(i)[2])))
b <- c(b,c(rep(i,maprim(i)[1]+maprim(i)[2])))
}
pred <- prediction(b,a,label.ordering=c(1,0))
performance(pred,"auc")@y.values[[1]]

# recherche de l'aire sous la courbe ROC maximale
library(ROCR)
aucprim <- matrix(NA,10,3)
j <- 1
for (n in seq(0.05,0.5,by=0.05)) {
prim <- prim.box(x, y, peel.alpha=n, paste.alpha=0.01,mass.min=0.05,
pasting=T, verbose=T,threshold.type=1,threshold=0.3)
aucprim[j,1] <- n
aucprim[j,2] <- prim$num.class
a <- NULL
b <- NULL
for (i in 1:prim$num.class) {
a <- c(a,c(rep(0,maprim(i)[1]),rep(1,maprim(i)[2])))
b <- c(b,c(rep(i,maprim(i)[1]+maprim(i)[2])))
}
pred <- prediction(b,a,label.ordering=c(1,0))
aucprim[j,3] <- performance(pred,"auc")@y.values[[1]]
j <- j+1
}
colnames(aucprim) <- c("alpha","nb boites","AUC")
aucprim


# mesure de l'AUC en test

formule <- as.formula(paste("y = ~ 0 +",paste(names(credit2[,varquali]),collapse="+")))

# échantillon d'apprentissage
x <- cbind(credit[id,varquanti],data.frame(model.matrix(formule,data=credit2[id,])))
y <- as.numeric(credit[id,"Cible"])
y[y == 1] <- 0 # crédits OK
y[y == 2] <- 1 # crédits KO

# échantillon de test
xt <- cbind(credit[-id,varquanti],data.frame(model.matrix(formule,data=credit2[-id,])))
yt <- as.numeric(credit[-id,"Cible"])
# recodage en 0/1
yt[yt == 1] <- 0 # crédits OK
yt[yt == 2] <- 1 # crédits KO

# modèle PRIM
prim <- prim.box(x, y, peel.alpha=0.5, paste.alpha=0.01,mass.min=0.05, pasting=T, verbose=T,threshold.type=1,threshold=0.3)

# utilisation de la fonction prim.which.box pour mesurer l'AUC
# cette fonction retourne le no de la boîte de x 
head(prim.which.box(x,prim))
tab <- table(prim.which.box(x,prim),y)
prop.table(tab,1)
library(pROC)
auc(y,prim.which.box(x,prim))
auc(yt,prim.which.box(xt,prim))

# recherche de l'aire sous la courbe ROC maximale
aucprim <- matrix(NA,10,4)
j <- 1
for (n in seq(0.05,0.5,by=0.05)) {
prim <- prim.box(x, y, peel.alpha=n, paste.alpha=0.04,mass.min=0.05,
pasting=T, verbose=F,threshold.type=1,threshold=0.3)
aucprim[j,1] <- n
aucprim[j,2] <- prim$num.class
aucprim[j,3] <- auc(y,prim.which.box(x,prim))
aucprim[j,4] <- auc(yt,prim.which.box(xt,prim))
j <- j+1
}
colnames(aucprim) <- c("alpha","nb boites","AUC train","AUC test")
aucprim



# ---------------------------------------------------------------------------------------------------------
# Analyse factorielle comme préalable à PRIM
# ---------------------------------------------------------------------------------------------------------

# chargement du package
library(FactoMineR)
library(prim)

# credit2 = jeu de données avec age, durée et montant de crédit discrétisés
summary(credit2)

# conserver age, durée et montant de crédit sous forme continue
credit3 <- credit2
credit3$Age <- credit$Age
credit3$Duree_credit <- credit$Duree_credit
credit3$Montant_credit <- credit$Montant_credit
summary(credit3)

# conserver age, durée et montant de crédit sous forme continue
credit4 <- credit2
credit4$Taux_effort <- NULL
credit4$Anciennete_domicile <- NULL
credit4$Nb_credits <- NULL
summary(credit4)

# AFDM avec la variable à expliquer en supplémentaire
AFDM <- FAMD (credit2, ncp = 34, graph = TRUE, sup.var = 17)
ACM <- MCA (credit4, ncp = 32, axes = c(1,2), graph = TRUE, quali.sup = 14)
summary(ACM)
library(lattice)
trellis.device(device="pdf", theme=standard.theme("postscript"), file="Ellipse.pdf", color=F)
trellis.device(device="pdf", file="Ellipse.pdf", color=F)
plotellipses(ACM,means=T,level = 0.95,keepvar="Cible",label='quali')
dev.off()
dimdesc(ACM)

# valeurs propres
AFDM$eig
ACM$eig
barplot(AFDM$eig[,2],names=paste("Dim",1:nrow(AFDM$eig)))
barplot(ACM$eig[,2],names=paste("Dim",1:nrow(ACM$eig)))

# individus, coloriés selon "Cible"
plot(AFDM,choix="ind",invisible="quali",habillage=17)
plot(ACM,choix="ind",invisible="var",habillage=14)

# individus coloriés selon leurs modalités sur les 2 premiers axes
plotellipses(ACM,keepvar=1:4)

# variables
plot(AFDM,choix="var")
plot(ACM,choix="var")
plot(ACM,choix="ind",invisible="ind",xlim=c(-1,2),autoLab="yes",cex=0.7)

# résultats pour les individus
AFDM$ind
AFDM$ind$coord[,1:4]
ACM$ind$coord[1:10,1:6]
ACM$var$coord[,1:3]

# ACM avec autres packages

library(MASS)
names(credit4)
ACM2 <- mca(credit4[,-14],nf=32)
ACM2$d
ACM2$d^2 # eigenvalues
ACM2$cs[,1:3]
ACM2$rs[1:10,1:6]
# plan factoriel
plot(ACM2,rows=F)
# ajout au plan factoriel des modalités d'un facteur supplémentaire
predict(ACM2,credit4[14],type="factor") # facteur supplémentaire
points(predict(ACM2,credit4[14],type="factor")[,1:2],cex=10)

# coord(FactoMineR) / (coord(MASS) * sqrt(eigenvalue)) = constante
# la même constante pour tous les axes et toutes les modalités
# coord_var(MASS) = coord_var(FactoMineR) / (p * sqrt(n) * sqrt(eigenvalue)) 
0.15717369/(7.698547e-04*0.4035070)/sqrt(1000)
0.15717369/(16*sqrt(1000)*0.4035070)
0.425655358/(2.453369e-03*0.3429062)
# coord_ind(MASS) = coord_ind(FactoMineR) / (p * sqrt(n)) 
-0.1384053/(16*sqrt(1000))

# number of categories per variable
cats = apply(credit4, 2, function(x) nlevels(as.factor(x)))

#install.packages('ade4')
library(ade4)
ACM4 <- dudi.acm(credit4[,-14], scannf = F, nf = 32)
ACM4$eig
ACM4$co[,1:3]


# PRIM sur les coordonnées factorielles

# échantillon d'apprentissage
x <- ACM$ind$coord[id,1:8]
y <- as.numeric(credit[id,"Cible"])
y[y == 1] <- 0 # crédits OK
y[y == 2] <- 1 # crédits KO

# échantillon de test
xt <- ACM$ind$coord[-id,1:8]
yt <- as.numeric(credit[-id,"Cible"])
# recodage en 0/1
yt[yt == 1] <- 0 # crédits OK
yt[yt == 2] <- 1 # crédits KO

library(prim)
prim <- prim.box(x, y, peel.alpha=0.05, paste.alpha=0.04,mass.min=0.05, pasting=T, verbose=T,threshold.type=1,threshold=0.3)
summary(prim,print.box=T)

# recherche de l'aire sous la courbe ROC maximale
library(pROC)
aucprim <- matrix(NA,10,4)
j <- 1
for (n in seq(0.05,0.5,by=0.05)) {
prim <- prim.box(x, y, peel.alpha=n, paste.alpha=0.06,mass.min=0.05,
pasting=T, verbose=F,threshold.type=1,threshold=0.3)
aucprim[j,1] <- n
aucprim[j,2] <- prim$num.class
#aucprim[j,3] <- auc(y,prim.which.box(x,prim))
#aucprim[j,4] <- auc(yt,prim.which.box(xt,prim))
aucprim[j,3] <- auc(y,prim$y.fun [prim.which.box(x,prim)])
aucprim[j,4] <- auc(yt,prim$y.fun [prim.which.box(xt,prim)])
j <- j+1
}
colnames(aucprim) <- c("alpha","nb boites","AUC train","AUC test")
aucprim


# calcul de l'AUC en test pour chaque combinaison (alpha,beta)
library(pROC)
f <- function(i,j)
{
prim <- prim.box(x, y, peel.alpha=i, paste.alpha=j,mass.min=0.05,
pasting=T, verbose=F,threshold.type=1,threshold=0.3)
auc(yt,prim.which.box(xt,prim))
}
f(0.5,0.03)

# utilisation de la fonction f précédente pour trouver
# le nb d'axes factoriels, et les valeurs de alpha et beta
# maximisant l'AUC en test
for (n in seq(2,12,by=1))
{
i <- seq(0.05,0.5,by=0.05)
j <- seq(0.01,0.1,by=0.01)
x <- ACM$ind$coord[id,1:n]
xt <- ACM$ind$coord[-id,1:n]
k <- outer(i,j,Vectorize(f))
cat("\n","nb facteurs = ",n," AUC max = ",max(k)," ")
}

for (n in seq(2,12,by=1))
{
i <- seq(0.05,0.5,by=0.05)
j <- seq(0.01,0.1,by=0.01)
x <- ACM$ind$coord[id,1:n]
xt <- ACM$ind$coord[-id,1:n]
k <- outer(i,j,Vectorize(f))
alpha1 <- i[which(k==max(k),arr.ind=TRUE)[1,1]]
alpha2 <- j[which(k==max(k),arr.ind=TRUE)[1,2]]
cat("\n","nb facteurs = ",n," AUC test max = ",max(k)," alpha peel = ",
alpha1," alpha paste = ",alpha2)
}

for (n in seq(2,12,by=1))
{
i <- seq(0.05,0.5,by=0.05)
j <- seq(0.01,0.1,by=0.01)
x <- AFDM$ind$coord[id,1:n]
xt <- AFDM$ind$coord[-id,1:n]
k <- outer(i,j,Vectorize(f))
alpha1 <- i[which(k==max(k),arr.ind=TRUE)[1,1]]
alpha2 <- j[which(k==max(k),arr.ind=TRUE)[1,2]]
cat("\n","nb facteurs = ",n," AUC test max = ",max(k)," alpha peel = ",
alpha1," alpha paste = ",alpha2)
}


# fonction implémentant les forêts aléatoires sur un modèle PRIM

RForestPRIM = function(y, x, yt, xt, nb_var, alpha, beta, nsimul)
{

nobs <- length(y)
nvar <- ncol(x)
predic   <- matrix(NA,length(yt),nsimul)
auc      <- matrix(0,nsimul,1)

for(i in (1:nsimul))
{

# tirage aléatoire simple des variables
sv <- sort(sample(nvar, nb_var, replace=F))
# tirage aléatoire avec remise de l'échantillon d'apprentissage
si <- sort(sample(nobs,nobs,replace=T))
alpha1 <- runif(1,0.3,0.5)
beta1 <- runif(1,0.03,0.07)

# formule avec l'ensemble des prédicteurs
prim <- prim.box(x[si,sv], y[si], peel.alpha=alpha1, paste.alpha=beta1,mass.min=0.05,
pasting=T, verbose=F,threshold.type=1,threshold=0.3)

# application du modèle logit à l'échantillon de validation
predic[,i] <- prim$y.fun [pmax(1,prim.which.box(xt[,sv],prim))]
if (i == 1) {moypred <- predic[,i]} else {moypred <- apply(predic[,1:i],1,mean)}
auc[i] <- auc(yt,moypred)

} # fin de la boucle

# résultats en retour
rf <- list(auc)

} # fin de la fonction

niter <- 1000

# échantillon d'apprentissage
x <- ACM$ind$coord[id,1:18]
y <- as.numeric(credit[id,"Cible"])
y[y == 1] <- 0 # crédits OK
y[y == 2] <- 1 # crédits KO

# échantillon de test
xt <- ACM$ind$coord[-id,1:18]
yt <- as.numeric(credit[-id,"Cible"])
# recodage en 0/1
yt[yt == 1] <- 0 # crédits OK
yt[yt == 2] <- 1 # crédits KO

# appel de la fonction
rf <- RForestPRIM(y, x, yt, xt, nb_var=4, alpha=0.5, beta=0.06, nsimul=niter)
plot(rf[[1]])



# ---------------------------------------------------------------------------------------------------------
# CHAPITRE 11
# Forêts aléatoires
# ---------------------------------------------------------------------------------------------------------

library(randomForest)
#credit$Cible <- factor(credit$Cible)

# remarque : la variable à expliquer doit avoir le type "factor" pour un classement
# ntree = 500 arbres agrégés
# mtry = 4 variables sélectionnées pour la scission de chaque noeud
# importance = TRUE calcule l'importance des variables
# replace = TRUE (valeur par défaut) pour tirage avec remise des individus
# keep.forest = TRUE pour conserver la forêt dans l'objet en sortie et pouvoir
#                    l'appliquer à un autre échantillon
# nodesize = effectif minimum de chaque feuille (par défaut : 1 en classement
#                    et 5 en régression)

set.seed(235)
rf <- randomForest(Cible ~ ., data=credit[id,vars], importance=TRUE,
proximity=TRUE,
ntree=500, mtry=3, replace=T, keep.forest=T,nodesize=5,ytest=test[,"Cible"],xtest=test[,1:19])
rf

set.seed(235)
ptm <- proc.time()
rf <- randomForest(Cible ~ ., data=credit[id,vars], importance=TRUE,
proximity=TRUE,
ntree=500, mtry=3, replace=T, keep.forest=T,nodesize=5,ytest=test[,"Cible"],xtest=test[,1:19],cutoff=c(0.5,0.5))
proc.time() - ptm
rf

stump <- randomForest(Cible ~ ., data=credit[id,vars], importance=TRUE,
ntree=500, mtry=1, replace=T, keep.forest=T,maxnodes=2,ytest=test[,"Cible"],xtest=test[,1:19])

# erreur
head(rf$votes) # votes pour chaque individu OOB des arbres de l'échantillon d'apprentissage
head(predict(rf,type='prob')) # commande identique à la précédente
mean(rf$oob.times) # nb de fois où chaque individu est out of bag
head(rf$oob.times)
n <- 100000000000000
rf$ntree*(1-(1/n))^n
head(rf$err.rate)
head(rf$test$err.rate)
tail(rf$err.rate)
tail(rf$test$err.rate)

# taux d'erreur en test
testrf <- predict(rf,test,type='response')
head(testrf)
head(test$Cible)
table(test$Cible,testrf)
sum(testrf != test$Cible) / nrow(test) # 0.2275281
# taux d'erreur en apprentissage
trainrf <- predict(rf,credit[id,],type='response')
table(credit[id,"Cible"],trainrf)
err_train <- sum(trainrf != credit[id,"Cible"]) / nrow(credit[id,])
# taux d'erreur en apprentissage OOB
#rfOOB <- ifelse(rf$votes[,2]>=0.5,1,0)
table(credit[id,"Cible"],rfOOB)
err_OOB <- sum(rf$pred != credit[id,"Cible"]) / nrow(credit[id,])
# estimateur .632 
0.368*err_train + 0.632*err_OOB # 0.1654286
0.368*0.02018634 + 0.632*0.25 # 0.1654286
# estimateur .632+
p <- table(credit[id,"Cible"])[2]/nrow(credit[id,])
p
q <- table(trainrf)[2] / nrow(credit[id,])
q
gamma <- p*(1-q) + (1-p)*q
gamma
txoverfit <- (err_OOB-err_train)/(gamma-err_train)
txoverfit
w <- 0.632/(1-0.368*txoverfit)
w
(1-w)*err_train + w*err_OOB # 0.2054516

# évolution du taux d'erreur (calculé sur individus OOB de l'échantillon d'apprentissage)
plot(rf$err.rate[1:rf$ntree,1],type='l',ylim=c(.2,.4),
xlab="nombre d'itérations",ylab='erreur')
lines(rf$test$err.rate[1:rf$ntree,1],type='l',lwd=2,col='red')
#abline(h=err_cart,col='green')
plot(rf$err.rate[,1],type='l',ylim=c(.2,.4),xlab="nombre d'itérations",ylab='erreur')
lines(rf$test$err.rate[,1],type='l',lwd=2,col='red')

# taux d'erreur minimum en apprentissage
(min.err <- min(rf$err.rate[,"OOB"]))
# nb d'itérations atteignant l'erreur minimale en apprentissage
(min.err.idx <- which(rf$err.rate[,"OOB"]== min.err))

# taux d'erreur minimum en test
(min.err <- min(rf$test$err.rate[,"Test"]))
# nb d'itérations atteignant l'erreur minimale en validation
(min.err.idx <- which(rf$test$err.rate[,"Test"]== min.err))

# prédiction
library(ROCR)
pred.rf <- predict(rf,test,type='prob')[,2]
pred <- prediction(pred.rf,test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# graphique de proximité
MDSplot(rf, credit[id,vars]$Cible)
MDSplot(rf, credit[id,vars]$Cible,pch=as.numeric(credit[id,vars]$Cible))

# graphique de l'importance des variables
importance(rf,type=1,scale=T)
importance(rf,type=1)
importance(rf,type=2)
varImpPlot(rf)
sd(as.vector(importance(rf,type=1,scale=T)))
sd(as.vector(importance(rf,type=2,scale=T)))
rf.imp <- importance(rf,type=1)[order(importance(rf,type=1), decreasing=TRUE),]
# la ligne précédente équivaut aux deux suivantes
#imp.tmp <- importance(rf, type = 1)
#rf.imp <- imp.tmp[order(imp.tmp, decreasing=TRUE),]
rf.imp # variables par ordre décroissant d'importance

par(mar = c(8, 4, 4, 0))
barplot(rf.imp, col = gray(0:nrow(rf$importance)/nrow(rf$importance)),
        ylab='Importance', ylim=c(0,30),cex.names = 0.8, las=3)
title(main = list("Importance des variables", font = 1, cex = 1.2))

# importance moyenne des variables
nvar <- ncol(credit[,vars])-1
vimp <- rep(0,nvar)
nsimul <- 30
for(i in 1:nsimul){
rf <- randomForest(Cible ~ ., data=credit[id,vars], importance=TRUE,
ntree=500, mtry=3, replace=T, keep.forest=T,nodesize=5)
vip <- importance(rf,type=1)
vimp <- vimp + vip[order(rownames(vip))] / nsimul
rm(rf)
}
#vimp
#vip
a1 <- sort(rownames(vip))
a2 <- order(vimp, decreasing = TRUE)
dotchart(vimp[a2][nvar:1], a1[a2][nvar:1], main="randomForest : Importance moyenne des variables")

# autre graphique de l'importance moyenne des variables
par(mar = c(8, 4, 4, 0))
barplot(vimp[a2], col = gray(0:nvar/nvar),names.arg = a1[a2],
        ylab='Importance', ylim=c(0,30), cex.names = 0.8, las=3)
title(main = list("randomForest : Importance moyenne des variables", font = 1, cex = 1.2))

# détermination du nb optimal de variables sélectionnées
set.seed(235)
mtry <- tuneRF(credit[id,-c(20,21)], credit[id,"Cible"], mtryStart = 1, 
ntreeTry=500,stepFactor=2, improve = 0.001)
mtry
mtry[,2] # taux d'erreur OOB pour chaque valeur testée de mtry
best.m <- mtry[mtry[, 2] == min(mtry[, 2]), 1]
best.m <- mtry[which.min(mtry[,2]),1]
best.m
# LIMITES :
# - cette détermination n'est faite qu'au vu de l'erreur OOB et non en test
# - deux calculs ne donnent pas la même valeur de mtry !

# simulations sur 30 random forests et calcul erreur moyenne
nsimul <- 30
rft  <- matrix(NA,3,nsimul+1)
for(i in 1:nsimul)
{
rf <- randomForest(Cible ~ ., data=credit[id,vars], importance=TRUE,
ntree=500, mtry=5, replace=T, keep.forest=T,nodesize=5,ytest=test[,"Cible"],
xtest=test[,2:20])
rft[1,i] <- min(data.frame(rf$err.rate)["OOB"])
rft[2,i] <- min(data.frame(rf$test$err.rate)["Test"])
rft[3,i] <- rf$test$err.rate[500,"Test"]
} # fin de la boucle
rft[,nsimul+1] <- apply(rft[,1:nsimul],1,mean) # moyenne des erreurs
rft

# simulations sur 100 random forests et calcul AUC moyenne
library(ROCR)
set.seed(235)
nsimul <- 100
nvarmin <- 1
nvarmax <- 6
auc  <- matrix(NA,nvarmax-nvarmin+1,2)
for(nvar in nvarmin:nvarmax)
{
auc[nvar-nvarmin+1,1] <- nvar
rft  <- matrix(NA,nrow(test),nsimul+1)
for(i in 1:nsimul)
{
rf <- randomForest(Cible ~ ., data=credit[id,vars], importance=F,
ntree=500, mtry=nvar, replace=T, keep.forest=T, nodesize=5)
# nodesize=5
# maxnodes=2
rft[,i] <- predict(rf,test,type='prob')[,2]
} # fin de la boucle intérieure
rft[,nsimul+1] <- apply(rft[,1:nsimul],1,mean) # moyenne des prédictions
pred <- prediction(rft[,nsimul+1],test$Cible,label.ordering=c(0,1))
auc[nvar-nvarmin+1,2] <- performance(pred,"auc")@y.values[[1]]
} # fin de la boucle extérieure
colnames(auc) <- c("nb variables","AUC test")
auc


# RF sans bootstrap (juste la randomisation des variables)
set.seed(235)
rf <- randomForest(Cible ~ ., data=credit[id,vars], importance=TRUE,
ntree=500, mtry=2, replace=F, sampsize=length(id), keep.forest=T,nodesize=5,ytest=test[,"Cible"],
xtest=test[,1:19])
test$rf <- predict(rf,test,type='prob')[,2]
pred <- prediction(test$rf,test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]


# RF en sélectionnant toutes les variables (= bagging : juste le bootstrap)
set.seed(235)
rf <- randomForest(Cible ~ ., data=credit[id,vars], importance=TRUE,
ntree=500, mtry=dim(credit)[2]-2, replace=T, keep.forest=T,nodesize=5,ytest=test[,"Cible"],
xtest=test[,1:19])
test$rf <- predict(rf,test,type='prob')[,2]
pred <- prediction(test$rf,test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]


# imputation de valeurs manquantes
summary(train)
logit <- glm(Cible~Comptes+Historique_credit+Duree_credit+Age+Epargne+Garanties+Autres_credits,data=train,family=binomial(link = "logit"))
summary(logit)
library(missForest)
train.mis <- prodNA(train, noNA = 0.1)
summary(train.mis)
logit.mis <- glm(Cible~Comptes+Historique_credit+Duree_credit+Age+Epargne+Garanties+Autres_credits,data=train.mis,family=binomial(link = "logit"))
summary(logit.mis)
train.imp <- missForest(train.mis)
train.imp$OOBerror
train.imp$error
logit.imp <- glm(Cible~Comptes+Historique_credit+Duree_credit+Age+Epargne+Garanties+Autres_credits,data=train.imp$ximp,family=binomial(link = "logit"))
summary(logit.imp)
# parallélisation
library(doSNOW)
cl <- makeCluster(4)
registerDoSNOW(cl)
ptm <- proc.time()
#train.imp <- missForest(train.mis,xtrue=train,parallelize='no')
train.imp <- missForest(train.mis,xtrue=train,parallelize='forests')
#train.imp <- missForest(train.mis,xtrue=train,parallelize='variables')
proc.time() - ptm
# fin parallélisation
stopCluster(cl)



# ---------------------------------------------------------------------------------------------------------
# Extra-trees
# ---------------------------------------------------------------------------------------------------------

install.packages('extraTrees')
library(extraTrees)
summary(credit)

# le vecteur de prédicteurs x de extraTrees doit être numérique => toutes
# les variables qui sont des facteurs sont remplacées par les indicatrices
# de leurs modalités (sauf la modalité de référence)
x <- model.matrix( ~ . -1 ,data=credit[id,!names(credit)%in%c("Cible","Cle")])
y <- credit[id,"Cible"]
xt <- model.matrix( ~ . -1 ,data=credit[-id,!names(credit)%in%c("Cible","Cle")])

set.seed(235)
ptm <- proc.time()
et <- extraTrees(x, y, ntree=1000, mtry=5, numRandomCuts=1, nodesize=5, numThreads=1)
proc.time() - ptm
et

# prédiction
#library(ROCR)
pred.et <- predict(et,xt,probability=T)[,2]
pred <- prediction(pred.et,test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# parallélisation
install.packages('microbenchmark')
library(microbenchmark)
microbenchmark(extraTrees(x, y, ntree=500, mtry=5, numRandomCuts=1, nodesize=5, numThreads=1), times = 30)
microbenchmark(extraTrees(x, y, ntree=500, mtry=5, numRandomCuts=1, nodesize=5, numThreads=4), times = 10)



# ---------------------------------------------------------------------------------------------------------
# Forêts aléatoires parallélisées
# ---------------------------------------------------------------------------------------------------------

library(randomForest)
summary(credit)
names(credit)

# version sans foreach
ptm <- proc.time()
rf <- randomForest(Cible ~ ., data=credit[id,vars], importance=TRUE,
ntree=10000, mtry=3, replace=T, keep.forest=T,nodesize=5,ytest=test[,"Cible"],
xtest=test[,1:19])
proc.time() - ptm
rf


# appel foreach sans parallélisation
library(foreach)
#rf <- foreach(ntree=rep(250, 4), .combine=combine) %do% randomForest(x, y, ntree=ntree)

ptm <- proc.time()
rf <- foreach(ntree=rep(2500, 4), .combine=combine) %do% 
randomForest(Cible ~ ., data=credit[id,vars], importance=TRUE,
ntree=ntree, mtry=3, replace=T, keep.forest=T,nodesize=5,ytest=test[,"Cible"],
xtest=test[,1:19])
proc.time() - ptm


# appel foreach avec parallélisation
library(parallel)
detectCores()
detectCores(logical=FALSE)
library(doSNOW)
cl <- makeCluster(2) #change the 2 to your number of physical CPU cores
registerDoSNOW(cl)

#rf <- foreach(ntree=rep(250, 4), .combine=combine, .packages='randomForest') %dopar% randomForest(x, y, ntree=ntree)
ptm <- proc.time()
rf <- foreach(ntree=rep(2500, 4), .combine=combine, .packages='randomForest') %dopar% 
randomForest(Cible ~ ., data=credit[id,vars], importance=TRUE,
ntree=ntree, mtry=3, replace=T, keep.forest=T,nodesize=5,ytest=test[,"Cible"],
xtest=test[,1:19])
proc.time() - ptm

# fin parallélisation
stopCluster(cl)



# ---------------------------------------------------------------------------------------------------------
# CHAPITRE 12
# Bagging
# ---------------------------------------------------------------------------------------------------------

#install.packages("ipred")
library(ipred)
library(rpart)
library(ROCR)
#library(adabag)

# bagging avec arbres de profondeur maximale
set.seed(235)
bag <- bagging(Cible ~ ., data=credit[id,vars],nbagg=200,coob=TRUE,
control= rpart.control(cp=0))
bag

# prédiction
test$bag <- predict(bag,type="prob",test)
pred <- prediction(test$bag[,2],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# mêmes analyses avec taille minimale des feuilles
set.seed(235)
bag1 <- bagging(Cible ~ ., data=credit[id,vars],nbagg=200,coob=TRUE,
control= rpart.control(minbucket=5))
bag1

# agrégation par moyenne des probabilités
test$bag1 <- predict(bag1, test, type="prob", aggregation="average")
head(test$bag1)
pred <- prediction(test$bag1[,2],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# agrégation par vote à la majorité
test$bag1 <- predict(bag1, test, type="prob", aggregation="majority")
head(test$bag1)
pred <- prediction(test$bag1[,2],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# bagging de stumps
set.seed(235)
stump <- bagging(Cible ~ ., data=credit[id,vars],nbagg=100,coob=TRUE,
control= rpart.control(maxdepth=1,cp=-1,minsplit=0))
bag1 <- stump

# rappel : AUC de stumps
cart <- rpart(Cible ~ . ,data = credit[id,vars],method="class",parms=list(split="gini"),
control=list(maxdepth=1,cp=-1,minsplit=0))
test$stump <- predict(cart, test, type="prob")
pred <- prediction(test$stump[,2],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# bagging d'arbres peu discriminants
library(rpart)
cart <- rpart(Cible ~ Anciennete_domicile+Telephone+Nb_pers_charge ,data = credit[id,vars],method="class",parms=list(split="gini"),cp=-1,maxdepth=1)
cart

plot(cart,branch=.2, uniform=T, compress=T, margin=.1)
text(cart, fancy=T,use.n=T,pretty=0,all=T,cex=.5)
predc <- predict(cart,type="prob",test)[,2]
pred <- prediction(predc,test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

set.seed(235)
bag1 <- bagging(Cible ~ Anciennete_domicile+Telephone+Nb_pers_charge, data=credit[id,vars],nbagg=200,coob=TRUE,
control= rpart.control(maxdepth=1,cp=-1))
test$bag1 <- predict(bag1, test, type="prob", aggregation="average")
pred <- prediction(test$bag1[,2],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.4924933
test$bag1 <- predict(bag1, test, type="prob", aggregation="majority")
pred <- prediction(test$bag1[,2],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.5

# utilisation du package randomForest
library(randomForest)
set.seed(235)
bag2 <- randomForest(Cible ~ ., data=credit[id,vars], importance=TRUE,
ntree=500, mtry=length(credit[id,vars])-1, replace=T, keep.forest=T,nodesize=5)
test$bag2 <- predict(bag2,test,type='prob')[,2]
pred <- prediction(test$bag2,test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]



# ---------------------------------------------------------------------------------------------------------
# CHAPITRE 13
# Forêts aléatoires logistiques
# ---------------------------------------------------------------------------------------------------------

library(ROCR)

# Attention : regrouper la modalité "Autres" de "Objet_credit" qui n'a que 12 observations
# La rareté de cette modalité fait qu'elle peut se trouver absente d'un échantillon
# d'apprentissage alors qu'elle sera présente dans l'échantillon de test

table(credit$Objet_credit)
credit$Objet_credit[credit$Objet_credit == "A410"] <- "A48"
table(credit$Objet_credit)
train  <- credit[id,]
valid  <- credit[-id,]
summary(train)

library(car)
credit2$Objet_credit <- recode(credit$Objet_credit,"'A40'='Voiture neuve';
'A41'='Voiture occasion';c('A42','A44','A45')='Intérieur';
'A43'='Vidéo - HIFI';c('A46','A48','A410')='Etudes';
'A47'='Vacances';'A49'='Business';else='Autres'")
train  <- credit2[id,]
valid  <- credit2[-id,]
summary(train)


# fonction implémentant les forêts aléatoires sur un modèle logit

RForestLog = function(apprent, validat, varY, varX, nb_var, crit="bic", nsimul)
{

nobs <- nrow(apprent)
y <- apprent[,varY]
formule <- as.formula(paste("y ~ ",paste(names(apprent[,varX]),collapse="+")))
sature   <- glm(formule,data=apprent[,c(varX,varY)], family=binomial(link = "logit"))
#cible <- apprent[,varY]
#sature   <- glm(Cible~.,data=apprent[,c(varX,varY)], family=binomial(link = "logit"))
coef     <- matrix(0,length(coef(sature)),nsimul+2)
rownames(coef) <- names(coef(sature))
predic   <- matrix(NA,nrow(validat),nsimul)
auc      <- matrix(0,nsimul,1)

for(i in (1:nsimul))
{

# tirage aléatoire simple des variables
size <- length(varX)
s <- sort(sample(size, nb_var, replace=F))
predicteurs <- varX[s]

# tirage aléatoire équiprobable avec remise des individus
s <- sample(nobs,nobs,replace=T)

# formule avec l'ensemble des prédicteurs sélectionnés aléatoirement
if (nb_var > 1) {
formule <- as.formula(paste("y ~ ",paste(names(apprent[,predicteurs]),collapse="+")))
} else {
formule <- as.formula(paste("y ~ ",predicteurs))
}
cible <- apprent[s,varY]
# modèle minimum
logit <- glm(cible ~ 1,data=apprent[s,], family=binomial(link = "logit"))

# sélection pas à pas sur la base du critère BIC ou AIC
if (crit == "bic") {
selection <- step(logit,direction="forward",trace=F,k=log(nobs),scope=list(upper=formule))
} else {
selection <- step(logit,direction="forward",trace=F,k=2,scope=list(upper=formule))
}

# application du modèle logit à l'échantillon de validation
predic[,i] <- predict(selection, newdata=validat, type="response")
# moyenne des prédictions des i premiers modèles agrégés
#if (i == 1) {moypred <- predic[,i]} else {moypred <- apply(predic[,1:i],1,mean)}
if (i == 1) {moypred <- predic[,i]} else {moypred <- rowMeans(predic[,1:i])}
# calcul de l’AUC des i premiers modèles agrégés
cible <- validat[,varY]
pred <- prediction(moypred,cible,label.ordering=c(0,1))
auc[i] <- performance(pred,"auc")@y.values[[1]]
# stockage des coefficients
index <- match(names(coef(selection)),names(coef(sature)))
coef[index,i] <- coef(selection)

} # fin de la boucle

coef[,nsimul+1] <- rowMeans(coef[,1:nsimul]) # moyenne des coefficients
coef[,nsimul+2] <- rowSums(coef[,1:nsimul]!=0) # nb de coefficients non nuls par modalité

# résultats en retour
rf <- list(auc,coef)

} # fin de la fonction

niter <- 300
varx <- c(varquali, varquanti)
#npred <- -grep('(Type_emploi|Telephone|Nb_pers_charge|Anciennete_domicile|Taux_effort|Nb_credits)', varx)
#varx <- varx[npred]
varx
ptm <- proc.time()
rf <- RForestLog(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=3, crit="aic", nsimul=niter)
proc.time() - ptm

# AUC
rf[[1]]
plot(rf[[1]])
rf[[2]]
# coefficients moyens
rf[[2]][,niter+1]
# nb de sélections de chaque variable
rf[[2]][,niter+2]
rf[[2]]

# fonction lissant les AUC et les coefficients de forêts aléatoires
# en effectuant plusieurs simulations

Ensemble = function(apprent, validat, varY, varX, nb_var, crit="bic", nsimul, nb)
{

y <- apprent[,varY]
formule <- as.formula(paste("y ~ ",paste(names(apprent[,varX]),collapse="+")))
sature   <- glm(formule,data=apprent[,c(varX,varY)], family=binomial(link = "logit"))
coe <- matrix(0,length(coef(sature)),nb+1)
rownames(coe) <- names(coef(sature))
occ <- matrix(0,length(coef(sature)),nb+1)
rownames(occ) <- names(coef(sature))
auc <- matrix(0,nsimul,nb+1)

for(i in (1:nb))
{
rf <- RForestLog(apprent, validat, varY, varX, nb_var, crit, nsimul)
auc[,i] <- rf[[1]]
coe[,i] <- rf[[2]][,nsimul+1]
occ[,i] <- rf[[2]][,nsimul+2]
}

auc[,nb+1] <- apply(auc[,1:nb],1,mean) # moyenne des AUC
coe[,nb+1] <- apply(coe[,1:nb],1,mean) # moyenne des coefficients
occ[,nb+1] <- apply(occ[,1:nb],1,mean) # moyenne des nb d'occurrences

# résultats en retour
ens <- list(auc,coe,occ)

} # fin de la fonction


niter <- 300
n     <- 30
varx <- c(varquali, varquanti)
varx

ptm <- proc.time()
ens <- Ensemble(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=2, crit="aic", nsimul=niter, nb=n)
proc.time() - ptm

ens[[1]][niter,n+1]
plot(ens[[1]][,n+1])



# ---------------------------------------------------------------------------------------------------------
# Random Forest sur Logit : version parallélisée
# ---------------------------------------------------------------------------------------------------------

library(foreach)
library(iterators)

# fonction implémentant les forêts aléatoires sur un modèle logit

pRForestLog = function(apprent, validat, varY, varX, nb_var, crit="bic", nsimul)
{

nobs <- nrow(apprent)
y <- apprent[,varY]
formule <- as.formula(paste("y ~ ",paste(names(apprent[,varX]),collapse="+")))
sature   <- glm(formule,data=apprent[,c(varX,varY)], family=binomial(link = "logit"))
coef     <- matrix(0,length(coef(sature)),nsimul+2)
rownames(coef) <- names(coef(sature))
predic   <- matrix(0,nrow(validat),1)
auc      <- matrix(0,nsimul,1)

rf <- foreach(icount(nsimul), .combine='cbind') %dopar% {

# tirage aléatoire simple des variables
size <- length(varX)
s <- sort(sample(size, nb_var, replace=F))
predicteurs <- varX[s]

# tirage aléatoire équiprobable avec remise des individus
s <- sample(nobs,nobs,replace=T)

# formule avec l'ensemble des prédicteurs sélectionnés aléatoirement
if (nb_var > 1) {
formule <- as.formula(paste("y ~ ",paste(names(apprent[,predicteurs]),collapse="+")))
} else {
formule <- as.formula(paste("y ~ ",predicteurs))
}
cible <- apprent[s,varY]
# modèle minimum
logit <- glm(cible ~ 1,data=apprent[s,], family=binomial(link = "logit"))

# sélection pas à pas sur la base du critère BIC ou AIC
if (crit == "bic") {
selection <- step(logit,direction="forward",trace=F,k=log(nobs),scope=list(upper=formule))
} else {
selection <- step(logit,direction="forward",trace=F,k=2,scope=list(upper=formule))
}

# application du modèle logit à l'échantillon de validation
pred <- predict(selection, newdata=validat, type="response")
return(pred)

} # fin du foreach

forest <- t(apply(rf, 1, cumsum))
aires <- as.vector(apply(forest,2,function(x) auc(valid$Cible,x)))
return(aires)

} # fin de la fonction


niter <- 300
varx <- c(varquali, varquanti)
varx
ptm <- proc.time()
rfo <- pRForestLog(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=3, crit="aic", nsimul=niter)
proc.time() - ptm

plot(rfo)
plot(names(rfo[,niter]),rfo[,niter])


# parallélisation sur plusieurs coeurs

library(parallel)
detectCores()
detectCores(logical=FALSE)

#install.packages('doSNOW')
library(doSNOW)
library(foreach)

cl <- makeCluster(4) #change the 2 to your number of physical CPU cores
registerDoSNOW(cl)

ptm <- proc.time()
rfo <- pRForestLog(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=3, crit="aic", nsimul=niter)
proc.time() - ptm

stopCluster(cl)



# ---------------------------------------------------------------------------------------------------------
# Random Forest sur Logit : autre version parallélisée
# ---------------------------------------------------------------------------------------------------------

library(parallel)
detectCores()
detectCores(logical=FALSE)

#install.packages('doSNOW')
library(doSNOW)
library(foreach)

cl <- makeCluster(4)
registerDoSNOW(cl)
getDoParWorkers()

pEnsemble = function(apprent, validat, varY, varX, nb_var, crit="bic", nsimul, nb)
{

y <- apprent[,varY]
formule <- as.formula(paste("y ~ ",paste(names(apprent[,varX]),collapse="+")))
sature   <- glm(formule,data=apprent[,c(varX,varY)], family=binomial(link = "logit"))
coe <- matrix(0,length(coef(sature)),nb+1)
rownames(coe) <- names(coef(sature))
occ <- matrix(0,length(coef(sature)),nb+1)
rownames(occ) <- names(coef(sature))
auc <- matrix(0,nsimul,nb+1)

rf <- foreach(i=1:nb,.export=c("RForestLog"),.packages='ROCR') %dopar%
{ RForestLog(apprent, validat, varY, varX, nb_var, crit, nsimul) }

roc <- function(x) {rf[[x]][[1]]}
auc <- Vectorize(roc)(1:nb)
aucm <- apply(auc[,1:nb],1,mean)

cf <- function(x) {rf[[x]][[2]][,nsimul+1]}
coe <- Vectorize(cf)(1:nb)
coem <- apply(coe[,1:nb],1,mean)

oc <- function(x) {rf[[x]][[2]][,nsimul+2]}
occ <- Vectorize(oc)(1:nb)
occm <- apply(occ[,1:nb],1,mean)

# résultats en retour
ens <- list(aucm,coem,occm)

} # fin de la fonction

niter <- 300
n     <- 30
varx <- c(varquali, varquanti)

ptm <- proc.time()
p.ens <- pEnsemble(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=2, crit="aic", nsimul=niter, nb=n)
proc.time() - ptm

p.ens[[1]][niter]
plot(p.ens[[1]])
p.ens[[2]]
p.ens[[3]]

# fin parallélisation
stopCluster(cl)

# grille de score
score <- rf[[2]]
sature   <- glm(Cible~.,data=train[,c(varx,"Cible")], family=binomial(link = "logit"))
sature$xlevels
names(unlist(sature$xlevels))
VARIABLE=c("",gsub("[0-9]", "", names(unlist(sature$xlevels))))
#VARIABLE = c(varquali, varquanti)
VARIABLE
MODALITE=c("",as.character(unlist(sature$xlevels)))
MODALITE
names=data.frame(VARIABLE,MODALITE,NOMVAR=c("(Intercept)",paste(VARIABLE,MODALITE,sep="")[-1]))
#regression=data.frame(NOMVAR=names(rf[[2]][,niter+1]),COEF=as.numeric(rf[[2]][,niter+1]))
#regression=data.frame(NOMVAR=names(ens[[2]][,nb+1]),COEF=as.numeric(ens[[2]][,nb+1]))
regression=data.frame(NOMVAR=names(score),COEF=as.numeric(score))
regression
param = merge(names,regression,all.x=TRUE)[-1]
param$COEF[is.na(param$COEF)] <- 0 
mini=aggregate(data.frame(min = param$COEF), by = list(VARIABLE = param$VARIABLE), min)
maxi=aggregate(data.frame(max = param$COEF), by = list(VARIABLE = param$VARIABLE), max)
total=merge(mini,maxi)
total$diff = total$max - total$min
poids_total = sum(total$diff)
grille = merge(param,mini,all.x=TRUE)
grille$delta = grille$COEF - grille$min
grille$POIDS = round((100*grille$delta) / poids_total)
grille[which(VARIABLE!=""),c("VARIABLE","MODALITE","POIDS")]



# ---------------------------------------------------------------------------------------------------------
# Random Forest sur Logit : tests divers
# ---------------------------------------------------------------------------------------------------------

niter <- 300
n     <- 30
varx
library(ROCR)
bic1 <- Ensemble(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=1, crit="bic", nsimul=niter, nb=n)
aic1 <- Ensemble(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=1, crit="aic", nsimul=niter, nb=n)
bic2 <- Ensemble(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=2, crit="bic", nsimul=niter, nb=n)
aic2 <- Ensemble(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=2, crit="aic", nsimul=niter, nb=n)
bic3 <- Ensemble(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=3, crit="bic", nsimul=niter, nb=n)
aic3 <- Ensemble(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=3, crit="aic", nsimul=niter, nb=n)
bic4 <- Ensemble(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=4, crit="bic", nsimul=niter, nb=n)
aic4 <- Ensemble(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=4, crit="aic", nsimul=niter, nb=n)
bic5 <- Ensemble(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=5, crit="bic", nsimul=niter, nb=n)
aic5 <- Ensemble(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=5, crit="aic", nsimul=niter, nb=n)
bic8 <- Ensemble(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=8, crit="bic", nsimul=niter, nb=n)
aic8 <- Ensemble(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=8, crit="aic", nsimul=niter, nb=n)
bic13 <- Ensemble(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=13, crit="bic", nsimul=niter, nb=n)
bic19 <- Ensemble(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=19, crit="bic", nsimul=niter, nb=n)
aic19 <- Ensemble(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=19, crit="aic", nsimul=niter, nb=n)

library(ggplot2)
trait <- 1.25
ggplot() + 
  geom_line(aes(seq(1,300), aic3[[1]][,n+1]), colour='black', size=trait, linetype=1, data.frame(seq(1,300),aic3[[1]][,n+1])) +
  geom_line(aes(seq(1,300), aic4[[1]][,n+1]), colour='black', size=trait, linetype=3, data.frame(seq(1,300),aic4[[1]][,n+1])) +
  geom_line(aes(seq(1,300), bic3[[1]][,n+1]), colour='grey60', size=trait, linetype=1, data.frame(seq(1,300),bic3[[1]][,n+1])) +  
  geom_line(aes(seq(1,300), bic4[[1]][,n+1]), colour='grey60', size=trait, linetype=5, data.frame(seq(1,300),bic4[[1]][,n+1])) +
  geom_line(aes(seq(1,300), bic8[[1]][,n+1]), colour='grey60', size=trait, linetype=4, data.frame(seq(1,300),bic8[[1]][,n+1])) +
  geom_line(aes(seq(1,300), bic19[[1]][,n+1]), colour='grey50', size=trait, linetype=3, data.frame(seq(1,300),bic19[[1]][,n+1])) +
  xlab("# agrégations") + ylab("AUC") +
  theme_bw() +
  annotate("segment", x=230, xend=250, y=.65, yend=.65, colour='black', size=trait, linetype=1) +
  annotate("text", x=260, y=.65, label="AIC 3 var", colour="black", size=4) +
  annotate("segment", x=230, xend=250, y=.64, yend=.64, colour='black', size=trait, linetype=3) +
  annotate("text", x=260, y=.64, label="AIC 4 var", colour="black", size=4) +
  annotate("segment", x=230, xend=250, y=.63, yend=.63, colour='grey60', size=trait, linetype=1) +
  annotate("text", x=260, y=.63, label="BIC 3 var", colour="black", size=4) +
  annotate("segment", x=230, xend=250, y=.62, yend=.62, colour='grey60', size=trait, linetype=5) +
  annotate("text", x=260, y=.62, label="BIC 4 var", colour="black", size=4) +
  annotate("segment", x=230, xend=250, y=.61, yend=.61, colour='grey60', size=trait, linetype=4) +
  annotate("text", x=260, y=.61, label="BIC 8 var", colour="black", size=4) +
  annotate("segment", x=230, xend=250, y=.6, yend=.6, colour='grey50', size=trait, linetype=3) +
  annotate("text", x=262, y=.6, label=" BIC bagging", colour="black", size=4) 

# grille de score
score <- aic3[[2]][,n+1]
sature   <- glm(Cible~.,data=train[,c(varx,"Cible")], family=binomial(link = "logit"))
sature$xlevels
names(unlist(sature$xlevels))
VARIABLE=c("",gsub("[0-9]", "", names(unlist(sature$xlevels))))
MODALITE=c("",as.character(unlist(sature$xlevels)))
names=data.frame(VARIABLE,MODALITE,NOMVAR=c("(Intercept)",paste(VARIABLE,MODALITE,sep="")[-1]))
regression=data.frame(NOMVAR=names(score),COEF=as.numeric(score))
param = merge(names,regression,all.x=TRUE)[-1]
param$COEF[is.na(param$COEF)] <- 0 
mini=aggregate(data.frame(min = param$COEF), by = list(VARIABLE = param$VARIABLE), min)
maxi=aggregate(data.frame(max = param$COEF), by = list(VARIABLE = param$VARIABLE), max)
total=merge(mini,maxi)
total$diff = total$max - total$min
poids_total = sum(total$diff)
grille = merge(param,mini,all.x=TRUE)
grille$delta = grille$COEF - grille$min
grille$POIDS = round((100*grille$delta) / poids_total)
grille[which(VARIABLE!=""),c("VARIABLE","MODALITE","POIDS")]



# ---------------------------------------------------------------------------------------------------------
# Boosting de modèles logistiques
# ---------------------------------------------------------------------------------------------------------

library(ROCR)

# Attention : regrouper la modalité "Autres" de "Objet_credit" qui n'a que 12 observations
# La rareté de cette modalité fait qu'elle peut se trouver absente d'un échantillon
# d'apprentissage alors qu'elle sera présente dans l'échantillon de test

table(credit$Objet_credit)
credit$Objet_credit[credit$Objet_credit == "A410"] <- "A48"
table(credit$Objet_credit)
train  <- credit[id,]
valid  <- credit[-id,]
summary(credit)

library(car)
credit2$Objet_credit <- recode(credit$Objet_credit,"'A40'='Voiture neuve';
'A41'='Voiture occasion';c('A42','A44','A45')='Intérieur';
'A43'='Vidéo - HIFI';c('A46','A48','A410')='Etudes';
'A47'='Vacances';'A49'='Business';else='Autres'")
train  <- credit2[id,]
valid  <- credit2[-id,]

# fonction implémentant le boosting sur un modèle logit

BoostLog = function(apprent, validat, varY, varX, crit="bic", shrink=1, nsimul)
{

nobs <- nrow(apprent)
w <- rep(1/nobs, nobs)
rep <- 0
y <- apprent[,varY]
formule <- as.formula(paste("y ~ ",paste(names(apprent[,varX]),collapse="+")))
sature   <- glm(formule,data=apprent[,c(varX,varY)], family=binomial(link = "logit"))
coef     <- matrix(0,length(coef(sature)),nsimul+2)
rownames(coef) <- names(coef(sature))
predic   <- matrix(NA,nrow(validat),nsimul)
auc      <- matrix(0,nsimul,1)

for(i in (1:nsimul))
{

# initialisation des poids
w <- w/sum(w,na.rm=TRUE) # gérer le cas où un poids w est si petit qu'il vaut NaN

nb_var <- length(varX)
predicteurs <- varX

# tirage aléatoire pondéré avec remise de l'échantillon d'apprentissage
s <- sort(sample(nobs,nobs,replace=T,prob=w))

# formule avec l'ensemble des prédicteurs
if (nb_var > 1) {
formule <- as.formula(paste("y ~ ",paste(names(apprent[,predicteurs]),collapse="+")))
} else {
formule <- as.formula(paste("y ~ ",predicteurs))
}
cible <- apprent[s,varY]
logit <- glm(cible ~ 1,data=apprent[s,], family=binomial(link = "logit"))
if (crit == "bic") {
selection <- step(logit,direction="forward",trace=F,k=log(nobs),scope=list(upper=formule))
} else {
selection <- step(logit,direction="forward",trace=F,k=2,scope=list(upper=formule))
}

# application du modèle logit à l'échantillon d'apprentissage et mise à jour des poids
#pred <- pmin(predict(selection, newdata=apprent, type="response"),0.999999)
pred <- predict(selection, newdata=apprent, type="response")
alpha <- log(pred/(1-pred))/2
signe = ifelse(apprent[,varY] == 1,1,-1)
w <- w * exp(-shrink*signe*alpha)

# application du modèle logit à l'échantillon de validation
#predic[,i] <- min(predict(selection, newdata=validat, type="response"),0.999999)
predic[,i] <- predict(selection, newdata=validat, type="response")
rep <- rep + predic[,i]
cible <- validat[,varY]
pred <- prediction(rep,cible,label.ordering=c(0,1))
auc[i] <- performance(pred,"auc")@y.values[[1]]

# stockage des coefficients
index <- match(names(coef(selection)),names(coef(sature)))
coef[index,i] <- coef(selection)

} # fin de la boucle

#coef <- apply(coef, 2, as.numeric)
coef[,nsimul+1] <- apply(coef[,1:nsimul],1,mean) # moyenne des coefficients
coef[,nsimul+2] <- apply((coef[,1:nsimul]!=0),1,sum) # nb de coefficients non nuls par modalité

# résultats en retour
rep <- rep/nsimul
rf <- list(auc,coef,rep)

} # fin de la fonction

seed <- 235
niter <- 1000
varx <- c(varquali, varquanti)

ptm <- proc.time()
bst <- BoostLog(apprent=train, validat=valid, varY="Cible", varX=varx, crit="bic", shrink=0.001, nsimul=niter)
proc.time() - ptm

# AUC
bst[[1]]
max(bst[[1]])
which.max(bst[[1]])
plot(bst[[1]])
# coefficients moyens
bst[[2]][,niter+1]
# nb de sélections de chaque variable
bst[[2]][,niter+2]
# prédictions du modèle boosté
length(bst[[3]])
library(pROC)
auc(valid$Cible,bst[[3]])
library(ROCR)
pred <- prediction(bst[[3]],valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]



# ---------------------------------------------------------------------------------------------------------
# Version parallélisée
# ---------------------------------------------------------------------------------------------------------

library(parallel)
detectCores()
detectCores(logical=FALSE)

install.packages('doSNOW')
library(doSNOW)
library(foreach)

cl <- makeCluster(4)
registerDoSNOW(cl)
getDoParWorkers()

foreach(i=1:10) %dopar% {

#loop contents here

} 

m <- matrix(rnorm(16), 4, 4)
f <- foreach(i=1:ncol(m),.combine=c) %dopar% {
mean(m[,i])
}
f

stopCluster(cl)


# fonction implémentant le boosting (version //) sur un modèle logit

ParallelBoostLog = function(apprent, validat, varY, varX, crit="bic", shrink=1, nsimul)
{

nobs <- nrow(apprent)
w <- rep(1/nobs, nobs)
rep <- 0
y <- apprent[,varY]
formule <- as.formula(paste("y ~ ",paste(names(apprent[,varX]),collapse="+")))
sature   <- glm(formule,data=apprent[,c(varX,varY)], family=binomial(link = "logit"))
coef     <- matrix(0,length(coef(sature)),nsimul+2)
rownames(coef) <- names(coef(sature))
predic   <- matrix(NA,nrow(validat),nsimul)
auc      <- matrix(0,nsimul,1)

#for(i in (1:nsimul))
#boucle <- foreach(i=1:nsimul,.combine=rbind) %dopar%
boucle <- foreach(i=1:nsimul,.packages='ROCR') %dopar%
{

# initialisation des poids
w <- w/sum(w,na.rm=TRUE) # gérer le cas où un poids w est si petit qu'il vaut NaN
nb_var <- length(varX)
predicteurs <- varX

# tirage aléatoire pondéré avec remise de l'échantillon d'apprentissage
s <- sort(sample(nobs,nobs,replace=T,prob=w))

# formule avec l'ensemble des prédicteurs
if (nb_var > 1) {
formule <- as.formula(paste("y ~ ",paste(names(apprent[,predicteurs]),collapse="+")))
} else {
formule <- as.formula(paste("y ~ ",predicteurs))
}
cible <- apprent[s,varY]
logit <- glm(cible ~ 1,data=apprent[s,], family=binomial(link = "logit"))
if (crit == "bic") {
selection <- step(logit,direction="forward",trace=F,k=log(nobs),scope=list(upper=formule))
} else {
selection <- step(logit,direction="forward",trace=F,k=2,scope=list(upper=formule))
}

# application du modèle logit à l'échantillon d'apprentissage et mise à jour des poids
#pred <- pmin(predict(selection, newdata=apprent, type="response"),0.999999)
pred <- predict(selection, newdata=apprent, type="response")
alpha <- log(pred/(1-pred))/2
signe = ifelse(apprent[,varY] == 1,1,-1)
w <- w * exp(-shrink*signe*alpha)

# stockage des coefficients
index <- match(names(coef(selection)),names(coef(sature)))
coef[index,i] <- coef(selection)

# application du modèle logit à l'échantillon de validation
#predic[,i] <- min(predict(selection, newdata=validat, type="response"),0.999999)
predic[,i] <- predict(selection, newdata=validat, type="response")
rep <- rep + predic[,i]
cible <- validat[,varY]
pred <- prediction(rep,cible,label.ordering=c(0,1))
auc[i] <- performance(pred,"auc")@y.values[[1]]

return(list(auc,coef,rep))

} # fin de la boucle

print(boucle)

auc <- sapply(boucle[[1]],mean)
coef <- boucle[[2]]
rep <- boucle[[3]]

print(auc)
print(coef)
print(rep)

#coef <- apply(coef, 2, as.numeric)
coef[,nsimul+1] <- apply(coef[,1:nsimul],1,mean) # moyenne des coefficients
coef[,nsimul+2] <- apply((coef[,1:nsimul]!=0),1,sum) # nb de coefficients non nuls par modalité

# résultats en retour
rep <- rep/nsimul
rf <- list(auc,coef,rep)

} # fin de la fonction

seed <- 235
niter <- 8
varx <- c(varquali, varquanti)
varx

ptm <- proc.time()
p.bst <- ParallelBoostLog(apprent=train, validat=valid, varY="Cible", varX=varx, crit="bic", shrink=0.001, nsimul=niter)
proc.time() - ptm

# AUC
p.bst[[1]]
class(p.bst[[1]])
max(p.bst[[1]])
which.max(p.bst[[1]])
plot(p.bst[[1]])
# coefficients moyens
p.bst[[2]][,niter+1]
# nb de sélections de chaque variable
p.bst[[2]][,niter+2]


# fin parallélisation
stopCluster(cl)



# ---------------------------------------------------------------------------------------------------------
# Fonction lissant les AUC et les coefficients de forêts aléatoires
# en effectuant plusieurs simulations
# ---------------------------------------------------------------------------------------------------------

Ensemble = function(apprent, validat, varY, varX, nb_var, crit="bic", nsimul, nb)
{

y <- apprent[,varY]
formule <- as.formula(paste("y ~ ",paste(names(apprent[,varX]),collapse="+")))
sature   <- glm(formule,data=apprent[,c(varX,varY)], family=binomial(link = "logit"))
coe <- matrix(0,length(coef(sature)),nb+1)
rownames(coe) <- names(coef(sature))
occ <- matrix(0,length(coef(sature)),nb+1)
rownames(occ) <- names(coef(sature))
auc <- matrix(0,nsimul,nb+1)

for(i in (1:nb))
{
rf <- RForestLog(apprent, validat, varY, varX, nb_var, crit, nsimul)
auc[,i] <- rf[[1]]
coe[,i] <- rf[[2]][,nsimul+1]
occ[,i] <- rf[[2]][,nsimul+2]
}

auc[,nb+1] <- apply(auc[,1:nb],1,mean) # moyenne des AUC
coe[,nb+1] <- apply(coe[,1:nb],1,mean) # moyenne des coefficients
occ[,nb+1] <- apply(occ[,1:nb],1,mean) # moyenne des nb d'occurrences

# résultats en retour
ens <- list(auc,coe,occ)

} # fin de la fonction


niter <- 300
n     <- 30
varx <- c(varquali, varquanti)
#npred <- -grep('(Type_emploi|Telephone|Nb_pers_charge|Anciennete_domicile|Taux_effort|Nb_credits)', varx)
#varx <- varx[npred]
varx

ens <- Ensemble(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=1, crit="aic", nsimul=niter, nb=n)

# évolution de l'AUC
ens[[1]][,nb+1]
plot(ens[[1]][,nb+1])
# coefficients agrégés
ens[[2]][,nb+1]
bic3[[3]][,n+1] 


# grille de score
score <- aic3[[2]][,n+1]
sature   <- glm(Cible~.,data=train[,c(varx,"Cible")], family=binomial(link = "logit"))
sature$xlevels
names(unlist(sature$xlevels))
VARIABLE=c("",gsub("[0-9]", "", names(unlist(sature$xlevels))))
#VARIABLE = c(varquali, varquanti)
VARIABLE
MODALITE=c("",as.character(unlist(sature$xlevels)))
MODALITE
names=data.frame(VARIABLE,MODALITE,NOMVAR=c("(Intercept)",paste(VARIABLE,MODALITE,sep="")[-1]))
names
#regression=data.frame(NOMVAR=names(rf[[2]][,niter+1]),COEF=as.numeric(rf[[2]][,niter+1]))
#regression=data.frame(NOMVAR=names(ens[[2]][,nb+1]),COEF=as.numeric(ens[[2]][,nb+1]))
regression=data.frame(NOMVAR=names(score),COEF=as.numeric(score))
regression
param = merge(names,regression,all.x=TRUE)[-1]
param$COEF[is.na(param$COEF)] <- 0 
param
mini=aggregate(data.frame(min = param$COEF), by = list(VARIABLE = param$VARIABLE), min)
maxi=aggregate(data.frame(max = param$COEF), by = list(VARIABLE = param$VARIABLE), max)
total=merge(mini,maxi)
total$diff = total$max - total$min
poids_total = sum(total$diff)
grille = merge(param,mini,all.x=TRUE)
grille$delta = grille$COEF - grille$min
grille$POIDS = round((100*grille$delta) / poids_total)
grille[which(VARIABLE!=""),c("VARIABLE","MODALITE","POIDS")]



# ---------------------------------------------------------------------------------------------------------
# CHAPITRE 14
# Boosting
# ---------------------------------------------------------------------------------------------------------

library(ada)

# discrete adaboost
names(test)
set.seed(235)
ada.boost <- ada(Cible ~ ., data=credit[id,vars], type="discrete",
loss="exponential", control = rpart.control(cp=0), iter = 5000, nu = 0.01,
test.y=test[,"Cible"], test.x=test[,1:19])
# iter = nb d'arbres agrégés

# real adaboost
set.seed(235)
real.boost.1 <- ada(Cible ~ ., data=credit[id,vars],type="real",
loss="exponential",iter = 5000, nu = 1,
#control = rpart.control(maxdepth = 10, minbucket = 5),
#control = rpart.control(maxdepth=1,cp=-1,minsplit=0),
control = rpart.control(cp=0),
test.y=test[,"Cible"],test.x=test[,1:19])

# pour les stumps : rpart.control(maxdepth=1,cp=-1,minsplit=0)
set.seed(235)
boost <- ada(Cible ~ ., data=credit[id,vars], type="discrete",
loss="exponential", control = rpart.control(maxdepth=1,cp=-1,minsplit=0), iter = 5000, nu = 1,
test.y=test[,"Cible"], test.x=test[,1:19])
boost <- ada(Cible ~ ., data=credit[id,vars], type="real",
loss="exponential", control = rpart.control(maxdepth=1,cp=-1,minsplit=0), iter = 5000, nu = 0.1,
test.y=test[,"Cible"], test.x=test[,1:19])

boost
summary(boost)

# affichage du taux d'erreur
plot(boost, kappa = F, test=T, tflag=F)
title(main = "Real Adaboost - Exponential Loss")

# comparaison des taux d'erreur boosting - random forests
head(boost$model$errs,30)
plot(rf$err.rate[,1],type='l',ylim=c(.22,.38),col='gray60',xlab="# itérations",ylab='erreur')
lines(rf$test$err.rate[,1],type='l',lwd=1,lty=1,col='red')
lines(boost$model$errs[,3],type='l',lwd=2,lty=1,pch=1,col='blue')
lines(ada.boost$model$errs[,3],type='l',lwd=2,lty=1,col='green')
#lines(real.boost.1$model$errs[,3],type='l',lwd=2,lty=1,col='blue')
#lines(real.boost$model$errs[,3],type='l',lwd=2,lty=1,col='green')
legend("topright",c("RF OOB","RF test","AdaBoost test (pénal = 1)","AdaBoost test (pénal = 0.01)"),col=c("gray60","red","blue","green"),lwd=2)

# prédiction
library(ROCR)
test$boost <- predict(boost,test,type='prob')[,2]
pred <- prediction(test$boost,test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

test$boost2000 <- predict(boost,test,type='prob',n.iter=2000)[,2]
pred <- prediction(test$boost2000,test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# simulations sur 30 boostings et calcul AUC moyenne
names(test)
set.seed(235)
nsimul <- 30
auc  <- matrix(NA,nsimul+1,1)
bst  <- matrix(NA,nrow(test),nsimul+1)
for(i in 1:nsimul)
{
rm(boost)
boost <- ada(Cible ~ ., data=credit[id,vars],type="real",
loss="logistic",
control = rpart.control(maxdepth = 10, minbucket = 5),iter = 2000, nu=0.01,
test.y=test[,"Cible"],test.x=test[,1:19])
bst[,i] <- predict(boost,test,type='prob')[,2]
pred <- prediction(bst[,i],test$Cible,label.ordering=c(0,1))
auc[i,1] <- performance(pred,"auc")@y.values[[1]]
} # fin de la boucle
bst[,nsimul+1] <- apply(bst[,1:nsimul],1,mean) # moyenne des prédictions
pred <- prediction(bst[,nsimul+1],test$Cible,label.ordering=c(0,1))
auc[nsimul+1,1] <- performance(pred,"auc")@y.values[[1]]
auc
mean(auc[1:30])

# importance des variables
varplot(boost,type="scores")
vip <- varplot(boost,type="scores")
vip
par(mar = c(8, 4, 4, 0))
barplot(vip, col = gray(0:length(vip)/length(vip)),
        ylab='Importance', ylim=c(0,0.03),cex.names = 0.8, las=3)
title(main = list("Ada : Importance des variables", font = 1, cex = 1.2))

# importance moyenne des variables
nvar <- ncol(credit[,vars])-1
vimp <- rep(0,nvar)
nsimul <- 10
t1 <- proc.time()
for(i in 1:nsimul){
boost <- ada(Cible ~ ., data=credit[id,vars],type="real",
control = rpart.control(maxdepth=1,cp=-1,minsplit=0),iter = 1000, nu = 0.01)
vip <- varplot(boost,type="scores")
vimp <- vimp + as.numeric(vip[order(names(vip))]) / nsimul
cat("i=", i, " time=", (proc.time()-t1)/60, "\n")
rm(boost)
}
a1 <- sort(names(vip))
a2 <- order(vimp, decreasing = TRUE)
dotchart(vimp[a2][nvar:1], a1[a2][nvar:1], main="Ada : Importance moyenne des variables")

# correctif de Mark Culp
# pondère l'importance des variables par le poids x$model$alpha des arbres
# dans l'agrégation finale
varplot2<-function (x, plot.it = TRUE, type = c("none", "scores"), max.var.show = 30, 
    ...) 
{
    if (class(x) != "ada") {
        stop("Object must be of type 'ada'")
    }
    if (missing(type)) {
        type = "none"
    }
    iter <- x$iter
    nm <- x$names
    vec <- rep(0, length(nm))
    p = x$dim[2]
    g1 <- function(i, obj,alp) {
        if (dim(obj[[i]]$frame)[1] < 2) {
            return(rep(0, p))
        }
        imp <- alp*obj[[i]]$splits[, 3]^2
        vals <- match(row.names(obj[[i]]$splits), nm)
        vec = rep(0, p)
        vec[vals] <- imp
        vec
    }
    alp=x$model$alpha
    vec <- 1/iter * sqrt(apply(sapply(1:iter, function(i) g1(i, 
        x$model$trees,alp[i])), 1, sum))
    vars <- order(vec, decreasing = TRUE)
    n <- length(vec)
    max.v = max.var.show
    if (p < max.v) 
        max.v = p
    if (plot.it == TRUE) {
        dotchart(vec[vars[max.v:1]], labels = nm[vars[max.v:1]], 
            xlab = "Score", main = "Variable Importance Plot")
    }
    if (type == "scores") {
        vars = vars[1:max.v]
        t1 <- vec[vars]
        attr(t1, "names") <- nm[vars]
        return(t1)
    }
}
varplot2(boost,plot.it=T,type="scores")



# ---------------------------------------------------------------------------------------------------------
# Boosting avec package gbm
# ---------------------------------------------------------------------------------------------------------

library(gbm)
library(ada)

# rappel : real adaboost avec "ada"
ptm <- proc.time()
mada <- ada(Cible ~ ., data=credit[id,vars],type="real",
loss="exponential",iter = 3000, nu = 0.01,
#control = rpart.control(maxdepth = 10, minbucket = 5),
#control = rpart.control(maxdepth=1,cp=-1,minsplit=0),
control = rpart.control(cp=0),
test.y=test[,"Cible"],test.x=test[,1:19])
proc.time() - ptm
plot(mada$model$alpha)

# affichage de la courbe du taux d'erreur
plot(mada, kappa = F, test=T, tflag=F)
title(main = "Real Adaboost avec package ada - Exponential Loss")

# prédiction avec ada
library(ROCR)
ada2000 <- predict(mada,test,type='prob',n.iter=2000)[,2]
pred <- prediction(ada2000,test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# utilisation gbm
ptm <- proc.time()
mgbm <- gbm(Cible ~ ., data=credit[id,vars],distribution="adaboost", n.trees = 3000, shrinkage = 0.01,
n.minobsinnode = 5,cv.folds = 5,verbose=T)
proc.time() - ptm
summary(mgbm,n.trees=2000,method=relative.influence) # plot variable influence
gbm.perf(mgbm,oobag.curve = T,overlay = TRUE)
# importance des variables
plot(mgbm,i.var=1:2)

pretty.gbm.tree(mgbm, i.tree = 100)

# check performance using an out-of-bag estimator
# OOB underestimates the optimal number of iterations
(best.iter <- gbm.perf(mgbm,method="OOB"))
# check performance using 5-fold cross-validation
(best.iter <- gbm.perf(mgbm,method="cv"))

# prédiction
gbm2000 <- predict(mgbm,test,n.trees=400)
# The adaboost method gives the predictions on logit scale. You can convert it to the 0-1 output:
gbm2000 <- plogis(gbm2000)
table(gbm2000)
pred <- prediction(gbm2000,test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# affichage taux d'erreur
B <- 3000
boucle <- seq(1,B,by=50)
errapp <- rep(0,length(boucle))
errtest <- errapp
k <- 0
for (i in boucle){
k <- k+1
prev_app <- predict(mgbm,newdata=credit[id,vars],n.trees=i)
errapp[k] <- sum(as.numeric(prev_app>0)!=credit[id,"Cible"])/nrow(credit[id,])
prev_test<- predict(mgbm,newdata=test,n.trees=i)
errtest[k] <- sum(as.numeric(prev_test>0)!=test$Cible)/nrow(test)
}

plot(boucle,errapp,type="l",col="blue",ylim=c(0,0.8),xlab="nombre iterations",ylab="erreur")
lines(boucle,errtest,col="red")



# ---------------------------------------------------------------------------------------------------------
# CHAPITRE 15
# SVM (package e1071)
# ---------------------------------------------------------------------------------------------------------

library(e1071)

# SVM linéaire
# -----------------------------------------------

ptm <- proc.time()
svmlin = svm(Cible ~ ., data=credit[id,vars],kernel="linear",probability=TRUE,cost=1)
proc.time() - ptm
print(svmlin)
summary(svmlin)
svmlin$coefs
svmlin$index # indices des vecteurs supports
svmlin$rho
svmlin$SV
table(svmlin$fitted,credit[id,"Cible"])

# détermination des coefficients des prédicteurs - avec e1071
names(credit)
x <- model.matrix( ~ . -1 ,data=credit[id,1:19])
y <- credit[id,20]
svmlin = svm(x, y,kernel="linear",probability=TRUE,cost=1,scale=T)
# w <- t(svmlin$coefs) %*% x[svmlin$index,]
# w <- t(svmlin$coefs) %*% scale(x[svmlin$index,])
# la ligne précédente ne convient pas car la standardisation ne s'effectue
# qu'en prenant en compte l'ensembles des points supports et pas l'ensemble des points
w <- t(svmlin$coefs) %*% svmlin$SV
# la ligne précédente s'appuie sur les vecteurs supports convenablement centrés-réduits
w <- t(svmlin$coefs) %*% scale(x[svmlin$index,],apply(x,2,mean),apply(x,2,sd))
# sinon on peut centrer-réduire avec la moyenne et l'écart-type calculés sur tous les points de x
# et pas seulement les points supports
# calcul du score pour données non centrées-réduites
p2 <- x %*% t(w) - svmlin$rho
# calcul du score pour données centrées-réduites
xscaled <- scale(x,apply(x,2,mean),apply(x,2,sd))
#xscaled <- scale(x,svmlin$x.scale[[1]],svmlin$x.scale[[2]])
p2 <- xscaled %*% t(w) - svmlin$rho
p2 <- x %*% t(w/apply(x,2,sd)) - svmlin$rho
svmlin$rho
# coefficients du score SVM sur prédicteurs initiaux (avant standardisation)
w/svmlin$x.scale[[2]]
w/apply(x,2,sd)
head(p2)
# score calculé par e1071
p1 <- predict(svmlin, newdata=x, decision.values=T)
p1 <- attr(p1, "decision.values")
head(p1)
# comparaison des scores
cor(p1,p2,method="spearman")
max(abs(p1-p2))
plot(p1,p2)

# détermination des coefficients des prédicteurs - avec LiblineaR - noyau linéaire
#install.packages('LiblineaR')
library('LiblineaR')
s <- scale(x)
(co <- heuristicC(s))
m <- LiblineaR(data=s,labels=y,type=3,cost=1,bias=TRUE,verbose=T)
m <- LiblineaR(data=s,labels=y,type=3,cost=1,bias=TRUE,verbose=T,epsilon=0.001)
m$W # coefficients du modèle SVM à noyau linéaire
(rho <- m$W[length(m$W)]) # constante (intercept)
# application du score SVM sur données centrées-réduites
#st <- scale(xt,apply(x,2,mean),apply(x,2,sd))
#st <- scale(xt,attr(s,"scaled:center"),attr(s,"scaled:scale"))
#p1 <- s %*% as.matrix(m$W[1:length(m$W)-1]) + rho
p1 <- cbind(s,1) %*% t(m$W)
head(p1)
cor(p1,p2,method="pearson")
cor(p1,p2,method="spearman")
plot(p1,p2)
# coefficients du score SVM sur prédicteurs initiaux (avant standardisation)
m$W/apply(x,2,sd)

# détermination des coefficients des prédicteurs - avec kernlab - noyau linéaire
library(kernlab)
ksvmlin <- ksvm(x,y,kernel="vanilladot",type = "C-svc",scaled=T,C=1,prob.model = TRUE)
ksvmlin
w <- t(unlist(ksvmlin@coef)) %*% xmatrix(ksvmlin)[[1]]
# xmatrix contient la matrice des vecteurs-supports utilisée pour les calculs 
# (après standardisation éventuelle)
# variante :
w <- t(unlist(ksvmlin@coef)) %*% scale(x[unlist(ksvmlin@alphaindex),],apply(x,2,mean),apply(x,2,sd))
xscaled <- scale(x,apply(x,2,mean),apply(x,2,sd))
p2 <- xscaled %*% t(w) - ksvmlin@b
head(p2)
# comparaison des prédictions
predksvm = predict(ksvmlin,type="prob",x)
head(predksvm)
cor(predksvm[,2],p2,method="spearman")
plot(predksvm[,2],p2)
# passage des coefficients sur prédicteurs centrés-réduits
# à des coefficients sur prédicteurs initiaux
w/apply(x,2,sd)

# calcul du taux d'erreur
# par resubstitution 1
svmlin$fitted
table(svmlin$fitted,credit[id,"Cible"])
mean(svmlin$fitted != credit[id,"Cible"]) # 0.2003106
# par validation croisée
set.seed(235)
svmlin = svm(Cible ~ ., data=credit[id,vars],kernel="linear",probability=TRUE,cost=1,cross=10)
summary(svmlin)
# AUC sur classe prédite
library(ROCR)
pred=prediction(as.numeric(svmlin$fitted),credit[id,"Cible"],label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.7324943
# AUC sur probabilités des classes (= score)
predsvmlin = predict(svmlin,type="prob",credit[id,vars],probability=TRUE)
pred=prediction(attr(predsvmlin,"probabilities")[,1],credit[id,"Cible"],label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.8466505
# taux d'erreur par resubstitution (2e calcul)
mean(predsvmlin != credit[id,"Cible"]) # 0.2034161
# écart entre valeurs "predict" et "fitted"
mean(predsvmlin != svmlin$fitted) # 0.05279503
table(predsvmlin,svmlin$fitted)

# Generates a scatter plot of the input data of a svm fit for classification models
# by highlighting the classes and support vectors. Optionally, draws a filled contour plot of the class regions.
plot(svmlin,credit,Duree_credit ~ Age, svSymbol="x",dataSymbol="o",grid=200)

# application sur échantillon test
library(ROCR)
testsvmlin = predict(svmlin,type="prob",test,probability=TRUE)
pred=prediction(attr(testsvmlin,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.7338138
# taux d'erreur en test
mean(testsvmlin != test$Cible) # 0.252809

# recherche meilleur paramètre de coût par validation croisée
set.seed(235)
svmlin = tune.svm(Cible ~ ., data=credit[id,vars],kernel="linear",cost=10^(-3:1))
summary(svmlin)
plot(svmlin)
svmlin$best.parameters
svmlin$best.performance
svmlin$train.ind

# recherche meilleur paramètre de coût sur données de test
names(test)
svmlin = tune.svm(Cible ~ ., data=credit[id,vars],kernel="linear",cost=10^(-3:1),
validation.x = test [,-20], validation.y = test [,20],
tunecontrol = tune.control(sampling = "fix", fix=1))
summary(svmlin)

plot(svmlin)
plot(svmlin, transform.x = log2, transform.y = log2)
plot(svmlin, type = "perspective", theta = 120, phi = 45)

# AUC en fonction du coût
svmlin = svm(Cible ~ ., data=credit[id,vars],kernel="linear",probability=TRUE,cost=0.1)
# apprentissage
pred=prediction(ifelse(svmlin$fitted=="1",1,0),credit[id,"Cible"],label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]
# test
test$svmlin = predict(svmlin,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmlin,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]


# SVM linéaire sur données discrétisées
# -----------------------------------------------

# AUC en fonction du coût (data frame credit avec Anciennete_emploi en factor)
svmlin = svm(Cible ~ ., data=credit[id,vars],kernel="linear",probability=TRUE,cost=0.1)
test$svmlin = predict(svmlin,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmlin,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# AUC en fonction du coût (data frame credit2 après discrétisation)
svmlin = svm(Cible ~ ., data=train,kernel="linear",probability=TRUE,cost=0.1)
valid$svmlin = predict(svmlin,type="prob",valid,probability=TRUE)
pred=prediction(attr(valid$svmlin,"probabilities")[,1],valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# recherche meilleur paramètre de coût sur données de test
svmlin = tune.svm(Cible ~ ., data=train,kernel="linear",cost=10^(-3:1),
validation.x = valid [,-c(17,21)], validation.y = valid [,17],
tunecontrol = tune.control(sampling = "fix", fix=1))
summary(svmlin)
svmlin$terms


# SVM avec noyau radial
# -----------------------------------------------

# ajustement des paramètres gamma et C d’une SVM  par validation croisée 10-fold,
# en faisant varier le paramètre dans (0.001, 0.01, 0.1, 1) et la constante de régularisation C de 1 à 10.
svmrad=tune.svm(Cible ~ ., data=credit[id,vars],kernel="radial",cost=seq(0.1,10,by=0.1),gamma=c(0.001,0.01,0.1,1),
validation.x = test [,-20], validation.y = test [,20],
tunecontrol = tune.control(sampling = "fix", fix=1))
summary(svmrad)

svmrad=tune.svm(Cible ~ ., data=credit[id,vars],kernel="radial",cost=seq(0.1,5,by=0.1),gamma=seq(0.01,0.5,by=0.01),
validation.x = test [,-20], validation.y = test [,20],
tunecontrol = tune.control(sampling = "fix", fix=1))
summary(svmrad)
plot(svmrad)

# recherche paramètres optimaux par validation croisée
set.seed(235)
svmrad=tune.svm(Cible ~ ., data=credit[id,vars],kernel="radial",cost=1:10,gamma=c(0.001,0.01,0.1,1))
summary(svmrad)

# recherche paramètres optimaux par validation croisée
set.seed(235)
svmrad=tune.svm(Cible ~ ., data=credit[id,vars],kernel="radial",cost=seq(0.01,5,by=0.01),gamma=c(0.01,0.1,1))
summary(svmrad)

# on retient les paramètres gamma = 0.08 et cost = 1.5
svm <- svm(Cible ~ ., data=credit[id,vars],kernel="radial",gamma=0.08,cost=1.5,probability=TRUE)
test$svmrad = predict(svm,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmrad,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.7689449

# on retient les paramètres gamma = 0.1 et cost = 1
svm <- svm(Cible ~ ., data=credit[id,vars],kernel="radial",gamma=0.1,cost=1,probability=TRUE)
svmrad = predict(svm,type="prob",test,probability=TRUE)
pred=prediction(attr(svmrad,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.7691326

# on retient les paramètres par défaut : gamma = 0.02083333 et cost = 1
svm <- svm(Cible ~ ., data=credit[id,vars],kernel="radial",probability=TRUE)
summary(svm)
test$svmrad = predict(svm,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmrad,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.7643471

# on retient les paramètres gamma = 0.1 et cost = 2
svm <- svm(Cible ~ ., data=credit[id,vars],kernel="radial",gamma=0.1,cost=2,probability=TRUE)
test$svmrad = predict(svm,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmrad,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.7591863

# on retient les paramètres gamma = 0.1 et cost = 1.52
svm <- svm(Cible ~ ., data=credit[id,vars],kernel="radial",gamma=0.1,cost=1.52,probability=TRUE)
test$svmrad = predict(svm,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmrad,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.7642157

# on retient les paramètres gamma = 0.0643 et cost = 1.2
svm <- svm(Cible ~ ., data=credit[id,vars],kernel="radial",gamma=0.0643,cost=1.2,probability=TRUE)
test$svmrad = predict(svm,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmrad,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.7720227


# SVM avec noyau radial sur données discrétisées
# -----------------------------------------------

# ajustement des paramètres gamma et C d’une SVM  par validation sur données de test,
# en faisant varier le paramètre dans (0.001, 0.01, 0.1, 1) et la constante de régularisation C de 1 à 10.
svmrad=tune.svm(Cible ~ ., data=train,kernel="radial",cost=seq(0.1,10,by=0.1),gamma=c(0.001,0.01,0.1,1),
validation.x = valid [,-c(17,21)], validation.y = valid [,17],
tunecontrol = tune.control(sampling = "fix", fix=1))
summary(svmrad)

svmrad=tune.svm(Cible ~ ., data=train,kernel="radial",cost=seq(0.1,5,by=0.1),gamma=seq(0.01,0.5,by=0.01),
validation.x = valid [,-c(17,21)], validation.y = valid [,17],
tunecontrol = tune.control(sampling = "fix", fix=1))
summary(svmrad)
plot(svmrad)

# on retient les paramètres gamma = 0.01 et cost = 2.8
svm <- svm(Cible ~ ., data=train,kernel="radial",gamma=0.01,cost=2.8,probability=TRUE)
valid$svmrad = predict(svm,type="prob",valid,probability=TRUE)
pred=prediction(attr(valid$svmrad,"probabilities")[,1],valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.7548699

# on retient les paramètres gamma = 0.04 et cost = 1.3
svm <- svm(Cible ~ ., data=train,kernel="radial",gamma=0.04,cost=1.3,probability=TRUE)
valid$svmrad = predict(svm,type="prob",valid,probability=TRUE)
pred=prediction(attr(valid$svmrad,"probabilities")[,1],valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.7520549


# SVM sigmoïde
# -----------------------------------------------

svmsig = svm(Cible ~ ., data=credit[id,vars],kernel="sigmoid",probability=TRUE)
summary(svmsig)

test$svmsig = predict(svmsig,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmsig,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.7440604

# ajustement des paramètres sur l'échantillon de test
svmsig = tune.svm(Cible ~ ., data=credit[id,vars],kernel="sigmoid",cost=seq(1,100,by=0.5),gamma=seq(0.001,0.05,by=0.001),
validation.x = test [,-20], validation.y = test [,20],
tunecontrol = tune.control(sampling = "fix", fix=1))
summary(svmsig)
svmsig$best.model

# ajustement des paramètres par validation croisée
set.seed(235)
tune.svm(Cible ~ ., data=credit[id,vars],kernel="sigmoid",cost=1:100,gamma=c(0.001,0.01,0.1,1))

# on retient les paramètres gamma = 0.01 et cost = 6.5
svm <- svm(Cible ~ ., data=credit[id,vars],kernel="sigmoid",gamma=0.01,cost=6.5,probability=TRUE)
test$svmsig = predict(svm,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmsig,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.743929

# on retient les paramètres gamma = 0.05 et cost = 1.5
svm <- svm(Cible ~ ., data=credit[id,vars],kernel="sigmoid",gamma=0.05,cost=1.5,probability=TRUE)
test$svmsig = predict(svm,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmsig,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.7319183

# on retient les paramètres gamma = 0.01 et cost = 19
svm <- svm(Cible ~ ., data=credit[id,vars],kernel="sigmoid",gamma=0.01,cost=19,probability=TRUE)
test$svmsig = predict(svm,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmsig,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.7370792


# SVM avec noyau sigmoïde sur données discrétisées
# -----------------------------------------------

names(valid)
valid$svmsig <- NULL

# ajustement des paramètres gamma et C d’une SVM  par validation sur données de test,
svmsig=tune.svm(Cible ~ ., data=train,kernel="sigmoid",cost=seq(0.5,100,by=0.5),gamma=c(0.0001,0.001,0.01,0.1,1),
validation.x = valid [,-17], validation.y = valid [,17],
tunecontrol = tune.control(sampling = "fix", fix=1))
summary(svmsig)
svmsig$best.model

# ajustement des paramètres gamma et C d’une SVM  par validation sur données de test,
svmsig=tune.svm(Cible ~ ., data=train,kernel="sigmoid",cost=seq(0.5,100,by=0.5),gamma=seq(0.001,0.01,by=0.001),
validation.x = valid [,-17], validation.y = valid [,17],
tunecontrol = tune.control(sampling = "fix", fix=1))
summary(svmsig)
svmsig$best.model

# on retient les paramètres gamma = 0.001 et cost = 50
svm <- svm(Cible ~ ., data=train,kernel="sigmoid",gamma=0.001,cost=50,probability=TRUE)
valid$svmsig = predict(svm,type="prob",valid,probability=TRUE)
pred=prediction(attr(valid$svmsig,"probabilities")[,1],valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.7585482

# on retient les paramètres gamma = 0.005 et cost = 10
svm <- svm(Cible ~ ., data=train,kernel="sigmoid",gamma=0.005,cost=10,probability=TRUE)
valid$svmsig = predict(svm,type="prob",valid,probability=TRUE)
pred=prediction(attr(valid$svmsig,"probabilities")[,1],valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.7584356


# SVM avec noyau polynomial de degré 2
# -----------------------------------------------

svmpol2=svm(Cible ~ ., data=credit[id,vars],kernel="polynomial",degree=2,cross=10)
summary(svmpol2)
svmpol2$coefs
plot(svmpol2,credit,Duree_credit ~ Age,svSymbol=16,dataSymbol=1,grid=200)

svmpol = tune.svm(Cible ~ ., data=credit[id,vars],kernel="polynomial",degree=2,gamma=seq(0.01,1,by=0.01),cost=seq(0.1,10,by=0.1),
validation.x = test [,-20], validation.y = test [,20],
tunecontrol = tune.control(sampling = "fix", fix=1))
svmpol$best.model

# affinage de la recherche en incluant la recherche du coefficient coef0
svmpol = tune.svm(Cible ~ ., data=credit[id,vars],kernel="polynomial",degree=2,gamma=seq(0.01,0.1,by=0.01),cost=seq(0.1,10,by=0.1),coef0=seq(0,5,by=0.1),
validation.x = test [,-20], validation.y = test [,20],
tunecontrol = tune.control(sampling = "fix", fix=1))
#summary(svmpol)
svmpol$best.model

# on retient les paramètres gamma = 0.03 et cost = 2.9
svm <- svm(Cible ~ ., data=credit[id,vars],kernel="polynomial",degree=2,gamma=0.03,cost=2.9,probability=TRUE)
test$svmpol2 = predict(svm,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmpol2,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# on retient les paramètres gamma = 0.03 et cost = 1.2 et coef0 = 1
svm <- svm(Cible ~ ., data=credit[id,vars],kernel="polynomial",degree=2,gamma=0.03,cost=1.2,coef0=1,probability=TRUE)
test$svmpol2 = predict(svm,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmpol2,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# on retient les paramètres gamma = 0.01 et cost = 0.6 et coef0 = 4.8
svm <- svm(Cible ~ ., data=credit[id,vars],kernel="polynomial",degree=2,gamma=0.01,cost=0.6,coef0=4.8,probability=TRUE)
test$svmpol2 = predict(svm,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmpol2,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# on retient les paramètres gamma = 0.02 et cost = 2.9 et coef0 = -0.005
svm <- svm(Cible ~ ., data=credit[id,vars],kernel="polynomial",degree=2,gamma=0.02,cost=2.9,coef0=-0.005,probability=TRUE)
test$svmpol2 = predict(svm,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmpol2,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]


# SVM avec noyau polynomial de degré 2 sur données discrétisées
# -----------------------------------------------

svmpol = tune.svm(Cible ~ ., data=train,kernel="polynomial",degree=2,gamma=seq(0.01,1,by=0.01),cost=seq(0.1,10,by=0.1),
validation.x = valid [,-c(17,21,22)], validation.y = valid [,17],
tunecontrol = tune.control(sampling = "fix", fix=1))
summary(svmpol)
svmpol$best.model

# on retient les paramètres gamma = 0.11 et cost = 0.9
svm <- svm(Cible ~ ., data=train,kernel="polynomial",degree=2,gamma=0.11,cost=0.9,probability=TRUE)
valid$svmpol2 = predict(svm,type="prob",valid,probability=TRUE)
pred=prediction(attr(valid$svmpol2,"probabilities")[,1],valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# affinage de la recherche en incluant la recherche du coefficient coef0
svmpol = tune.svm(Cible ~ ., data=train,kernel="polynomial",degree=2,gamma=seq(0.05,0.2,by=0.01),cost=seq(0.1,2,by=0.1),coef0=seq(0,1,by=0.1),
validation.x = valid [,-17], validation.y = valid [,17],
tunecontrol = tune.control(sampling = "fix", fix=1))
svmpol$best.model

# on retient les paramètres gamma = 0.09 et cost = 0.3 et coef0=1
svm <- svm(Cible ~ ., data=train,kernel="polynomial",degree=2,gamma=0.09,cost=0.3,coef0=1,probability=TRUE)
valid$svmpol2 = predict(svm,type="prob",valid,probability=TRUE)
pred=prediction(attr(valid$svmpol2,"probabilities")[,1],valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] # 0.7565965


# SVM avec noyau polynomial de degré 3
# -----------------------------------------------

svmpol = tune.svm(Cible ~ ., data=credit[id,vars],kernel="polynomial",degree=3,gamma=seq(0.01,1,by=0.01),cost=seq(0.1,10,by=0.1),
validation.x = test [,-20], validation.y = test [,20],
tunecontrol = tune.control(sampling = "fix", fix=1))
svmpol$best.model

# on retient les paramètres gamma = 0.07 et cost = 1.9
# à noter que ces paramètres sont proches de ceux du noyau de degré 2
svm <- svm(Cible ~ ., data=credit[id,vars],kernel="polynomial",degree=3,gamma=0.07,cost=1.9,probability=TRUE)
test$svmpol3 = predict(svm,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmpol3,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# on retient les paramètres gamma = 0.07 et cost = 1.4
# à noter que ces paramètres sont proches de ceux du noyau de degré 2
svm <- svm(Cible ~ ., data=credit[id,vars],kernel="polynomial",degree=3,gamma=0.07,cost=1.4,probability=TRUE)
test$svmpol3 = predict(svm,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmpol3,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# affinage de la recherche en incluant la recherche du coefficient coef0
svmpol = tune.svm(Cible ~ ., data=credit[id,vars],kernel="polynomial",degree=3,gamma=seq(0.01,0.1,by=0.01),cost=seq(0.1,10,by=0.1),coef0=seq(0,1,by=0.1),
validation.x = test [,-20], validation.y = test [,20],
tunecontrol = tune.control(sampling = "fix", fix=1))
svmpol$best.model

# on retient les paramètres gamma = 0.03 et cost = 0.6 et coef0 = 0.9
svm <- svm(Cible ~ ., data=credit[id,vars],kernel="polynomial",degree=3,gamma=0.03,cost=0.6,coef0=0.9,probability=TRUE)
test$svmpol3 = predict(svm,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmpol3,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# on retient les paramètres gamma = 0.049 et cost = 0.99 et coef0 = 0.076
svm <- svm(Cible ~ ., data=credit[id,vars],kernel="polynomial",degree=3,gamma=0.049,cost=0.99,coef0=0.076,probability=TRUE)
test$svmpol3 = predict(svm,type="prob",test,probability=TRUE)
pred=prediction(attr(test$svmpol3,"probabilities")[,1],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# risque estimé
mean(test$Cible!=test$svmpol3)


# SVM avec noyau polynomial de degré 3 sur données discrétisées
# -----------------------------------------------

svmpol = tune.svm(Cible ~ ., data=train,kernel="polynomial",degree=3,gamma=seq(0.01,1,by=0.01),cost=seq(0.1,10,by=0.1),
validation.x = valid [,-c(17,21,22,23)], validation.y = valid [,17],
tunecontrol = tune.control(sampling = "fix", fix=1))
#summary(svmpol)
svmpol$best.model

# on retient les paramètres gamma = 0.2 et cost = 0.1
svm <- svm(Cible ~ ., data=train,kernel="polynomial",degree=3,gamma=0.2,cost=0.1,probability=TRUE)
valid$svmpol3 = predict(svm,type="prob",valid,probability=TRUE)
pred=prediction(attr(valid$svmpol3,"probabilities")[,1],valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# affinage de la recherche en incluant la recherche du coefficient coef0
svmpol = tune.svm(Cible ~ ., data=train,kernel="polynomial",degree=3,gamma=seq(0.01,0.5,by=0.01),cost=seq(0.05,1,by=0.05),coef0=seq(0,1,by=0.1),
validation.x = valid[, -c(17, 21,22)], validation.y = valid[,17],
tunecontrol = tune.control(sampling = "fix", fix=1))
svmpol$best.model

# on retient les paramètres gamma = 0.05 et cost = 0.3 et coef0 = 1
svm <- svm(Cible ~ ., data=train,kernel="polynomial",degree=3,gamma=0.05,cost=0.3,coef0 =1, probability=TRUE)
valid$svmpol3 = predict(svm,type="prob",valid,probability=TRUE)
pred=prediction(attr(valid$svmpol3,"probabilities")[,1],valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# on retient les paramètres gamma = 0.04 et cost = 0.4 et coef0 = 1
svm <- svm(Cible ~ ., data=train,kernel="polynomial",degree=3,gamma=0.04,cost=0.4,coef0 =1, probability=TRUE)
valid$svmpol3 = predict(svm,type="prob",valid,probability=TRUE)
pred=prediction(attr(valid$svmpol3,"probabilities")[,1],valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]



# ---------------------------------------------------------------------------------------------------------
# SVM (package kernlab)
# ---------------------------------------------------------------------------------------------------------

library(kernlab)

# kernlab - noyau linéaire
# -----------------------------------------------

ptm <- proc.time()
ksvmlin <- ksvm(Cible ~ ., data=credit[id,vars],kernel="vanilladot",
type = "C-svc",probability=TRUE,C=1,prob.model = TRUE)
proc.time() - ptm
ksvmlin
kernelf(ksvmlin)
ksvmlin@alpha
ksvmlin@coef
ksvmlin@prob.model
# prédiction
predksvm = predict(ksvmlin,type="prob",test)
head(predksvm)
library(ROCR)
pred=prediction(predksvm[,2],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]
# corrélation entre prédiction par kernlab et e1071
cor(predksvm[,2],attr(testsvmlin,"probabilities")[,1])


# kernlab - noyau personnalisé
# -----------------------------------------------

k <- function(x, y) {
 (sum(x * y) + 1) * exp(0.001 * sum((x - y)^2))
 }
class(k) <- "kernel"
ptm <- proc.time()
ksvmuser = ksvm(Cible ~ .,data=credit[id,vars],kernel=k,
type = "C-svc",probability=TRUE,C=1,prob.model = TRUE)
proc.time() - ptm
ksvmuser
predkuk <- outer(i,j,Vectorize(f))
alpha1 <- i[which(k==max(k),arr.ind=TRUE)[1,1]]
alpha2 <- j[which(k==max(k),arr.ind=TRUE)[1,2]]
cat("\n","nb facteurs = ",n," AUC test max = ",max(k)," alpha peel = ",
alpha1," alpha paste = ",alpha2)ser = predict(ksvmuser,type="prob",test)
pred=prediction(predkuser[,2],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]


# kernlab - noyau radial gaussien
# -----------------------------------------------

library(kernlab)
ptm <- proc.time()
ksvmrbf = ksvm(Cible ~ ., data=credit[id,vars],kernel="rbfdot",
type = "C-svc",probability=TRUE,C=1,prob.model = TRUE,kpar=list(sigma=0.1))
ksvmrbf = ksvm(Cible ~ ., data=credit[id,vars],kernel="rbfdot",
type = "C-svc",probability=TRUE,C=1.2,prob.model = TRUE,kpar=list(sigma=0.0643))
proc.time() - ptm
ksvmrbf
predkrbf = predict(ksvmrbf,type="prob",test)
head(predkrbf)
library(ROCR)
pred=prediction(predkrbf[,2],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]
# corrélation entre prédiction par kernlab et e1071
cor(predkrbf[,2],attr(svmrad,"probabilities")[,1])

# maximisation AUC en test
maxauc = function(x)
{
ksvmlap = ksvm(Cible ~ ., data=credit[id,vars],kernel="rbfdot",
type = "C-svc",C=x[1],prob.model = TRUE,kpar=list(sigma=x[2]))
predklap = predict(ksvmlap,type="prob",test)
pred=prediction(predklap[,2],test$Cible,label.ordering=c(0,1))
-performance(pred,"auc")@y.values[[1]]
}

algo <- "Nelder-Mead"
est = optim(c(1,0.05),maxauc,method=algo)
est$convergence
est$par
est$value


# kernlab - noyau radial laplacien
# -----------------------------------------------

library(kernlab)
ptm <- proc.time()
ksvmlap = ksvm(Cible ~ ., data=credit[id,vars],kernel="laplacedot",
type = "C-svc",prob.model = TRUE)
ksvmlap = ksvm(Cible ~ ., data=credit[id,vars],kernel="laplacedot",
type = "C-svc",C=1.79,prob.model = TRUE,kpar=list(sigma=0.314))
proc.time() - ptm
ksvmlap
predklap = predict(ksvmlap,type="prob",test)
head(predklap)
library(ROCR)
pred=prediction(predklap[,2],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]
# corrélation entre prédiction par kernlab et e1071
cor(predklap[,2],attr(svmrad,"probabilities")[,1])

# maximisation AUC en test
maxauc = function(x)
{
ksvmlap = ksvm(Cible ~ ., data=credit[id,vars],kernel="laplacedot",
type = "C-svc",C=x[1],prob.model = TRUE,kpar=list(sigma=x[2]))
predklap = predict(ksvmlap,type="prob",test)
pred=prediction(predklap[,2],test$Cible,label.ordering=c(0,1))
-performance(pred,"auc")@y.values[[1]]
}

algo <- "Nelder-Mead"
algo <- "BFGS"
algo <- "L-BFGS-B"
algo <- "CG"
est = optim(c(1,0.05),maxauc,method=algo)
est$convergence
est$par
est$value


# kernlab - noyau polynomial
# -----------------------------------------------

library(kernlab)
ptm <- proc.time()
ksvmpol = ksvm(Cible ~ ., data=credit[id,vars],kernel="polydot",
type = "C-svc",probability=TRUE,C=2.9,prob.model = TRUE,
kpar=list(degree=2,scale=0.03,offset=0))
proc.time() - ptm
ksvmpol
predkpol = predict(ksvmpol,type="prob",test)
pred=prediction(predkpol[,2],test$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# maximisation AUC
maxauc = function(x)
{
ksvmpol = ksvm(Cible ~ ., data=credit[id,vars],kernel="polydot",
type = "C-svc",C=x[1],prob.model = TRUE,
kpar=list(degree=3,scale=x[2],offset=x[3]))
predkpol = predict(ksvmpol,type="prob",test)
pred=prediction(predkpol[,2],test$Cible,label.ordering=c(0,1))
-performance(pred,"auc")@y.values[[1]]
}

algo <- "Nelder-Mead"
algo <- "BFGS"
algo <- "L-BFGS-B"
algo <- "CG"
est = optim(c(2.9,0.03,0),maxauc,method=algo)
est = optim(c(1,0.05,0),maxauc,method=algo)
est = optim(c(0.6,0.03,0.9),maxauc,method=algo)
est$convergence
est$par
est$value



# ---------------------------------------------------------------------------------------------------------
# SVM sur coordonnées factorielles
# ---------------------------------------------------------------------------------------------------------

library(e1071)
library(kernlab)
library(ROCR)
library(FactoMineR)

summary(credit2)
credit4 <- credit2
credit4$Taux_effort <- NULL
credit4$Anciennete_domicile <- NULL
credit4$Nb_credits <- NULL
credit4$Nb_pers_charge <- as.factor(ifelse(credit4$Nb_pers_charge==1,"NbP1","NbP2"))
str(credit4)
head(credit4)
names(credit4)


# ACM
# -----------------------------------------------

ACM <- MCA(credit4, ncp = 32, axes = c(1,2), graph = TRUE, quali.sup = 14)
ACM$eig
ACM$var$coord

# nombre de modalités
modal = apply(credit4[,-14], 2, function(x) nlevels(as.factor(x)))
# nb de valeurs propres non triviales
# = nb modalités - nb variables
sum(modal)-ncol(credit4[,-14])
# inertie totale = (# modalités / # variables) - 1
(sum(modal)/ncol(credit4[,-14]))-1
# inertie totale = 2 = somme des valeurs propres sur ACM sur 
# le tableau disjonctif complet
sum(ACM$eig[,1])


# Test d'une solution (avec kernlab)
# -----------------------------------------------

naxes <- 2
x <- ACM$ind$coord[id,1:naxes]
y <- credit[id,"Cible"]
xt <- ACM$ind$coord[-id,1:naxes]
yt <- credit[-id,"Cible"]
s <- scale(x)

#ksvmlin <- ksvm(y~x,data=cbind(x,y),type="C-svc",kernel="vanilladot",C=1,scaled=c())
ksvmlin <- ksvm(s,y,type="C-bsvc",kernel="vanilladot",C=1,scaled=F,prob.model = TRUE)
#ksvmlin <- ksvm(cbind(x[,1],x[,2]),y,type="C-svc",kernel="vanilladot",probability=TRUE,C=1,prob.model = TRUE)
ksvmlin

# affichage pour 2 axes factoriels
plot(ksvmlin,data=cbind(x,y))
plot(ksvmlin,data=x)

# On extrait w and b du modèle
(w <- colSums(coef(ksvmlin)[[1]] * x[SVindex(ksvmlin),]))
#(w <- colSums(coef(ksvmlin)[[1]] * x[unlist(alphaindex(ksvmlin)),]))
(b <- b(ksvmlin))

# autre représentation graphique
(ys <-ymatrix(ksvmlin))
(xs <-as.vector(xmatrix(ksvmlin)))
head(xs)
as.numeric(y)
all.equal(ymatrix(ksvmlin),as.numeric(y))
# affichage des points supports en caractères plus petits
# les positifs (dossiers risqués) sont des cercles rouges
# les négatifs (dossiers non risqués) sont des + bleus
# à noter que les positifs sont tous des points supports
# cad mal classés ou dans la marge

# représentation avec Axe 2 en abscisse
plot(x[-SVindex(ksvmlin),2], x[-SVindex(ksvmlin),1], 
col = ifelse(as.numeric(y)[-SVindex(ksvmlin)]>1,"red","blue"), 
pch = ifelse(as.numeric(y)[-SVindex(ksvmlin)]>1, 1, 3),lwd=3, 
xlim=c(min(x[,2]), max(x[,2])), ylim=c(min(x[,1]), 
max(x[,1])),xlab=colnames(x)[2],ylab=colnames(x)[1])
points(x[SVindex(ksvmlin),2], x[SVindex(ksvmlin),1], 
col = ifelse(as.numeric(y)[SVindex(ksvmlin)]>1,"red","blue"), 
pch = ifelse(as.numeric(y)[SVindex(ksvmlin)]>1, 1, 3))
# tracé des frontières
-b/w[1]
-w[2]/w[1]
-w[1]/w[2]
abline(b/w[1],w[2]/w[1])
abline((b+1)/w[1],w[2]/w[1],lty=2)
abline((b-1)/w[1],w[2]/w[1],lty=2)

# représentation avec Axe 1 en abscisse
plot(x[-SVindex(ksvmlin),1], x[-SVindex(ksvmlin),2], 
col = ifelse(as.numeric(y)[-SVindex(ksvmlin)]>1,"red","blue"), 
pch = ifelse(as.numeric(y)[-SVindex(ksvmlin)]>1, 1, 3),lwd=3, 
xlim=c(min(x[,1]), max(x[,1])), ylim=c(min(x[,2]), 
max(x[,2])),xlab=colnames(x)[1],ylab=colnames(x)[2])
points(x[SVindex(ksvmlin),1], x[SVindex(ksvmlin),2], 
col = ifelse(as.numeric(y)[SVindex(ksvmlin)]>1,"red","blue"), 
pch = ifelse(as.numeric(y)[SVindex(ksvmlin)]>1, 1, 3))
# tracé des frontières
abline(b/w[2],w[1]/w[2])
abline(-(b+1)/w[2],w[1]/w[2],lty=2)
abline(-(b-1)/w[2],w[1]/w[2],lty=2)

# v2
plot(s[-SVindex(ksvmlin),1], s[-SVindex(ksvmlin),2], 
col = ifelse(as.numeric(y)[-SVindex(ksvmlin)]>1,"red","blue"), 
pch = ifelse(as.numeric(y)[-SVindex(ksvmlin)]>1, 1, 3),lwd=3, 
xlim=c(min(s[,1]), max(s[,1])), ylim=c(min(s[,2]), 
max(s[,2])),xlab=colnames(x)[1],ylab=colnames(x)[2])
points(s[SVindex(ksvmlin),1], s[SVindex(ksvmlin),2], 
col = ifelse(as.numeric(y)[SVindex(ksvmlin)]>1,"red","blue"), 
pch = ifelse(as.numeric(y)[SVindex(ksvmlin)]>1, 1, 3))
# coefficients du modèle
(w <- colSums(coef(ksvmlin)[[1]] * s[SVindex(ksvmlin),]))
(b <- b(ksvmlin))
(norme <- sqrt(w[1]^2+w[2]^2))
(norme2 <- sqrt(w[1]^2+w[2]^2+b^2))
# tracé des frontières
abline(b/norme2,norme/norme2)
abline((b+norme2)/norme2,norme/norme2,lty=2)
abline((b-norme2)/norme2,norme/norme2,lty=2)

# prédiction
library(ROCR)
#predksvm = predict(ksvmlin,type="prob",cbind(xt[,1],xt[,2]))
predksvm = predict(ksvmlin,type="prob",xt)
pred=prediction(predksvm[,2],yt,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]


# Coefficients de la solution (avec e1071)
# -----------------------------------------------

# détermination des coefficients des prédicteurs
svmlin = svm(x, y,kernel="linear",probability=TRUE,cost=1,scale=F)
# score calculé par e1071
p1 <- predict(svmlin, newdata=xt, decision.values=T)
p1 <- attr(p1, "decision.values")
head(p1)
pred=prediction(p1,yt,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

svmlin$index
SVindex(ksvmlin)

plot(x)
plot(x[-svmlin$index,1], x[-svmlin$index,2], 
col = ifelse(as.numeric(y[-svmlin$index])>1,"red","blue"), 
pch = ifelse(as.numeric(y[-svmlin$index])>1, 1, 3),lwd=3, 
xlim=c(min(x[,1]), max(x[,1])), ylim=c(min(x[,2]), 
max(x[,2])),xlab=colnames(x)[1],ylab=colnames(x)[2])
points(x[svmlin$index,1], x[svmlin$index,2], 
col = ifelse(as.numeric(y[svmlin$index])>1,"red","blue"), 
pch = ifelse(as.numeric(y[svmlin$index])>1, 1, 3))

w <- t(svmlin$coefs) %*% x[svmlin$index,]
w <- t(svmlin$coefs) %*% scale(x[svmlin$index,],apply(x,2,mean),apply(x,2,sd))
# calcul du score pour données non centrées-réduites
p2 <- xt %*% t(w) - svmlin$rho
# calcul du score pour données centrées-réduites
xscaled <- scale(xt,svmlin$x.scale[[1]],svmlin$x.scale[[2]])
p2 <- xscaled %*% t(w) - svmlin$rho
svmlin$rho
# coefficients du score SVM sur prédicteurs initiaux (avant standardisation)
# à la constante près
w/svmlin$x.scale[[2]]
head(p2)
# comparaison des scores
cor(p1,p2,method="spearman")
max(abs(p1-p2))

# passage des coefficients sur coordonnées factorielles
# aux coefficients sur prédicteurs initiaux
ACM$var$coord[,1:naxes]
t(ACM$var$coord[,1:naxes])
# coefficient de CC [0-200 euros[
(0.4962056*0.15717369)+(1.658412*0.42565536)+(0.4342413*0.08898376)+(0.5120888*0.55367832)-(0.2799403*0.05418345)-(0.6950526*0.71052167)
coef <- w %*% t(ACM$var$coord[,1:naxes])
coef

# utilisation de la fonction "model.matrix" avec une option de contraste
# pour avoir tous les niveaux (y compris celui de référence)
X <- model.matrix( ~ .-1,data=credit4[,-which(names(credit4)=="Cible")],
contrasts.arg = lapply(credit4[,-which(names(credit4)=="Cible")], contrasts, contrasts=FALSE))
head(X,1)
# utilisation d'un autre package pour générer les indicatrices avec les noms
# des modalités sans les noms de leurs variables, comme FactoMineR
# (afin de permettre le rapprochement)
library(caret)
X <- dummyVars( ~ ., data=credit4[-id,-which(names(credit4)=="Cible")],levelsOnly = T)
X <- predict(X,newdata=credit4[-id,-which(names(credit4)=="Cible")])
p3 <- X %*% t(coef)
head(p3)
# comparaison des scores
cor(p1,p3,method="pearson")
cor(p1,p3,method="spearman")
plot(p1,p3)
pred=prediction(p3,yt,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]


# variante : ACM avec package MASS
# -----------------------------------------------

library(MASS)
str(credit4)
ACM2 <- mca(credit4[,-14],nf=32)
ACM2$d # valeurs singulières
ACM2$d^2 # valeurs propres

# inertie totale = 2 = somme des valeurs propres sur ACM sur 
# le tableau disjonctif complet
sum(ACM2$d^2)

ACM2$cs[,1:3]
ACM2$rs[1:5,1:5]

naxes <- 6
# on multiplie les coordonnes factorielles des individus par p*sqrt(n)
# pour se ramener aux coordonnées factorielles usuelles
x <- ACM2$rs[id,1:naxes]*(ncol(credit4)-1)*sqrt(nrow(credit4))
x <- ACM2$rs[id,1:naxes]
y <- credit[id,"Cible"]
xt <- ACM2$rs[-id,1:naxes]*(ncol(credit4)-1)*sqrt(nrow(credit4))
xt <- ACM2$rs[-id,1:naxes]
yt <- credit[-id,"Cible"]
# on multiplie les coordonnes factorielles des variables par p*sqrt(n)*sqrt(eigenvalue)
# pour se ramener aux coordonnées factorielles usuelles
ACM2$cs <- ACM2$cs*(ncol(credit4)-1)*sqrt(nrow(credit4))*ACM2$d
# en effet : coord_var(MASS) = coord_var(FactoMineR) / (p * sqrt(n) * sqrt(eigenvalue)) 

# détermination des coefficients des prédicteurs
svmlin = svm(x, y,kernel="linear",probability=TRUE,cost=1,scale=F)
# score calculé par e1071
p1 <- predict(svmlin, newdata=xt, decision.values=T)
p1 <- attr(p1, "decision.values")
head(p1)
pred=prediction(p1,yt,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]
# coefficients si prédicteurs non centrés-réduits
w <- t(svmlin$coefs) %*% x[svmlin$index,]
# coefficients si prédicteurs centrés-réduits
w <- t(svmlin$coefs) %*% scale(x[svmlin$index,],apply(x,2,mean),apply(x,2,sd))
# calcul du score pour données non centrées-réduites
p2 <- xt %*% t(w) - svmlin$rho
head(p2)
# calcul du score pour données centrées-réduites
xscaled <- scale(xt,apply(x,2,mean),apply(x,2,sd))
p2 <- xscaled %*% t(w) - svmlin$rho
# coefficients du score SVM sur prédicteurs initiaux (avant standardisation)
# à la constante près
w/apply(x,2,sd)

# comparaison des scores
cor(p1,p2,method="spearman")
plot(p1,p2)

# passage des coefficients sur coordonnées factorielles
# aux coefficients sur prédicteurs initiaux
coef <- (w/apply(x,2,sd)) %*% t(ACM2$cs[,1:naxes])
coef <- w %*% t(scale(ACM2$cs[,1:naxes],apply(x,2,mean),apply(x,2,sd)))
# facteurs non centrés-réduits
coef <- w %*% t(ACM2$cs[,1:naxes])
# utilisation d'un autre package pour générer les indicatrices avec les noms
# des modalités sans les noms de leurs variables, comme FactoMineR
# (afin de permettre le rapprochement)
library(caret)
X <- dummyVars( ~ ., data=credit4[-id,-which(names(credit4)=="Cible")],levelsOnly = F)
X <- predict(X,newdata=credit4[-id,-which(names(credit4)=="Cible")])
p3 <- X %*% t(coef)
head(p3)
# comparaison des scores
cor(p1,p3,method="pearson")
cor(p1,p3,method="spearman")
plot(p1,p3)
pred=prediction(p3,yt,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# construction de la grille de score
var <- as.vector(unlist(strsplit(colnames(coef),".",fixed=T)))
param <- data.frame(VARIABLE=var[seq(1,2*length(coef),by=2)],MODALITE=var [seq(2,2*length(coef),by=2)],COEF=as.numeric(coef))
param

# calcul du poids total pour normalisation
mini=aggregate(data.frame(min = param$COEF), by = list(VARIABLE = param$VARIABLE), min)
maxi=aggregate(data.frame(max = param$COEF), by = list(VARIABLE = param$VARIABLE), max)
total=merge(mini,maxi)
total$diff = total$max - total$min
poids_total = sum(total$diff)
# calcul des poids par modalité
grille = merge(param,mini,all.x=TRUE)
grille$delta = grille$COEF - grille$min
grille$POIDS = round((100*grille$delta) / poids_total)
#grille[which(grille$VARIABLE!=""),c("VARIABLE","MODALITE","COEF","POIDS")]
grille[order(grille$VARIABLE,grille$MODALITE),c("VARIABLE","MODALITE","POIDS")]


# Recherche d'une solution optimale en termes d'AUC
# -----------------------------------------------

# noyau linéaire avec package kernlab
f <- function(i)
{
ksvmlin <- ksvm(x,y,type="C-svc",kernel="vanilladot",C=i,prob.model = TRUE,verbose=F)
predksvm = predict(ksvmlin,type="prob",xt)
pred=prediction(predksvm[,2],yt,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]
}
f(1)

# noyau linéaire avec package e1071
library(e1071)
f <- function(i)
{
ksvm <- svm(x,y,kernel="linear",cost=i,probability=TRUE)
predksvm = predict(ksvm,type="prob",xt,probability=TRUE)
pred=prediction(attr(predksvm,"probabilities")[,1],yt,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]
}
f(0.33)

# noyau radial avec package e1071
library(e1071)
f <- function(i,j)
{
ksvm <- svm(x,y,kernel="radial",cost=i,gamma=j,probability=TRUE)
predksvm = predict(ksvm,xt,probability=TRUE)
pred=prediction(attr(predksvm,"probabilities")[,1],yt,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]
}
f(0.05,0.01)

# noyau polynomial avec package e1071
library(e1071)
f <- function(i,j)
{
ksvm <- svm(x,y,kernel="polynomial",degree=2,cost=i,gamma=j,probability=TRUE)
predksvm = predict(ksvm,type="prob",xt,probability=TRUE)
pred=prediction(attr(predksvm,"probabilities")[,1],yt,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]
}
f(3,0.03)


# utilisation de la fonction f précédente pour trouver
# le nb d'axes factoriels et le coût maximisant l'AUC en test
# noyau linéaire
for (n in seq(2,12,by=1))
{
i <- seq(0.01,2,by=0.01)
x <- ACM$ind$coord[id,1:n]
xt <- ACM$ind$coord[-id,1:n]
k <- Vectorize(f)(i)
cout <- i[which(k==max(k),arr.ind=TRUE)[1]]
cat("\n","nb facteurs = ",n," AUC test max = ",max(k)," coût = ",
cout)
}

# utilisation de la fonction f précédente pour trouver
# le nb d'axes factoriels et le coût maximisant l'AUC en test
# noyau radial ou polynomial
for (n in seq(2,12,by=1))
{
i <- seq(0.1,5,by=0.1)
j <- seq(0.01,0.5,by=0.01)
x <- ACM$ind$coord[id,1:n]
xt <- ACM$ind$coord[-id,1:n]
k <- outer(i,j,Vectorize(f))
cout <- i[which(k==max(k),arr.ind=TRUE)[1,1]]
gamma <- j[which(k==max(k),arr.ind=TRUE)[1,2]]
cat("\n","nb facteurs = ",n," AUC test max = ",max(k)," coût = ",
cout," gamma = ",gamma)
}


# modèle logit sur coordonnées factorielles (pour comparaison avec SVM)
# -----------------------------------------------

naxes <- 6
x <- ACM$ind$coord[id,1:naxes]
y <- credit[id,"Cible"]
xt <- ACM$ind$coord[-id,1:naxes]
yt <- credit[-id,"Cible"]
library(ROCR)

# modèle logit sur coordonnées factorielles (pour comparaison avec SVM)
acm <- as.data.frame(cbind(x,y))
acm$y <- acm$y - 1
#head(acm)
logit <- glm(y~., data=acm,family=binomial(link = "logit"))
summary(logit)
predlogit <- predict(logit, newdata=as.data.frame(xt), type="response")
pred=prediction(predlogit,yt,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]
# le modèle logit sur 6 axes factoriels est légèrement plus discriminant
# que le modèle SVM à noyau linéaire

# régression logistique avec sélection pas à pas
# il faut renommer les variables à cause des " " dans leurs noms
names(acm)
colnames(acm) <- c(paste("Dim",1:naxes,sep=""),'y')
(predicteurs <- -grep('y', names(acm)))
(formule <- as.formula(paste("y ~ ",paste(names(acm[,predicteurs]),collapse="+"))))

# sélection ascendante sur les axes factoriels
modcst <- glm(y~1,data=acm,family=binomial(link = "logit"))
summary(modcst)
(formule <- as.formula(paste("y ~ ",paste("Dim",1:6,sep="",collapse="+"))))
selection <- glm(formule, data=acm,family=binomial(link = "logit"))
selection <- step(modcst,direction="forward",trace=TRUE,k = log(nrow(acm)),
selection <- step(modcst,direction="forward",trace=TRUE,k = 2,
scope=list(upper=formule))
summary(selection)
xtest <- as.data.frame(xt)
names(xtest)
colnames(xtest) <- paste("Dim",1:naxes,sep="")
predlogit <- predict(selection, newdata=xtest, type="response")
pred=prediction(predlogit,yt,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]

# passage des coefficients sur coordonnées factorielles
# aux coefficients sur prédicteurs initiaux
naxes <- 6
t(ACM$var$coord[,1:naxes])
logit$coefficients[1:naxes+1]
(coef <- logit$coefficients[1:naxes+1] %*% t(ACM$var$coord[,1:naxes]))

# utilisation d'un autre package pour générer les indicatrices avec les noms
# des modalités sans les noms de leurs variables, comme FactoMineR
# (afin de permettre le rapprochement)
summary(logit)
library(caret)
X <- dummyVars( ~ ., data=credit4[-id,-which(names(credit4)=="Cible")],levelsOnly = T)
X <- predict(X,newdata=credit4[-id,-which(names(credit4)=="Cible")])
p3 <- X %*% t(coef)
head(p3)
p1 <- predict(logit, newdata=as.data.frame(xt), type="response")
# comparaison des scores
cor(p1,p3,method="pearson")
cor(p1,p3,method="spearman")
plot(p1,p3)
pred=prediction(p3,yt,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]



# ---------------------------------------------------------------------------------------------------------
# CHAPITRE 16
# Réseaux de neurones
# ---------------------------------------------------------------------------------------------------------

# On charge le package nnet pour les réseaux de neurones (perceptron à une couche cachée)
library(nnet)
# paramètres importants :
# size : nb noeuds dans couche cachée
# decay : paramètre de régularisation analogue à celui utilisé en régression ridge
# decay pénalise la norme du vecteurs des paramètres et contraint ainsi la flexibilité du modèle

library(pROC)
library(ROCR)

# valeurs possibles de size :
(ncol(train)-1)*2
2*sqrt((ncol(train)-1)*2)

# données brutes
train  <- credit[id,]
train$Cle <- NULL
valid  <- credit[-id,]
valid$Cle <- NULL
# données travaillées
train  <- credit2[id,]
valid  <- credit2[-id,]
summary(train)

# réseau avec sortie sigmoïde
seed <- 235
set.seed(seed)
rn <- nnet(Cible~., data=train, size=10, decay=2.2, maxit=200, softmax=F)
#summary(rn)
# application du réseau
valid$rn <- predict(rn, newdata = valid, type = "raw")
head(valid$rn)
auc(valid$Cible,valid$rn)

# réseau avec sortie softmax
seed <- 235
set.seed(seed)
rn <- nnet(class.ind(Cible)~., data=train, size=10, decay=2.2, maxit=200, softmax=T)
#summary(rn)
# application du réseau
valid$rn <- predict(rn, newdata = valid, type = "raw")[,2]
head(valid$rn)
auc(valid$Cible,valid$rn)
# l'AUC est parfois plus basse avec softmax


# optimisation des paramètres size et decay
library(e1071)
# optimisation sur échantillon de test
names(valid)
valid$rn <- NULL
optirn <- tune.nnet(Cible ~ ., data=train,size=1:20,decay=seq(0.1,5,by=0.1),
validation.x = valid [,-17], validation.y = valid [,17],
tunecontrol = tune.control(sampling = "fix", fix=1))
# optimisation par validation croisée
optirn <- tune.nnet(Cible~., data=train,size=1:10,decay=0:10)
# paramètres optimaux
optirn
plot(optirn)


# calcul de l'AUC en test pour chaque combinaison (size,decay)
f <- function(i,j)
{
set.seed(seed)
#cat("\n","size = ",i," decay = ",j," ")
rn <- nnet(Cible~., data=train, size=i, decay=j, maxit=100, trace=FALSE)
auc(valid$Cible,predict(rn, newdata = valid, type = "raw"))
#pred <- prediction(predict(rn, newdata = valid, type = "raw"),valid$Cible,label.ordering=c(0,1))
#performance(pred,"auc")@y.values[[1]]
}
f(2,1)
f(1,1.5)

# calcul de l'AUC en test pour chaque combinaison (size,decay) - variante softmax
f <- function(i,j)
{
set.seed(seed)
rn <- nnet(class.ind(Cible)~., data=train, size=i, decay=j, maxit=200, softmax=T, trace=FALSE)
auc(valid$Cible,predict(rn, newdata = valid, type = "raw")[,2])
#pred <- prediction(predict(rn, newdata = valid, type = "raw"),valid$Cible,label.ordering=c(0,1))
#performance(pred,"auc")@y.values[[1]]
}
f(2,0.1)

# utilisation de la fonction f précédente
i <- seq(1,20,by=1)
j <- seq(0.1,5,by=0.1)
ptm <- proc.time()
k <- outer(i,j,Vectorize(f))
proc.time() - ptm
max(k)
which.max(k)
which(k==max(k),arr.ind=TRUE)
(size <- i[which(k==max(k),arr.ind=TRUE)[1]])
(decay <- j[which(k==max(k),arr.ind=TRUE)[2]])
which(k==min(k),arr.ind=TRUE)

persp(i,j,k)
library(rgl)
persp3d(i,j,k)
contour(i,j,k)
# heatmap
image(i,j,k, col=gray((0:32)/32))
image(i,j,k, col=gray((64:18)/64))
box() # pour entourer le graphique



# ---------------------------------------------------------------------------------------------------------
# Boosting de réseaux de neurones
# ---------------------------------------------------------------------------------------------------------

library(nnet)
library(ROCR)

# discrete adaboost - avec correctif pour erreurs > 0,5
nnet.adaboost = function(apprent, validat, M, seed, formule, noeuds=2, penal=1)
{
  n <- nrow(apprent)
  w <- rep(1/n, n)
  rep <- 0    
  auc <- matrix(0,M,1)
  
  for(i in 1:M)
  {
    w <- w/sum(w)
    set.seed(seed)
    s <- sort(sample(n,n,replace=T,prob=w))
    rn <- nnet(formule, data = apprent[s,], size=noeuds, decay=penal, softmax = TRUE, maxit=100, trace=F)
    pred <- predict(rn, apprent, type = "class")
    badPred <- 1 - as.numeric(pred == apprent$Cible)
    tauxErr <- sum(badPred * w) / sum(w)
    if(tauxErr < 0.5 & tauxErr != 0)
    {
      alpha <- log((1-tauxErr) / tauxErr)
      w <- w * exp(alpha*badPred)
    }
    if(tauxErr >= 0.5)
    {
      alpha <- 0
      w <- rep(1/n, n)
    }
    rep <- rep + (alpha * predict(rn, newdata = validat, type = "raw")[,2])
    #rep <- rep + (alpha * predict(rn, newdata = validat, type = "raw"))
    pred <- prediction(rep,validat$Cible,label.ordering=c(0,1))
    auc[i] <- performance(pred,"auc")@y.values[[1]]
  }
  
  return(list(auc,rep))
}

# discrete adaboost
discrete.adaboost = function(apprent, validat, M, seed, formule, noeuds=2, penal=1)
{
  n <- nrow(apprent)
  w <- rep(1/n, n)
  rep <- 0    
  auc <- matrix(0,M,1)
  
  for(i in 1:M)
  {
    w <- w/sum(w)
    set.seed(seed)
    s <- sort(sample(n,n,replace=T,prob=w))
    rn <- nnet(formule, data = apprent[s,], size=noeuds, decay=penal, softmax = TRUE, maxit=100, trace=F)
    pred <- predict(rn, apprent, type = "class")
    badPred <- 1 - as.numeric(pred == apprent$Cible)
    tauxErr <- sum(badPred * w) / sum(w)
    if(tauxErr != 0)
    {
      alpha <- log((1-tauxErr) / tauxErr)
      w <- w * exp(alpha*badPred)
    }
    rep <- rep + (alpha * predict(rn, newdata = validat, type = "raw")[,2])
    #rep <- rep + (alpha * predict(rn, newdata = validat, type = "raw"))
    pred <- prediction(rep,validat$Cible,label.ordering=c(0,1))
    auc[i] <- performance(pred,"auc")@y.values[[1]]
  }
  
  rep <- rep/M
  return(list(auc,rep))
}

# real adaboost
real.adaboost = function(apprent, validat, M, seed, formule, noeuds=2, penal=1)
{
  n <- nrow(apprent)
  w <- rep(1/n, n)
  rep <- 0    
  auc <- matrix(0,M,1)
  
  for(i in 1:M)
  {
    w <- w/sum(w,na.rm=TRUE) # gérer le cas où un poids w est si petit qu'il vaut NaN
    set.seed(seed)
    s <- sort(sample(n,n,replace=T,prob=w))
    rn <- nnet(formule, data = apprent[s,], size=noeuds, decay=penal, softmax = TRUE, maxit=100, trace=F)
    pred <- min(predict(rn, apprent, type = "raw")[,2],0.999999)
    alpha <- log(pred/(1-pred))/2
    signe = ifelse(apprent$Cible == 1,1,-1)
    w <- w * exp(-signe*alpha)
    rep <- rep + predict(rn, newdata = validat, type = "raw")[,2]
    auc[i] <- performance(prediction(rep,validat$Cible,label.ordering=c(0,1)),"auc")@y.values[[1]]
  }
  
  rep <- rep/M
  return(list(auc,rep))
}

valid$Cle <- NULL

f <- function(i,j)
{
seed <- 235
formule = class.ind(Cible)~.
adaboost <- real.adaboost(apprent=train, validat=valid, M=3000, seed=seed, formule = formule, 
noeuds=i, penal=j)
# résultats en retour
return(list(max(adaboost[[1]]),which.max(adaboost[[1]])))
}
f(15,0.01)

# utilisation de la fonction f précédente
i <- c(2,5,10,15,20)	
j <- c(1,0.1,0.01,0.001)
k <- outer(i,j,Vectorize(f))

seed <- 235
formule = class.ind(Cible)~.
adaboost <- discrete.adaboost(apprent=train, validat=valid, M=5000, seed=seed, formule = formule, 
noeuds=2, penal=1)
summary(adaboost[[2]])
max(adaboost[[1]])
which.max(adaboost[[1]])



# ---------------------------------------------------------------------------------------------------------
# Forêts aléatoires de réseaux de neurones
# ---------------------------------------------------------------------------------------------------------

#install.packages('gplots')
library(nnet)
library(ROCR)

table(credit$Objet_credit)
credit$Objet_credit[credit$Objet_credit == "A410"] <- "A48"
table(credit$Objet_credit)
train  <- credit[id,]
valid  <- credit[-id,]
summary(valid)
str(train)

library(car)
credit2$Objet_credit <- recode(credit$Objet_credit,"'A40'='Voiture neuve';
'A41'='Voiture occasion';c('A42','A44','A45')='Intérieur';
'A43'='Vidéo - HIFI';c('A46','A48','A410')='Etudes';
'A47'='Vacances';'A49'='Business';else='Autres'")
train  <- credit2[id,]
valid  <- credit2[-id,]

varX <- c(varquali, varquanti)
varY <- "Cible"
nsimul <- 10
nb_var <- 2
nobs <- nrow(train)
y <- train[,varY]
predic   <- matrix(NA,nrow(valid),nsimul)
auc      <- matrix(0,nsimul,1)
varsel   <- matrix(0,nsimul,nb_var)
sv <- sample(length(varX), nb_var, replace=F)
si <- sample(nobs,nobs,replace=T)
i <- 1
(varsel[i,] <- sv)
#sv <- c(1,2,17)
cible <- train[si,varY]
varex <- train[si,varX[sv]]
head(train[si,varx[sv]])
head(train[,varX[sv]])
train1 <- as.data.frame(train[si,varX[sv]])
colnames(train1)=varX[sv]
names(train1)
head(train1)
rn <- nnet(cible~., data = train1, size=2, decay=1, softmax = F, maxit=100, trace=F)
rn <- nnet(cible~., data = train[,c(varY,varX[sv])], size=2, decay=1, softmax = F, maxit=100, trace=F)

rn <- nnet(cible~varex, data = train[si,c(varY,varX[sv])], size=2, decay=1, softmax = F, maxit=100, trace=F)
rn <- nnet(class.ind(cible)~., data = train[si,varx[sv]], size=2, decay=1, softmax = T, maxit=100, trace=F)
summary(rn)
predict(rn, newdata = valid, type = "raw")
predic[,i] <- predict(rn, newdata = valid, type = "raw")
#predic[,i] <- predict(rn, newdata = valid, type = "raw")[,2] # softmax
head(predic[,1])
if (i == 1) {moypred <- predic[,i]} else {moypred <- apply(predic[,1:i],1,mean)}
moypred
# calcul de l’AUC des i premiers modèles agrégés
cible <- valid[,varY]
pred <- prediction(moypred,cible,label.ordering=c(0,1))
(auc[i] <- performance(pred,"auc")@y.values[[1]])


nnet.forest = function(apprent, validat, varY, varX, nb_var, nsimul, seed, noeuds=2, penal=1)
{

nobs <- nrow(apprent)
y <- apprent[,varY]
predic   <- matrix(NA,nrow(validat),nsimul)
auc      <- matrix(0,nsimul,1)
varsel   <- matrix(0,nsimul,nb_var)

for(i in (1:nsimul))
{

# tirage aléatoire simple des variables
sv <- sort(sample(length(varX), nb_var, replace=F))
varsel[i,] <- sv

# tirage aléatoire avec remise de l'échantillon d'apprentissage
si <- sort(sample(nobs,nobs,replace=T))

# réseau de neurones
cible <- apprent[si,varY]
explvar <- as.data.frame(apprent[si,varX[sv]])
colnames(explvar) <- varX[sv]
#rn <- nnet(cible~., data = apprent[si,varX[sv]], size=noeuds, decay=penal, softmax = F, maxit=200, trace=F)
# variante : tirage aléatoire des paramètres
#noeuds <- sample(seq(1,20,1),1)
#penal  <- sample(c(1,0.1,0.01,0.001),1)
#rn <- nnet(cible~., data = explvar, size=noeuds, decay=penal, softmax = F, maxit=100, trace=F)
# variante : randomisation de la largeur de la plage des poids initiaux
#plage <- runif(1,min=0.1,max=0.9)
plage <- 1
#plage <- rnorm(1,0.5,0.17)
rn <- nnet(cible~., data = explvar, size=noeuds, decay=penal, softmax = F, maxit=100, trace=F, rang=plage)

# application du modèle à l'échantillon de test
predic[,i] <- predict(rn, newdata = validat, type = "raw")
# calcul de l’AUC des i premiers modèles agrégés
if (i == 1) {moypred <- predic[,i]} else {moypred <- rowMeans(predic[,1:i])}
cible <- validat[,varY]
pred <- prediction(moypred,cible,label.ordering=c(0,1))
auc[i] <- performance(pred,"auc")@y.values[[1]]

} # fin de la boucle

# résultats en retour
rf <- list(auc,predic,varsel)

} # fin de la fonction

niter <- 500
# appel de la fonction
#(varx <- c(varquali, varquanti))
ptm <- proc.time()
rf <- nnet.forest(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=3, nsimul=niter, seed=seed, noeuds=2, penal=0)
proc.time() - ptm

# AUC
which.max(rf[[1]])
max(rf[[1]])
tail(rf[[1]])
plot(rf[[1]],xlab = '# agrégations', ylab = 'AUC')
# prédiction moyenne
rowMeans(rf[[2]][,1:100])
# AUC
pred <- prediction(rowMeans(rf[[2]][,1:100]),valid$Cible,label.ordering=c(0,1))
pred <- prediction(rf[[2]][,19],valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]]
# AUC de chacun des modèles neuronaux agrégés
roc <- function(i) { 
pred <- prediction(rf[[2]][,i],valid$Cible,label.ordering=c(0,1))
performance(pred,"auc")@y.values[[1]] }
auctest <- Vectorize(roc)(1:niter)
length(which(auctest>0.7))
length(which(auctest>0.75))
barplot(Vectorize(roc)(1:niter),names.arg=seq(1,niter),cex.names = 0.8,las=3,ylim=c(0,0.8),ylab='AUC en test')
abline(h=c(0.7,0.75),lty=2)
# liste des variables pour chacun des modèles agrégés
rf[[3]]
# liste des variables donnant un modèle avec AUC > 0.7
rf[[3]][which(auctest>0.6),]
table(rf[[3]][which(auctest>0.6),])

# balayage des valeurs du nb d'unités cachées et du weight decay
f <- function(i,j)
{
rf <- nnet.forest(apprent=train, validat=valid, varY="Cible", varX=varx, nb_var=nbpred, nsimul=niter, noeuds=i, penal=j)
# résultats en retour
#return(list(max(rf[[1]]),which.max(rf[[1]]),rf[[1]]))
return(max(rf[[1]]))
}
rfnn <- f(1,0)
rfnn[[1]]

# utilisation de la fonction f précédente
niter <- 500
i <- c(2,5,10,15,20)	
j <- c(1,0.1,0.01,0.001,0)
i <- c(5,10,15)	
j <- c(0.1,0.01)
nbpred <- 1
k <- outer(i,j,Vectorize(f))
k
nbpred <- 2
k <- outer(i,j,Vectorize(f))
k

persp(i,-j,k)
library(rgl)
persp3d(i,-j,k)



# ---------------------------------------------------------------------------------------------------------
# ANNEXE 1
# Détection des règles d'association
# ---------------------------------------------------------------------------------------------------------

library(arules)

# transformation de l'ensemble des variables qualitatives en facteurs
vardiscr <- c("Cible","Comptes", "Epargne", "Historique_credit", "Objet_credit",
 "Situation_familiale", "Garanties", "Biens", "Autres_credits",
 "Statut_domicile", "Type_emploi", "Anciennete_emploi", "Telephone",
 "Nb_pers_charge","Taux_effort","Anciennete_domicile","Nb_credits")
indices <- names(credit) %in% vardiscr
for (i in (1:ncol(credit))) { if (indices[i] == 1) { credit[,i] <- factor(credit[,i]) } }

# discrétisation des variables continues
Age <- cut(credit$Age,c(0,25,Inf),right=TRUE) #intervalle semi-fermé à droite
Duree_credit <- cut(credit$Duree_credit,c(0,15,36,Inf),right=TRUE)
Montant_credit <- cut(credit$Montant_credit,c(0,4000,Inf),right=TRUE)

npred <- -grep('(Cle|Duree_credit|Montant_credit|Age)', names(credit))
credit2 <- credit[,npred]
credit2$Age <- Age
credit2$Duree_credit <- Duree_credit
credit2$Montant_credit <- Montant_credit

rules <- apriori(credit2,parameter = list(supp = 0.1, conf = 0.8,maxlen = 9, target = "rules"))
rules
rules <- apriori(credit2,parameter = list(supp = 0.06, conf = 0.6,maxlen = 5,
target = "rules"),appearance = list(rhs = "Cible=2",default="lhs"))

rules <- apriori(credit2,parameter = list(supp = 0.06, conf = 0.6,maxlen = 5,
target = "rules"),appearance = list(rhs = c("Cible=1","Cible=2"),default="lhs"))

summary(rules)
inspect(rules[1:10])
inspect(sort(rules, by = "confidence", decreasing = TRUE)[1:30])
inspect(sort(sort(rules, by = "support", decreasing = TRUE), by = "confidence", decreasing = TRUE)[1:5])
inspect(sort(sort(rules, by = "confidence", decreasing = TRUE), by = "lift", decreasing = TRUE)[1:30])
inspect(sort(sort(rules, by = "lift", decreasing = TRUE), by = "confidence", decreasing = TRUE)[1:30])

# graphe
#install.packages("arulesViz")
library(arulesViz)
plot(rules)
inspect(sort(rules, by = "lift", decreasing = TRUE)[1:10])
plot(rules,interactive=TRUE)
plot(rules, method="grouped", shading="confidence",interactive=F)
plot(rules, method="grouped",interactive=T)
plot(rules, shading="order", control=list(main = "Two-key plot"))
plot(rules, method="matrix3D", measure=c("lift", "confidence"))
plot(rules, method="graph") # calcul très long => il ne faut garder que qq dizaines de règles à afficher
plot(rules[quality(rules)$confidence > 0.999], method = "graph")



# ---------------------------------------------------------------------------------------------------------
# ANNEXE 2
# Régression logistique avec les packages speedglm, biglm et ff
# ---------------------------------------------------------------------------------------------------------

logit <- glm(Cible~Comptes+Historique_credit+Duree_credit+Age+Epargne+Garanties+Autres_credits,
data=train,family=binomial('logit'))
summary(logit)
predglm <- predict(logit, newdata=valid, type="response")
head(predglm)

# modèle logit avec speedglm
#install.packages("speedglm")
library(speedglm)
logits <- speedglm(Cible~Comptes+Historique_credit+Duree_credit+Age+Epargne+Garanties+Autres_credits,
data=train,family=binomial('logit'))
summary(logits)
w <- as.matrix(coef(logits))
library(caret)
xt <- dummyVars(Cible~Comptes+Historique_credit+Duree_credit+Age+Epargne+Garanties+Autres_credits,
data=valid,levelsOnly = F,sep = NULL,fullRank = T)
xt <- predict(xt,newdata=valid)
p2 <- (xt %*% w[-1]) + w[1]
predspeed <- 1/(1+exp(-p2))
cor(predglm,predspeed)
all.equal(predglm,as.numeric(predspeed),check.attributes = F)
identical(predglm,predspeed)

# modèle logit avec biglm
library(biglm)
logitb <- bigglm(as.numeric(Cible)-1~Comptes+Historique_credit+Duree_credit+Age+Epargne+Garanties+Autres_credits,
data=train,family=binomial('logit'),chunksize=10000)
summary(logitb)
predbig <- predict(logitb, newdata=valid, type="response")
all.equal(predbig,predspeed)

# comparaison des temps de calcul
install.packages('microbenchmark')
library(microbenchmark)
microbenchmark(glm(Cible~Comptes+Historique_credit+Duree_credit+Age+Epargne+Garanties+Autres_credits,
data=train,family=binomial('logit')), times = 100)
microbenchmark(speedglm(Cible~Comptes+Historique_credit+Duree_credit+Age+Epargne+Garanties+Autres_credits,
data=train,family=binomial('logit')), times = 100)
microbenchmark(bigglm(as.numeric(Cible)-1~Comptes+Historique_credit+Duree_credit+Age+Epargne+Garanties+Autres_credits,
data=train,family=binomial('logit'),chunksize=10000), times = 100)

# régression logistique avec le package ff
library(ff)
library(ffbase)
library(biglm)
trainff <- as.ffdf(train)
validff <- as.ffdf(valid)
logitf <- bigglm(as.numeric(Cible)-1~Comptes+Historique_credit+Duree_credit+Age+Epargne+Garanties+Autres_credits,data=trainff)
summary(logitf)
object.size(logitf)
predff <- predict(logitf, newdata=validff, type="response")
predff <- 1/(1+exp(-predff))
all.equal(predbig,predff)
cor(predbig,predff)
microbenchmark(bigglm(as.numeric(Cible)-1~Comptes+Historique_credit+Duree_credit+Age+Epargne+Garanties+Autres_credits,
data=trainff,family=binomial('logit'),chunksize=1000000), times = 100)