ghoulaymen · January 17, 2017 21:22
diff --git a/MODELISATION PREDICTIVE ET APPRENTISSAGE STATISTIQUE AVEC R b/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([email protected][[1]][email protected][[1]])
 ks <- [email protected][[1]][email protected][[1]]
 (seuil <- pred@cutoffs[[1]][which.max(ks)])
 segments(seuil,[email protected][[1]][which.max(ks)],seuil,[email protected][[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 derreur (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 derreur (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 lAUC 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 dune 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 dune 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 dune 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 dune 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
 [email protected]
 # 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 lAUC 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 lAUC 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)